GenotypeGvcfs mangles the input phase
open | Created 2019-02-26 | Last updated 2019-03-13| Posted by nh13 | See in Github


VariantCalling bug tools


From my review of the code, we have:

  1. GenotypeGvcfs first merges overlapping variant contexts, resulting in genotypes with no-calls.
  2. Next, the samples are re-genotyped.

The problem is that phasing is lost in (1), where the alleles are converted to no-calls, and then the phase is never properly reconstructed in (2).

I tried diving into the code, and made it to GATKVariantContextUtils.makeGenotypeCall to try to figure out how to insert phase information, but I went cross-eyed. @ldgauthier @davidbenjamin @droazen what do you think?

HaplotypeCaller output:

chr7	92147580	.	C	G,<NON_REF>	1909.77	.	BaseQRankSum=-0.304;ClippingRankSum=-2.049;DP=108;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQandDP=388800,108;REF_BASES=CTAGAACCAACAGACGAAAAG;ReadPosRankSum=1.577	GT:AD:DP:F1R2:F2R1:GQ:PGT:PID:PL:PS:SB	0|1:53,55,0:108:24,25,0:29,30,0:99:0|1:92147580_C_G:1938,0,2058,2097,2222,4319:92147580:24,29,25,30
chr7	92147581	.	A	<NON_REF>	.	.	END=92147583	GT:DP:GQ:MIN_DP:PL	0/0:108:99:108:0,120,1800
chr7	92147584	.	C	T,<NON_REF>	1908.77	.	BaseQRankSum=0.609;ClippingRankSum=1.815;DP=109;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQandDP=392400,109;REF_BASES=AACCAACAGACGAAAAGATCA;ReadPosRankSum=-1.456	GT:AD:DP:F1R2:F2R1:GQ:PGT:PID:PL:PS:SB	1|0:56,53,0:109:26,24,0:30,29,0:99:1|0:92147580_C_G:1937,0,2138,2106,2297,4403:92147580:26,30,24,29

This gives the genotype CAGAT/GAGAC

GenotypeGvcfs output:

chr7	92147580	.	C	G	1930.60	.	AC=1;AF=0.500;AN=2;BaseQRankSum=-3.040e-01;ClippingRankSum=-2.049e+00;DP=108;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=17.88;ReadPosRankSum=1.58;SOR=0.687	GT:AD:DP:F1R2:F2R1:GQ:PGT:PID:PL:PS	0|1:53,55:108:24,25,0:29,30,0:99:0|1:92147580_C_G:1938,0,2058:92147580
chr7	92147584	.	C	T	1929.60	.	AC=1;AF=0.500;AN=2;BaseQRankSum=0.609;ClippingRankSum=1.82;DP=109;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=17.70;ReadPosRankSum=-1.456e+00;SOR=0.738	GT:AD:DP:F1R2:F2R1:GQ:PGT:PID:PL:PS	0|1:56,53:109:26,24,0:30,29,0:99:1|0:92147580_C_G:1937,0,2138:92147580

This gives the genotype CAGAC/GAGAT


Return to top