Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Unexpected heterozygous genotype produced by GenotypeGVCFs

mlindermmlinderm Member
edited January 2015 in Ask the GATK team

I am seeing an unexpected heterozygous genotype in the GenotypeGVCFs output. The region is at the edge of a long homopolymer, which presumably plays a role, but it is still not clear where the unexpected heterozygote is coming from. The final VCF is below. The unexpected call is 2:152402516A>T in Sample1. The call is 0/1, but I don't see any evidence for where that alternate allele is coming from. The input VCFs are below. Any suggestions? How is GenotypeGVCFs interpreting this situation? I am running the latest GATK 3.3. Jar.

Final VCF

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1        Sample2
2       152402515       .       T       TA      86.23   .       AC=2;AF=0.500;AN=4;BaseQRankSum=-5.700e-02;ClippingRankSum=0.321;DP=131;FS=0.000;GQ_MEAN=166.50;GQ_STDDEV=44.55;MLEAC=1;MLEAF=0.250;MQ=59.89;MQ0=0;MQRankSum=2.42;NCC=0;QD=1.31;ReadPosRankSum=0.649       GT:AD:DP:GQ:PL  0/1:10,8:34:99:203,0,198        0/1:37,11:56:99:135,0,1429
2       152402516       .       A       T       380.44  .       AC=2;AF=0.500;AN=4;BaseQRankSum=-3.647e+00;ClippingRankSum=-2.597e+00;DP=130;FS=3.866;GQ_MEAN=204.50;GQ_STDDEV=136.47;MLEAC=2;MLEAF=0.500;MQ=59.73;MQ0=0;MQRankSum=0.754;NCC=0;QD=5.76;ReadPosRankSum=1.43 GT:AD:DP:GQ:PL  0/1:10,0:34:99:108,0,157        0/1:43,23:66:99:301,0,1142

Sample1 VCF

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
2       152402511       .       G       <NON_REF>       .       .       END=152402514   GT:DP:GQ:MIN_DP:PL      0/0:38:57:33:0,51,765
2       152402515       .       TA      T,TAA,TAAA,TAAAA,<NON_REF>      0.10    .       BaseQRankSum=-0.057;ClippingRankSum=0.321;DP=55;MLEAC=0,0,1,0,0;MLEAF=0.00,0.00,0.500,0.00,0.00;MQ=59.89;MQ0=0;MQRankSum=-1.002;ReadPosRankSum=-0.646      GT:AD:DP:GQ:PL:SB       0/3:10,8,8,7,1,0:34:21:224,189,625,21,87,219,0,109,141,388,143,293,250,447,708,116,220,156,201,330,273:0,10,0,7
2       152402517       .       A       <NON_REF>       .       .       END=152402517   GT:DP:GQ:MIN_DP:PL      0/0:55:0:55:0,0,1431
2       152402518       .       A       <NON_REF>       .       .       END=152402621   GT:DP:GQ:MIN_DP:PL      0/0:49:99:37:0,101,1442

