Bug Bulletin: The GenomeLocPArser error in SplitNCigarReads has been fixed; if you encounter it, use the latest nightly build.

Strange GT and GQ when using sample_ploidy set to 1

aunderwoaunderwo Posts: 3Member
edited October 2012 in Ask the GATK team

I am experiencing a problem with the ploidy 1 option. Having used GATK2 unified genotyper with the params --sample_ploidy 1 --genotype_likelihoods_model BOTH -rf BadCigar I get the following line in a vcf file (see sample in bold)

Staphylococcus 1553115 . A G 24454.01 . AC=13;AF=0.813;AN=16;BaseQRankSum=1.072;DP=1040;Dels=0.00;FS=32.822;HaplotypeScore=3.3463;MLEAC=13;MLEAF=0.813;MQ=40.20;MQ0=47;MQRankSum=-10.543;QD=32.13;ReadPosRankSum=-1.148;SB=-9.076e+03 GT:AD:DP:GQ:MLPSAC:MLPSAF:PL 1:0,29:29:99:1:1.00:1015,0 1:0,62:62:99:1:1.00:2053,0 1:0,106:106:99:1:1.00:3210,0 1:0,102:102:99:1:1.00:3305,0 1:0,88:88:99:1:1.00:2750,0 1:0,41:41:99:1:1.00:1324,0 1:0,76:76:99:1:1.00:2448,0 1:0,39:39:99:1:1.00:1303,0 0:64,40:104:99:0:0.00:0,1334 1:0,41:41:99:1:1.00:1373,0 1:0,49:49:99:1:1.00:1668,0 0:72,50:122:99:0:0.00:0,1258 1:0,59:59:99:1:1.00:1852,0 1:0,38:38:99:1:1.00:1192,0 1:0,31:31:99:1:1.00:961,0 0:53,0:53:99:0:0.00:0,1633

The sample in bold is called as WT (genotype 0) with a high GQ despite there being 72 reads of genotype 0 and 50 of genotype 1. Examining the bam file suggests that this is a mapping error in a repetitive phage region

If I set ploidy to be 2 the equivalent line in the resulting vcf file is

Staphylococcus 1553115 . A G 24788.02 . AC=28;AF=0.875;AN=32;BaseQRankSum=0.947;DP=1040;Dels=0.00;FS=32.822;HaplotypeScore=3.3463;InbreedingCoeff=0.4286;MLEAC=28;MLEAF=0.875;MQ=40.20;MQ0=47;MQRankSum=-10.096;QD=25.11;ReadPosRankSum=-1.177;SB=-9.871e+03 GT:AD:DP:GQ:PL 1/1:0,29:29:81:986,81,0 1/1:0,62:62:99:1895,156,0 1/1:0,106:106:99:2992,247,0 1/1:0,102:102:99:3169,268,0 1/1:0,88:88:99:2452,193,0 1/1:0,41:41:99:1243,99,0 1/1:0,76:76:99:2283,193,0 1/1:0,39:39:99:1233,105,0 0/1:64,40:104:99:886,0,1706 1/1:0,41:41:99:1298,108,0 1/1:0,49:49:99:1581,129,0 0/1:72,50:122:99:1235,0,2126 1/1:0,59:59:99:1649,132,0 1/1:0,38:38:87:1065,87,0 1/1:0,31:31:69:821,69,0 0/0:53,0:53:99:0,138,1588

As can be seen from the bold text, the same position is called as heterozygote which based on the number of the reads mapping would be likley except for the fact this is a bacterial haploid genome. Previously I would have discarded this since the heterozygous call indicates mis-mapping as the bam file confirms. I had been hoping to use the sample_polidy option set to 1 for bacterial genomes but this result concerns me. I could obviously filter based on AD but the wonder why the sample was given a high GQ when the polidy is set to 1 and the AD suggests the call of genotype 0 should be doubted.

Any suggestions on what is going on here?? Many thanks

Anthony

Anthony

Post edited by Geraldine_VdAuwera on

Answers

  • Mark_DePristoMark_DePristo Posts: 153Administrator, GATK Developer admin

    This is exactly the behavior to expect by telling the UG to use ploidy 1. Even at sites where there's multiple haplotypes -- perhaps because of alignment artifacts or perhaps polyclonality in your samples -- it must choose the best of 0 or 1 genotypes.

    -- Mark A. DePristo, Ph.D. Co-Director, Medical and Population Genetics Broad Institute of MIT and Harvard

  • aunderwoaunderwo Posts: 3Member
    edited October 2012

    Thanks for the reply. I understand that it must call either 0 or 1 but why does it assign a high GQ? Surely it can't be sure it's 0 in this case when 72 reads are 0 and 50 are 1? I would have hoped that the GQ would have been lower so a filter could be applied.

    Post edited by aunderwo on
  • Mark_DePristoMark_DePristo Posts: 153Administrator, GATK Developer admin

    GQ is the diff between ref and alt genotype likelihoods, and you can see they are 0 and 1258 here, so it's very high quality. It's because there's a much better than of being ref than alt in this case. Personally I prefer to use the diploid model, and then use the het states as an error sink.

    -- Mark A. DePristo, Ph.D. Co-Director, Medical and Population Genetics Broad Institute of MIT and Harvard

  • aunderwoaunderwo Posts: 3Member

    Thanks for the explanation. If it's not too much trouble please can you tell me how the likelihoods are calculated and how these then make the final GQ? Would you think filtering those positions where the fraction of reads matching the chosen genotype is >0.9 (using the values in AD) would be a possible way of using a haploid model.

    Thanks again for your responses.

  • Mark_DePristoMark_DePristo Posts: 153Administrator, GATK Developer admin

    That's too detailed for the forum. Have a look at our papers and, if you want, at the source code.

    -- Mark A. DePristo, Ph.D. Co-Director, Medical and Population Genetics Broad Institute of MIT and Harvard

  • ebanksebanks Posts: 683GATK Developer mod

    I just want to add a bit to Mark's excellent answer above: remember that the AD gives the raw unfiltered counts, but it is very possible that the reads showing the G allele are extremely low quality (which is why the likelihoods are skewed). This is addressed directly in the documentation for the genotyper.

    Eric Banks, PhD -- Senior Group Leader, MPG Analysis, Broad Institute of Harvard and MIT

Sign In or Register to comment.