Invert the phase of diploid genotypes with +setGT -t a -n i

Author: rwk-unil

plugins/setGT.c

@@ -82,6 +82,7 @@ args_t *args = NULL;
 #define GT_CUSTOM   (1<<10)
 #define GT_X_VAF    (1<<11)
 #define GT_RAND     (1<<12)
+#define GT_INVPHASE (1<<13)
 
 #define MINOR_ALLELE   -1
 #define MAJOR_ALLELE   -2
@@ -113,6 +114,7 @@ const char *usage(void)
         "           X       .. allele with bigger read depth as determined from FMT/AD\n"
         "           p       .. phase genotype (0/1 becomes 0|1)\n"
         "           u       .. unphase genotype and sort by allele (1|0 becomes 0/1)\n"
+        "           i       .. invert the genotype phase (0|1 becomes 1|0)\n"
         "Usage: bcftools +setGT [General Options] -- [Plugin Options]\n"
         "Options:\n"
         "   run \"bcftools plugin\" for a list of common options\n"
@@ -235,6 +237,7 @@ int init(int argc, char **argv, bcf_hdr_t *in, bcf_hdr_t *out)
                 if ( strchr(optarg,'X') ) args->new_mask |= GT_X_VAF;
                 if ( strchr(optarg,'p') ) args->new_mask |= GT_PHASED;
                 if ( strchr(optarg,'u') ) args->new_mask |= GT_UNPHASED;
+                if ( strchr(optarg,'i') ) args->new_mask |= GT_INVPHASE;
                 if ( !strncmp(optarg,"c:",2) ) { args->new_mask |= GT_CUSTOM; args->custom.gt_str = optarg; }
                 if ( args->new_mask==0 ) error("Unknown parameter to --new-gt: %s\n", optarg);
                 break;
@@ -364,6 +367,21 @@ static inline int unphase_gt(int32_t *ptr, int ngts)
     return changed;
 }
 
+// invert the phase for a single sample, ngts is the ploidy
+static inline int invert_phase_gt(int32_t *ptr, int ngts)
+{
+    int gt;
+    if ( ngts!=2 ) return 0;    // don't invert the phase for sample ploidy != 2
+    if ( ptr[0]==bcf_int32_vector_end || ptr[1]==bcf_int32_vector_end ) return 0;
+    gt = bcf_gt_allele(ptr[1]);
+    if (bcf_gt_is_phased(ptr[1]))
+        ptr[1] = bcf_gt_phased(bcf_gt_allele(ptr[0]));
+    else
+        ptr[1] = bcf_gt_unphased(bcf_gt_allele(ptr[0]));
+    ptr[0] = bcf_gt_unphased(gt);    // first gt doesn't carry the phased flag per BCF specifications
+    return 2;
+}
+
 // sets GT for a single sample, ngts is the ploidy, allele
 static inline int set_gt(int32_t *ptr, int ngts, int allele)
 {
@@ -543,6 +561,8 @@ bcf1_t *process(bcf1_t *rec)
             }
             else if ( args->new_mask & GT_X_VAF )
                 changed += set_gt(ptr, ngts, args->xarr[i]);
+            else if ( args->new_mask & GT_INVPHASE )
+                changed += invert_phase_gt(ptr, ngts);
             else
                 changed += set_gt(ptr, ngts, args->new_gt);
         }
@@ -584,6 +604,8 @@ bcf1_t *process(bcf1_t *rec)
             }
             else if ( args->new_mask & GT_X_VAF )
                 changed += set_gt(args->gts + i*ngts, ngts, args->xarr[i]);
+            else if ( args->new_mask & GT_INVPHASE )
+                changed += invert_phase_gt(args->gts, ngts);
             else
                 changed += set_gt(args->gts + i*ngts, ngts, args->new_gt);
         }
@@ -620,6 +642,8 @@ bcf1_t *process(bcf1_t *rec)
             }
             else if ( args->new_mask & GT_X_VAF )
                 changed += set_gt(ptr, ngts, args->xarr[i]);
+            else if ( args->new_mask & GT_INVPHASE )
+                changed += invert_phase_gt(ptr, ngts);
             else
                 changed += set_gt(ptr, ngts, args->new_gt);
         }

test/setGT.2.1.out

@@ -0,0 +1,17 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=20,length=81195210>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depth">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	HG00100	HG00101	HG00102	HG00200	HG00201	HG00202
+20	302	.	T	TA	999	.	.	GT:AD	1|0:55,5	0|1:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	828	.	T	C	999	.	.	GT:AD	1|0:55,5	1|1:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	834	.	G	A	999	.	.	GT:AD	1|0:55,5	1|0:55,5	1|1:5,55	0|1:55,5	0|1:55,5	1|1:5,55
+20	1665	.	T	C	999	.	.	GT:AD	0|1:55,5	1|0:55,5	1|0:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	1869	.	A	T	999	.	.	GT:AD	1|0:55,5	0|0:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	2041	.	G	A	999	.	.	GT:AD	1|0:55,5	1|1:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	2220	.	G	A	999	.	.	GT:AD	1|0:55,5	0|1:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	2564	.	A	G	999	.	.	GT:AD	1|0:55,5	1|1:55,5	1|1:5,55	1|0:55,5	0|1:55,5	1|1:5,55
+20	3104	.	C	T	999	.	.	GT:AD	0|0:5,5	0|0:5,5	1|0:5,5	1|0:5,5	0|1:5,5	1|1:5,5
+20	3587	.	G	A	999	.	.	GT:AD	1|0:5,5	1|0:5,5	1|1:5,5	1|0:5,5	0|1:5,5	1|1:5,5
+20	3936	.	A	G	999	.	.	GT:AD	1|0:1,2	1|1:3,4	1|1:5,6	1|0:7,8	0|1:9,1	1|1:2,3

test/test.pl

@@ -633,6 +633,7 @@
 run_test(\&test_vcf_plugin,$opts,in=>'setGT',out=>'setGT.1.out',cmd=>'+setGT --no-version',args=>'-- -t q -n 0 -i \'GT~"." && FMT/DP=30 && GQ=150\'');
 run_test(\&test_vcf_plugin,$opts,in=>'setGT.2',out=>'setGT.2.out',cmd=>'+setGT --no-version',args=>'-- -t q -n . -i \'GT[@{QPATH}/setGT.samples.txt]="het"\'');
 run_test(\&test_vcf_plugin,$opts,in=>'setGT.2',out=>'setGT.3.out',cmd=>'+setGT --no-version',args=>'-- -t q -n . -i \'GT[@{QPATH}/setGT.samples.txt]="het" & binom(AD[@{QPATH}/setGT.samples.txt])<0.1\'');
+run_test(\&test_vcf_plugin,$opts,in=>'setGT.2',out=>'setGT.2.1.out',cmd=>'+setGT --no-version',args=>'-- -t a -n i');
 run_test(\&test_vcf_plugin,$opts,in=>'setGT.3',out=>'setGT.3.1.out',cmd=>'+setGT --no-version',args=>'-- -t a -n pM');
 run_test(\&test_vcf_plugin,$opts,in=>'setGT.3',out=>'setGT.3.2.out',cmd=>'+setGT --no-version',args=>'-- -t a -n pm');
 run_test(\&test_vcf_plugin,$opts,in=>'setGT.3',out=>'setGT.3.3.out',cmd=>'+setGT --no-version',args=>'-- -t a -n c:1');