Sample2 VCF

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample2
2       152402513       .       C       CT,<NON_REF>    0       .       BaseQRankSum=-0.168;ClippingRankSum=0.437;DP=48;MLEAC=0,0;MLEAF=0.00,0.00;MQ=60.00;MQ0=0;MQRankSum=-0.772;ReadPosRankSum=0.034  GT:AD:DP:GQ:PL:SB  0/0:35,2,0:37:99:0,103,1228,143,1234,1275:6,29,0,0
2       152402514       .       T       <NON_REF>       .       .       END=152402514   GT:DP:GQ:MIN_DP:PL      0/0:63:99:63:0,114,1710
2       152402515       .       TA      T,TAA,TAAA,<NON_REF>    78.73   .       BaseQRankSum=-1.332;ClippingRankSum=-0.242;DP=76;MLEAC=0,1,0,0;MLEAF=0.00,0.500,0.00,0.00;MQ=59.74;MQ0=0;MQRankSum=2.423;ReadPosRankSum=0.649      GT:AD:DP:GQ:PL:SB       0/2:37,4,11,4,0:56:99:135,256,1876,0,804,1429,116,1061,1435,1759,157,933,805,963,890:2,35,0,11
2       152402516       rs199791504     A       T,<NON_REF>     272.77  .       BaseQRankSum=-3.647;ClippingRankSum=-2.597;DB;DP=75;MLEAC=1,0;MLEAF=0.500,0.00;MQ=59.73;MQ0=0;MQRankSum=0.754;ReadPosRankSum=1.426GT:AD:DP:GQ:PL:SB        0/1:43,23,0:66:99:301,0,1142,429,1210,1639:5,38,1,22
2       152402517       .       A       T,<NON_REF>     0       .       BaseQRankSum=-0.810;ClippingRankSum=-0.661;DP=79;MLEAC=0,0;MLEAF=0.00,0.00;MQ=59.75;MQ0=0;MQRankSum=-1.160;ReadPosRankSum=1.259 GT:AD:DP:GQ:PL:SB  0/0:67,4,0:71:99:0,122,2156,201,2168,2247:6,61,0,0
2       152402518       .       A       <NON_REF>       .       .       END=152402621   GT:DP:GQ:MIN_DP:PL      0/0:67:99:48:0,105,1710
Tagged:

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mlinderm‌

    Hi,

    This looks like an issue that has been fixed in the nightly builds. Please try again with the latest nightly build, and let us know if it still gives this error.
    https://www.broadinstitute.org/gatk/nightly

    Thanks,
    Sheila

  • mlindermmlinderm Member

    Thanks @Sheila‌ but unfortunately using GenotypeGVCFs from the nightly build (nightly-2015-01-13-g22d6966) did not fix the problem. I am still seeing the unexpected heterozygous call. Does the fix also require regenerating the GVCF files themselves? I am still getting the following:

    #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1        Sample2
    2       152402515       .       T       TA      330.44  .       AC=2;AF=0.500;AN=4;BaseQRankSum=-5.700e-02;ClippingRankSum=0.321;DP=131;FS=0.000;GQ_MEAN=166.50;GQ_STDDEV=44.55;MLEAC=1;MLEAF=0.250;MQ=59.89;MQ0=0;MQRankSum=2.42;NCC=0;QD=5.01;ReadPosRankSum=0.649;SOR=0.930     GT:AD:DP:GQ:PL  0/1:10,8:34:99:203,0,198        0/1:37,11:56:99:135,0,1429
    2       152402516       .       A       T       380.44  .       AC=2;AF=0.500;AN=4;BaseQRankSum=-3.647e+00;ClippingRankSum=-2.597e+00;DP=130;FS=3.866;GQ_MEAN=204.50;GQ_STDDEV=136.47;MLEAC=2;MLEAF=0.500;MQ=59.73;MQ0=0;MQRankSum=0.754;NCC=0;QD=5.76;ReadPosRankSum=1.43;SOR=1.476       GT:AD:DP:GQ:PL  0/1:10,0:34:99:108,0,157        0/1:43,23:66:99:301,0,1142
    
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mlinderm‌

    Hi,

    Can you please post IGV screenshots of the region from the input bam files and the bamout files?

    Thanks,
    Sheila

  • mlindermmlinderm Member

    @Sheila‌ Can do... But what role do the BAMs play in this issue? Everything is happening at the GVCF level.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mlinderm‌

    Yes, I am wondering if this is related to another issue I have to report. It seems like the DP/AD is getting messed up in the transfer from GVCF to final vcf. Notice that in the GVCF it shows the AD to be 10,7 for the genotype, but in the final vcf, it shows an AD of 10,0 which is weird because the genotype is 0/1. I realize the alternate allele is not correct, but I am curious what the bam out files looks like compared to the input bam file.

    Thanks for being patient! I will probably ask you to submit a bug report after I see the bam file.

    -Sheila

  • mlindermmlinderm Member

    Sample1 is AD=10,7 for the insertion, but not for the SNV at issue (at 152402516). When that SNV is "brought over" from Sample2 GVCF, I don't think Sample1 should be called heterozygous at that site. The AD=10,0 at that site in Sample1 seems to make sense (or as much sense as possible given the homopolymer).

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mlinderm‌

    Right, that is why I wanted to see the bamout files for the position. Sometimes reads get moved during reassembly. And, it seems like GVCF mode is not noting that change. I agree the call for the alternate allele is incorrect.

    In any case, please submit your files for us to debug locally if you can. Instructions are here: http://gatkforums.broadinstitute.org/discussion/1894/how-do-i-submit-a-detailed-bug-report

    Thanks,
    Sheila

  • SteveLSteveL BarcelonaMember ✭✭

    Hi @Sheila and @mlinderm - did you ever get to the bottom of this issue -- I am seeing some similar issues with GenotypeGVCFs - I am not sure whether such issues should disappear following VQSR? Thanks, Steve

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @SteveL
    Hi,

    I did put in a bug report, but it has not been fixed yet. You can follow this thread for updates: http://gatkforums.broadinstitute.org/discussion/5044/genotypegvcfs-unexpected-homozygous-heterozygous-genotypes-near-indels

    -Sheila

Sign In or Register to comment.