We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Missing DP field in Haplotype Caller produced vcf file

Hi, so I have some data that I ran through both Unified Genotype AND Haplotype Caller (because I wanted to see the difference in the number and quality of SNPs returned). I have merged both vcf files (using a simple awk command) to produce one that contains concordant SNPs between the two. However, I am now about to filter on GQ and DP except that I don't have a DP field in the file containing the SNPs that both programs called. Neither does it appear in the Haplotype Caller file. It DOES appear in the Unified Genotyper file though. I really would like to continue using my concordant file: is there a way that I can force the addition of the DP field into the Haplotype Caller vcf? I have the AD field showing the read numbers for both alleles, and the DP field is included in the INFO column but I can't filter depth at present so I'd really appreciate some guidance with this. PS. I'm a newbie to bioinformatics processing and computing. Thanks for any help given.

Tagged:

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Rachel,

    Merging the files using awk is not a good way to achieve what you want. Look into the CombineVariants tool instead.

    Coverage annotation can be added to a VCF using VariantAnnotator if you have the original bam files available.

  • rachelwiltshirerachelwiltshire IN, USAMember

    Thanks Geraldine. I'll give this a go and let you know whether it solves my issue.

  • rachelwiltshirerachelwiltshire IN, USAMember

    Hi Geraldine,

    I tried VariantAnnotator with my .bam file but the DP field is still not being emitted in the format (as it is with Unified Genotyper using the same .bam file). I also saw on the forum that only v2.6 and above allows for this so I downloaded v3.4 and reran HaplotypeCaller. My output remains the same: GT:AD:GQ:PL. Do you have any other ideas please?

    Thanks.

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @rachelwiltshire
    Hi,

    Can you tell me if all of the records in the VCF do not have DP, or if only some of the records do not have DP? Have you tried running HaplotypeCaller with the -A Coverage option? What is the exact command you ran with Haplotype Caller?

    Thanks,
    Sheila

  • rachelwiltshirerachelwiltshire IN, USAMember

    HI Sheila,

    It is every line in the vcf. There is no DP field just GT:AD:GQ:PL. This is the call I ran:

    java -Xmx12G -jar ~/GenomeAnalysisTK.jar \
    -T VariantAnnotator \
    -R ~/AGam_scaffolds/AGam_PEST_Scaffolds_ref.fa \
    --fix_misencoded_quality_scores \
    -I ~/AGam_scaffolds/aln_sorted.bam \
    --variant ~/AGam_scaffolds/AGam_HapCall_no_DP/AGam_HapCall_merged.vcf \
    -A AlleleBalance -A BaseCounts -A BaseQualityRankSumTest -A ChromosomeCounts \
    -A Coverage -A FisherStrand -A GCContent -A HaplotypeScore \
    -A HomopolymerRun -A InbreedingCoeff -A LowMQ -A MappingQualityRankSumTest \
    -A MappingQualityZero -A NBaseCount -A QualByDepth -A RMSMappingQuality \
    -A ReadPosRankSumTest -A SpanningDeletions -A TandemRepeatAnnotator \
    -o AGam_HapCall_VariantAnnotator.vcf \
    -rf BadCigar \
    -nt 8

    Thanks.

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @rachelwiltshire
    Hi,

    What happens if you run Haplotype Caller with the latest version of GATK? Haplotype Caller should by default output the DP. You can also try adding -A Coverage to the Haplotype Caller command.

    Thanks,
    Sheila

  • rachelwiltshirerachelwiltshire IN, USAMember

    This is part of my output from HaplotypeCaller in GATK v3.4-46:

    QRankSum=-1.141;DP=2502;FS=0.000;MLEAC=6;MLEAF=0.032;MQ=59.54;MQRankSum=-0.644;QD=10.06;ReadPosRankSum=-0.954;SOR=0.626 GT:AD:GQ:PL 0/0:42,0:99:0,126,1039 0/0:33,0:98
    =0.000;MLEAC=6;MLEAF=1.00;MQ=53.95;QD=31.92;SOR=3.258 GT:AD:GQ:PL ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./. ./.
    FS=0.000;MLEAC=12;MLEAF=1.00;MQ=48.50;QD=32.41;SOR=4.174 GT:AD:GQ:PL ./. ./. ./. ./. ./. ./. ./. ./. 1/1:0,1:3:45,3,0 ./. ./.
    FS=0.000;MLEAC=12;MLEAF=1.00;MQ=48.50;QD=29.99;SOR=4.174 GT:AD:GQ:PL ./. ./. ./. ./. ./. ./. ./. ./. 1/1:0,1:3:45,3,0 ./. ./.
    QRankSum=6.435;DP=908;FS=4.883;MLEAC=9;MLEAF=0.048;MQ=58.93;MQRankSum=2.042;QD=18.49;ReadPosRankSum=4.607;SOR=4.580 GT:AD:GQ:PL 0/0:4,0:12:0,12,180 0/0:10,0:30:0,30,44
    .376e-03;AN=186;BaseQRankSum=7.936;DP=1295;FS=0.000;MLEAC=1;MLEAF=5.376e-03;MQ=59.16;MQRankSum=-8.045;QD=19.98;ReadPosRankSum=0.425;SOR=0.677 GT:AD:GQ:PL 0/0:7,0:21:0,21,315
    QRankSum=1.145;DP=656;FS=0.000;MLEAC=8;MLEAF=0.044;MQ=58.74;MQRankSum=4.080;QD=11.27;ReadPosRankSum=-2.721;SOR=0.083 GT:AD:GQ:PL 0/0:3,0:9:0,9,135 0/0:9,0:30:0,30,450
    QRankSum=-0.608;DP=653;FS=0.000;MLEAC=8;MLEAF=0.044;MQ=58.88;MQRankSum=1.670;QD=11.27;ReadPosRankSum=4.941;SOR=0.082 GT:AD:GQ:PL 0/0:3,0:9:0,9,135 0/0:10,0:30:0,30,45
    QRankSum=2.949;DP=653;FS=0.000;MLEAC=7;MLEAF=0.038;MQ=58.88;MQRankSum=3.588;QD=15.69;ReadPosRankSum=4.754;SOR=0.232 GT:AD:GQ:PL 0/0:3,0:9:0,9,128 1/1:0,10:30:404,30,

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @rachelwiltshire
    Hi,

    Can you please post the command you ran with Haplotype Caller? Is this the output with -A Coverage?

    Thanks,
    Sheila

  • rachelwiltshirerachelwiltshire IN, USAMember

    java -Xmx2g -jar ~/AGam_scaffolds/AGam_HapCall_v3.4/GenomeAnalysisTK.jar -T HaplotypeCaller -R ~/AGam_scaffolds/AGam_PEST_Scaffolds_ref.fa --downsampling_type NONE --annotation AlleleBalance -rf BadCigar -nct 8 --fix_misencoded_quality_scores -I ~/AGam_scaffolds/aln_sorted.bam -o AGam_HapCall_v3.4.vcf

    This is NOT with -A coverage. That job is currently running and I'll let you know if DP emits in the format field if so.

  • rachelwiltshirerachelwiltshire IN, USAMember

    Hello again.

    So I finished running HaplotypeCaller v3.4 with the flag -A coverage: my output remains the same; there is no DP field emitted. Here is some output from the file:
    20.00;MQRankSum=-0.359;QD=1.37;ReadPosRankSum=-0.399;SOR=0.000 GT:AD:GQ:PL 0/0:63,0:99:0,190,2959 0/0:25,0:75:0,75,1196 0/0:56,0:99:0,169,2680
    20.00;MQRankSum=0.553;QD=1.37;ReadPosRankSum=-0.021;SOR=0.000 GT:AD:GQ:PL 0/0:63,0:99:0,190,2959 0/0:25,0:75:0,75,1196 0/0:56,0:99:0,169,2680
    ;MQRankSum=0.728;QD=2.68;ReadPosRankSum=0.449;SOR=0.030 GT:AD:GQ:PL 0/1:43,20:99:306,0,1178 0/0:25,0:75:0,75,760 0/1:45,11:99:104,0,1127 0/1:7,5:
    MQRankSum=0.274;QD=1.72;ReadPosRankSum=-3.304;SOR=0.002 GT:AD:GQ:PL 0/0:61,2:99:0,140,2128 0/0:22,3:1:0,1,591 0/0:51,5:44:0,44,1442 0/0:12,0
    MQRankSum=1.852;QD=1.59;ReadPosRankSum=-1.075;SOR=0.014 GT:AD:GQ:PL 0/1:50,13:99:132,0,1493 0/1:21,4:24:24,0,546 0/0:52,4:69:0,69,1511 0/1:8,4:
    20.00;MQRankSum=0.012;QD=0.65;ReadPosRankSum=-2.939;SOR=0.000 GT:AD:GQ:PL 0/1:52,11:83:83,0,1358 0/0:25,0:75:0,75,699 0/0:56,0:99:0,168,1466

    This was the command I used:
    java -Xmx2g -jar ~/AGam_scaffolds/AGam_HapCall_v3.4/GenomeAnalysisTK.jar -T HaplotypeCaller -R ~/AGam_scaffolds/AGam_PEST_Scaffolds_ref.fa --downsampling_type NONE --annotation AlleleBalance -rf BadCigar -nct 8 --fix_misencoded_quality_scores -A Coverage -I ~/AGam_scaffolds/aln_sorted.bam -o AGam_HapCall_v3.4_-ACoverage.vcf

    So, to recap: none of my jobs using HaplotypeCaller v3.4-46 are producing output with the DP field whether it is a) run alone, b) run with -T VariantAnnotator or c) run with -A Coverage.

    Hope you can help.

    PS. UnifiedGenotyper works fine.

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @rachelwiltshire
    Hi,

    It seems that adding -A "annotation" turns off the default annotations output in the VCF. So, for example, if you don't add any -A "annotation" to your command, DP shows up.

    However, you can work around this by using -G Coverage -G AlleleBalancePerSample -G "otherAnnotations" https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php#--group

    I am going to put in a ticket to make this behavior more straightforward.

    -Sheila

    Issue · Github
    by Sheila

    Issue Number
    1126
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @rachelwiltshire
    Hi,

    It turns out this is a known issue, and we are working on fixing it.

    -Sheila

  • rachelwiltshirerachelwiltshire IN, USAMember

    Thanks Sheila. I look forward to hearing about the fix. In the meantime, I'll continue with Unified Genotyper..

Sign In or Register to comment.