Bug Bulletin: we have identified a bug that affects indexing when producing gzipped VCFs. This will be fixed in the upcoming 3.2 release; in the meantime you need to reindex gzipped VCFs using Tabix.

Forced Calls Homozygous Reference don't have GQ

torstensontorstenson Posts: 2Member
edited January 2013 in Ask the team

I am implementing a tool that uses reads to identify potential sites for cancer, but I need a reliable genotype call from non-cancerous tissue. Currently, I'm trying to use GATK to make those calls, filter out the bad, and then use the remaining calls when I look over the raw reads in the cancer tissue.

The problem I'm finding is that GATK doesn't provide a GQ for variants that are homozygous for the reference, and I don't know how to correlate the QUAL with the GQs from those homozygous non-ref. Is there a way that I'm overlooking to have GATK provide the GQ for homozygous reference non-variant sites?

Post edited by Geraldine_VdAuwera on

Comments

  • Mark_DePristoMark_DePristo Posts: 153Administrator, GSA Member admin

    You need to jointly call the tumor and normal samples together to get what you want. You cannot correctly calculate the GQ without the alternative allele from the tumor.

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

  • torstensontorstenson Posts: 2Member

    Thank you for the response we are trying to decide if that will interfere with our algorithm in a significant way.

    If, however, we were to not want to call them jointly, is the QUAL on the same scale as the GQ and would we be able to filter those homozygous references using QUAL? We have observed that a recent version of GATK will use QUAL=GQ-30 for heterozygous, and wonder if that the homozygous/ref might have a similar correlation.

    Thank you very much for any thoughts on this.

  • Mark_DePristoMark_DePristo Posts: 153Administrator, GSA Member admin

    This relationship only holds for the simple single sample case for SNPs. You will immediately have problems with anything more complex if you don't jointly call your samples.

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

  • evansbevansb Posts: 4Member

    Does anyone have a suggestion for how one can make calls of homozygous genotypes jointly across samples (if not with GATK, perhaps with another software)? It seems that if there is not a known alternative allele, one could simply consider each of the 3 non-reference nucleotides as an equally probably SNP (?). Thanks.

  • ebanksebanks Posts: 671GSA Member mod

    You can run in EMIT_ALL_SITES mode.

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

  • evansbevansb Posts: 4Member

    Thanks for getting back to me Eric. I did try using "EMIT_ALL_SITES" but it appears I am still getting only variable sites in the vcf file. I've pasted a sample of the file below. I'd like each homozygous site to be output as well (for examples positions on chr1 22115, 22116 ... 22127 are probably homozygous but are not listed). Am I missing something? Thanks again.

    chr1 22114 . A G 32304.02 PASS AC=18;AF=1.00;AN=18;ActiveRegionSize=263;DP=323;EVENTLENGTH=0;FS=0.000;MLEAC=18;MLEAF=1.00;MQ=58.98;NVH=5;NumHapAssembly=14;NumHapEval=13;QD=10 0.01;QDE=20.00;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 1/1:99:3789,468,0 1/1:99:2935,181,0 1/1:99:3747,376,0 1/1:99:2721,238,0 1/1:99:3665,615,0 1/1:99:3444,205,0 1/1:99:3757,510,0 1/1:99:2816,155,0 1/1:99:5430,794,0 chr1 22128 . C T 1228.47 PASS AC=1;AF=0.056;AN=18;ActiveRegionSize=263;ClippingRankSum=-1.665;DP=322;EVENTLENGTH=0;FS=0.000;MLEAC=1;MLEAF=0.056;MQ=59.12;MQRankSum=-0.065;NVH=4;NumHa pAssembly=14;NumHapEval=13;QD=26.14;QDE=6.54;ReadPosRankSum=-1.097;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/0:99:0,465,2305 0/0:99:0,179,1777 0/0:99:0,373,2402 0/0:99:0,237,17 52 0/0:99:0,612,2480 0/0:99:0,205,2278 0/0:99:0,507,2573 0/0:99:0,155,1578 0/1:99:1268,0,652 chr1 22129 . T C 32304.02 PASS AC=18;AF=1.00;AN=18;ActiveRegionSize=263;DP=322;EVENTLENGTH=0;FS=0.000;MLEAC=18;MLEAF=1.00;MQ=59.12;NVH=5;NumHapAssembly=14;NumHapEval=13;QD=10 0.32;QDE=20.06;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 1/1:99:3789,468,0 1/1:99:2935,181,0 1/1:99:3747,376,0 1/1:99:2721,238,0 1/1:99:3665,615,0 1/1:99:3444,205,0 1/1:99:3757,510,0 1/1:99:2816,155,0 1/1:99:5430,794,0 chr1 22179 . C T 1009.75 PASS AC=5;AF=0.278;AN=18;ActiveRegionSize=263;ClippingRankSum=0.496;DP=135;EVENTLENGTH=0;FS=1.904;MLEAC=5;MLEAF=0.278;MQ=59.80;MQRankSum=-0.882;NVH=5;NumHap Assembly=14;NumHapEval=13;QD=13.83;QDE=2.77;ReadPosRankSum=-0.634;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/1:99:301,0,625 0/0:99:0,185,1267 0/1:99:200,0,641 0/1:76:76,0,787 0/1:99:197,0,1269 0/0:99:0,139,1817 0/1:99:297,0,394 0/0:92:0,92,962 0/0:99:0,797,2680 chr1 22265 . C A 1109.16 PASS AC=4;AF=0.222;AN=18;ActiveRegionSize=47;ClippingRankSum=-0.082;DP=159;EVENTLENGTH=0;FS=0.000;MLEAC=4;MLEAF=0.222;MQ=59.56;MQRankSum=-0.153;NVH=2;NumHap Assembly=4;NumHapEval=4;QD=22.64;QDE=11.32;ReadPosRankSum=-1.997;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/0:99:0,314,1010 0/1:99:241,0,278 0/0:99:0,331,1136 0/0:99:0,110,38 3 0/1:99:260,0,412 1/1:56:673,56,0 0/0:57:0,57,1343 0/0:54:0,54,658 0/0:69:0,69,847 chr1 22266 . C T 1537.59 PASS AC=6;AF=0.333;AN=18;ActiveRegionSize=47;ClippingRankSum=0.208;DP=157;EVENTLENGTH=0;FS=0.000;MLEAC=6;MLEAF=0.333;MQ=59.56;MQRankSum=1.480;NVH=2;NumHapAs sembly=4;NumHapEval=4;QD=18.98;QDE=9.49;ReadPosRankSum=0.688;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/1:99:254,0,364 0/0:99:0,292,864 0/1:99:363,0,294 0/1:86:86,0,153 0/1:99: 233,0,476 0/0:60:0,60,1351 1/1:57:670,57,0 0/0:51:0,51,623 0/0:68:0,68,810 chr1 22294 . G A 6698.02 PASS AC=18;AF=1.00;AN=18;ActiveRegionSize=47;DP=113;EVENTLENGTH=0;FS=0.000;MLEAC=18;MLEAF=1.00;MQ=58.56;NVH=2;NumHapAssembly=4;NumHapEval=4;QD=59.27;QDE=29. 64;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 1/1:99:709,205,0 1/1:78:589,78,0 1/1:99:914,272,0 1/1:99:339,109,0 1/1:99:710,260,0 1/1:66:1039,66,0 1/1:75:1315,75, 0 1/1:35:416,35,0 1/1:54:667,54,0 chr1 22379 . C T 196.42 PASS AC=2;AF=0.111;AN=18;ActiveRegionSize=165;ClippingRankSum=-0.970;DP=39;EVENTLENGTH=0;FS=0.000;MLEAC=2;MLEAF=0.111;MQ=59.41;MQRankSum=0.091;NVH=2;NumHapA ssembly=4;NumHapEval=4;QD=28.06;QDE=14.03;ReadPosRankSum=0.165;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/0:9:0,9,92 0/0:6:0,6,59 0/0:24:0,24,267 0/0:18:0,18,187 0/0:6:0,6,73 0/0:3:0,3,36 0/0:12:0,12,145 0/0:18:0,18,214 1/1:21:249,21,0 chr1 22460 . C G 53.26 PASS AC=2;AF=0.333;AN=6;ActiveRegionSize=165;ClippingRankSum=-0.937;DP=7;EVENTLENGTH=0;FS=0.000;MLEAC=2;MLEAF=0.333;MQ=56.62;MQRankSum=0.937;NVH=2;NumHapAss embly=4;NumHapEval=4;QD=17.75;QDE=8.88;ReadPosRankSum=-0.937;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL ./. ./. 1/1:9:93,9,0 0/0:6:0,6,67 ./. ./. ./. ./. 0/0:6:0,6,64 chr1 33823 . T G 4151.02 PASS AC=18;AF=1.00;AN=18;ActiveRegionSize=214;DP=104;EVENTLENGTH=0;FS=0.000;MLEAC=18;MLEAF=1.00;MQ=37.03;NVH=3;NumHapAssembly=10;NumHapEval=10;QD=39.91;QDE= 13.30;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 1/1:41:449,41,0 1/1:45:692,45,0 1/1:33:490,33,0 1/1:33:342,33,0 1/1:24:296,24,0 1/1:36:519,36,0 1/1:24:258,24,0 1/1:32:521,32,0 1/1:45:584,45,0 chr1 33836 . G . 13.72 LowQual AN=18;ActiveRegionSize=214;DP=124;EVENTLENGTH=0;MQ=42.14;NVH=0;NumHapAssembly=10;NumHapEval=10;TYPE=INDEL;extType=NO_VARIATION;set=FilteredInAll GT 0/0 0/0 0/0 0/0 0/0 0/0 0/0 0/0 0/0 chr1 33962 . T C 4828.02 PASS AC=16;AF=0.889;AN=18;ActiveRegionSize=214;DP=151;EVENTLENGTH=0;FS=0.000;MLEAC=16;MLEAF=0.889;MQ=44.84;NVH=3;NumHapAssembly=10;NumHapEval=10;QD=31.97;QD E=10.66;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/1:72:389,0,72 1/1:54:589,54,0 1/1:48:534,48,0 1/1:51:582,51,0 1/1:42:459,42,0 0/1:40:606,0,40 1/1:39:428,39,0 1/1:42:450,42,0 1/1:69:791,69,0 chr1 34052 . G A 127.73 PASS AC=2;AF=0.111;AN=18;ActiveRegionSize=240;ClippingRankSum=-0.386;DP=88;EVENTLENGTH=0;FS=0.000;MLEAC=2;MLEAF=0.111;MQ=55.73;MQRankSum=1.460;NVH=3;NumHapA ssembly=10;NumHapEval=10;QD=7.51;QDE=2.50;ReadPosRankSum=0.248;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/1:99:158,0,390 0/0:99:0,189,811 0/0:78:0,78,665 0/0:99:0,189,997 0/0:99: 0,173,658 0/1:18:18,0,700 0/0:99:0,114,671 0/0:99:0,204,723 0/0:99:0,154,1116 chr1 34106 . A G 1714.70 PASS AC=9;AF=0.500;AN=18;ActiveRegionSize=240;DP=49;EVENTLENGTH=0;FS=0.000;MLEAC=9;MLEAF=0.500;MQ=52.76;NVH=3;NumHapAssembly=10;NumHapEval=10;QD=34.99;QDE=1 1.66;TYPE=SNP;extType=SNP;set=variant GT:GQ:PL 0/1:99:254,0,390 0/1:99:200,0,400 0/1:69:69,0,467 0/1:99:259,0,513 0/1:99:202,0,380 0/1:99:170,0,481 0/1:99:129,0,36 4 0/1:99:263,0,377 0/1:99:228,0,702

  • ebanksebanks Posts: 671GSA Member mod

    This doesn't look like the output of EMIT_ALL_SITES. Which version of the GATK are you using and what's your command-line?

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

  • evansbevansb Posts: 4Member
    edited November 2012

    I'm using version 2.1-13-g1706365.

    Here's my command-line:

    java -Xmx512m -jar GenomeAnalysisTK.jar -T HaplotypeCaller -L chr1 
        -R /work/ben/macaque_RAD_TAGs/macaque_genome/hardMask/macaque_masked_chromosomes.fasta 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PF515sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM561sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM565sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM566sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM567sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM582sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM584sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM592sorted_dedup.realigned.bam 
        -I /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/PM602sorted_dedup.realigned.bam 
        -out_mode EMIT_ALL_SITES -fullHaplotype 
        -o /work/ben/macaque_RAD_TAGs/johnomics-RADtools-a814a6b/fuzzymacaques/chr1_1.vcf 
        -maxAltAlleles 4
    
    Post edited by Geraldine_VdAuwera on
  • ebanksebanks Posts: 671GSA Member mod

    You need to use the Unified Genotyper for EMIT_ALL_SITES to work. This is an issue with the documentation for the Haplotype Caller that we already know about and will aim to fix soon.

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

Sign In or Register to comment.