The current GATK version is 3.8-0
Examples: Monday, today, last week, Mar 26, 3/26/04

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Get notifications!

You can opt in to receive email notifications, for example when your questions get answered or when there are new announcements, by following the instructions given here.

Got a problem?

1. Search using the upper-right search box, e.g. using the error message.
2. Try the latest version of tools.
3. Include tool and Java versions.
4. Tell us whether you are following GATK Best Practices.
5. Include relevant details, e.g. platform, DNA- or RNA-Seq, WES (+capture kit) or WGS (PCR-free or PCR+), paired- or single-end, read length, expected average coverage, somatic data, etc.
6. For tool errors, include the error stacktrace as well as the exact command.
7. For format issues, include the result of running ValidateSamFile for BAMs or ValidateVariants for VCFs.
8. For weird results, include an illustrative example, e.g. attach IGV screenshots according to Article#5484.
9. For a seeming variant that is uncalled, include results of following Article#1235.

Did we ask for a bug report?

Then follow instructions in Article#1894.

Formatting tip!

Wrap blocks of code, error messages and BAM/VCF snippets--especially content with hashes (#)--with lines with three backticks ( ``` ) each to make a code block as demonstrated here.

Jump to another community
Download the latest Picard release at
GATK version 4.beta.3 (i.e. the third beta release) is out. See the GATK4 beta page for download and details.

Forced Calls Homozygous Reference don't have GQ

torstensontorstenson Member
edited January 2013 in Ask the GATK 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


  • Mark_DePristoMark_DePristo Broad InstituteMember

    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.

  • 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 Broad InstituteMember

    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.

  • 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 Broad InstituteMember, Broadie, Dev

    You can run in EMIT_ALL_SITES mode.

  • 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 Broad InstituteMember, Broadie, Dev

    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?

  • evansbevansb Member
    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
  • ebanksebanks Broad InstituteMember, Broadie, Dev

    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.

Sign In or Register to comment.