Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

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.

GATK 3.7 HaplotypeCaller out vcf format

Hi everyone,

I just want to confirm what I got back from HaplotypeCaller (GATK 3.7) is correct format, it is VCF4.2 format, however, after following the best practice, I have a vcf file with the GT field with "/" not "|", according to the VCF4.2 format, the "/" means unphased genotype, however, haplotypeCaller should provide haplotype calls instead of genotypes right? Should I use the results as unphased or phased?

Here are couple lines from the vcf file generated.

fileformat=VCFv4.2

...

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878

1 10492 . C T 39.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.385;ClippingRankSum=0.000;DP=7;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=50.37;MQRankSum=-1.282;QD=6.63;ReadPosRankSum=0.000;SOR=0.693 GT:AD:DP:GQ:PL 0/1:3,3:6:66:68,0,66

Thanks a lot!

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @yzharold
    Hi,

    The genotype will have | only if it is in phase with another variant. In your example, that variant is not phased with another variant. Do you have an example of variants you think should be phased but are not?

    Can you post an example? Please include the BAM file and VCF records.

    Thanks,
    Sheila

  • yzharoldyzharold Member

    Thanks for your explanation, the bam file was downloaded from 1000 Genomes Project, the famous NA12878 sample, it is a huge bam file from exome sequencing and mapped reads using illumina technology. I have used the following command to call the variants using HaplotypeCaller. I also attached the example vcf file to show multiple variants were called.

    Command line:

    java -Xmx10g -Djava.io.tmpdir=/tmp -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R /pathto/reference -I /pathto/NA12878.bam -o raw.snps.indels.SIngleSample.vcf

    My assumption is that the HaplotypeCaller will not only call all the variants, but also provide the phased genotypes for each loci with "|" instead of "/" for this single sample?

    Sorry for the inconvenience, due to some limited uploading capability, I have to copy all the lines of variants called to give you an example, if it is not clear, I could email you the example I created (vcf) file with your the provided.

    fileformat=VCFv4.2
    ...

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878

    1 10492 . C T 39.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.385;ClippingRankSum=0.000;DP=7;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=50.37;MQRankSum=-1.2\
    82;QD=6.63;ReadPosRankSum=0.000;SOR=0.693 GT:AD:DP:GQ:PL 0/1:3,3:6:66:68,0,66
    1 14907 . A G 409.77 . AC=2;AF=1.00;AN=2;DP=16;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=26.34;QD=25.61;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,16:16:48\
    :438,48,0
    1 14930 . A G 580.77 . AC=2;AF=1.00;AN=2;DP=23;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=25.54;QD=25.25;SOR=0.963 GT:AD:DP:GQ:PL 1/1:0,23:23:69\
    :609,69,0
    1 15118 . A G 27.75 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=23.55;QD=13.87;SOR=2.303 GT:AD:DP:GQ:PL 1/1:0,2:2:6:55,6,0
    1 50481 . GGT G 28.71 . AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=36.67;QD=14.35;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:65,6,0
    1 55299 . C T 57.28 . AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=46.36;QD=19.09;SOR=1.179 GT:AD:DP:GQ:PL 1/1:0,3:3:9:85,9,0
    1 104176 . C G 62.74 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=33.24;QD=31.37;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:90,6,0
    1 104186 . T C 62.74 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=33.24;QD=31.37;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:90,6,0
    1 104208 . C G 18.59 . AC=2;AF=1.00;AN=2;DP=1;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=29.00;QD=18.59;SOR=1.609 GT:AD:DP:GQ:PL 1/1:0,1:1:3:45,3,0
    1 234481 . T A 36.74 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=34.48;QD=18.37;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:64,6,0
    1 255909 . A ATG 33.71 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=24.02;QD=16.85;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:70,6,0
    1 528040 . C A 12.08 . AC=1;AF=0.500;AN=2;BaseQRankSum=-0.967;ClippingRankSum=0.000;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=37.96;MQRankSum=0.9\
    67;QD=4.03;ReadPosRankSum=0.000;SOR=0.223 GT:AD:DP:GQ:PL 0/1:1,2:3:21:40,0,21
    1 565286 . C T 2515.77 . AC=2;AF=1.00;AN=2;BaseQRankSum=-0.498;ClippingRankSum=0.000;DP=53;ExcessHet=3.0103;FS=30.881;MLEAC=2;MLEAF=1.00;MQ=22.45;MQRankSum=-6.\
    926;QD=34.24;ReadPosRankSum=0.164;SOR=7.784 GT:AD:DP:GQ:PL 1/1:2,48:50:76:2544,76,0
    1 565295 . A G 4468.77 . AC=2;AF=1.00;AN=2;BaseQRankSum=-0.614;ClippingRankSum=0.000;DP=142;ExcessHet=3.0103;FS=7.360;MLEAC=2;MLEAF=1.00;MQ=25.79;MQRankSum=-1.\
    442;QD=32.62;ReadPosRankSum=-1.599;SOR=0.676 GT:AD:DP:GQ:PL 1/1:2,135:137:99:4497,282,0
    1 565406 . C T 42.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.126;ClippingRankSum=0.000;DP=6;ExcessHet=3.0103;FS=3.979;MLEAC=1;MLEAF=0.500;MQ=36.79;MQRankSum=-0.6\
    74;QD=7.13;ReadPosRankSum=0.189;SOR=2.800 GT:AD:DP:GQ:PL 0/1:3,3:6:71:71,0,85
    1 565419 . C G 21.80 . AC=1;AF=0.500;AN=2;BaseQRankSum=1.282;ClippingRankSum=0.000;DP=6;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=43.00;MQRankSum=-0.2\
    53;QD=3.63;ReadPosRankSum=1.834;SOR=0.223 GT:AD:DP:GQ:PL 0/1:4,2:6:50:50,0,113
    1 567242 . G A 21.77 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=10.88;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:49,6,0
    1 569933 . G A 52.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.652;ClippingRankSum=0.000;DP=7;ExcessHet=3.0103;FS=3.136;MLEAC=1;MLEAF=0.500;MQ=45.61;MQRankSum=0.84\
    2;QD=7.54;ReadPosRankSum=-0.484;SOR=0.916 GT:AD:DP:GQ:PL 0/1:4,3:7:81:81,0,102
    1 610242 . G A 23.79 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.060;ClippingRankSum=0.000;DP=8;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=37.00;MQRankSum=0.00\
    0;QD=3.40;ReadPosRankSum=0.366;SOR=0.223 GT:AD:DP:GQ:PL 0/1:5,2:7:52:52,0,157
    1 663097 . G C 39.74 . AC=2;AF=1.00;AN=2;DP=3;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=29.85;QD=19.87;SOR=2.303 GT:AD:DP:GQ:PL 1/1:0,2:2:6:67,6,0
    1 672209 . G A 21.77 . AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=39.00;QD=10.88;SOR=0.693 GT:AD:DP:GQ:PL 1/1:0,2:2:6:49,6,0

    Thanks a lot!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @yzharold
    Hi,

    Can you post the BAM file and bamout file for sites you think should be phased but are not?

    Thanks,
    Sheila

  • yzharoldyzharold Member

    Hi Sheila,

    I think I figure out the issue, I found http://gatkforums.broadinstitute.org/gatk/discussion/3061/haplotype-score-and-phasing helpful explaining the HC not emit haplotypes, but try to best estimate the genotypes with high confidence. And you need extra steps to phase your sample, such as family information. As my understanding, haplotypeCaller will take advantage of the haplotype information for calling variants, but not necessary phasing the individual sample by default. Am I correct?

    Do you think it is a good explanation or I am off? The bam file is huge and I could not upload from my work.

    Sorry for the inconvenience.

    Thanks,

    Harold

Sign In or Register to comment.