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.

DP=0 in VCF

Hi,
i've run the HaplotypeCaller on single sample, and got a few variants in the VCF file that had only 3 fields in FORMAT:
GT:GQ:PL (instead of five: GT:AD:DP:GQ:PL). when i look at the the INFO field, i see that their "DP=0", and indeed there was no coverage at their positions (by looking at teh BAM file in IGV).
if the DP=0, then why and how these variants were called?
what can i learn about them?

Thanks for your help
Moran

Answers

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    Can you please post the vcf record at the site and the bamout IGV screenshot?

    Thanks,
    Sheila

  • galmorangalmoran Member

    Thanks for your help
    here are two examples:

    example 1:
    chr7 107336463 . C T 457.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=5.440;ClippingRankSum=0.320;DP=126;FS=0.000;MLEAC=1;MLEAF=0.500;MQ

    example 2:
    chr8 2131251 . A C 332.78 . AC=2;AF=1.00;AN=2;DP=0;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=0.00;MQ0=0 GT:GQ:PL 1/1:25:361,25,0

    the IGV screenshot are attached

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    It looks like the IGV screenshot of the first example does not match the example vcf record you posted. Can you please confirm that these IGV screenshots are of the bamout file and not of the original bam file?

    Thanks,
    Sheila

  • galmorangalmoran Member

    sorry, the first example correspond to this record:
    chr7 71912741 . TGTACACAC T 99.14 . AC=2;AF=1.00;AN=2;DP=0;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=0.00;MQ0=0 GT:GQ:PL 1/1:10:136,10,0

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    Thanks. Haplotype Caller does a local reassembly step that may shift reads from their original positions in the bam file. I think that is what is happening here. The reads after the local reassembly step can be viewed by using the bamout argument. https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php#--bamOutput If you can post the IGV screenshot after the reassembly step, that would be very helpful.

    -Sheila

  • galmorangalmoran Member

    Hi Sheila, i still don't understand. shoudl i re-run the Haplotype Caller on the bams with bamout argument? i see that it would cost lots of time and not recommended. did you mean i'll turn te bamout flag on only for those variants? if so, how do i do that?
    does the HC local assembly step is done in adiition to the local re-alignment step ? ( i used the recoomended GATK pipeline taht included oacl re-alignment and base recalibrator steps).

    but the most important thing - should i treat this variants are true calls? and if so- i would need their DP and AD info.
    Thanks for your help
    Moran

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    I just need to see an IGV screenshot of the bamout file for the region you posted above. You can simply run Haplotype Caller on that region. Please include ~500 bases before and after the position of interest when you run Haplotype Caller.

    Thanks,
    Sheila

  • galmorangalmoran Member

    Hi
    i re-run the HaplotypeCaller again in the intervals of the 2 examples listed above with the bamout argument. these are the IGV screenshots i got from thr bamout file
    Hope it will help
    Many thanks
    Moran

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    Can you confirm you are using GATK's latest version (3.3)? It looks like the AD is reflecting pre-reassembly values. This has been fixed in the latest version.

    After reassembly, Haplotype Caller is calling a deletion at the site. You can look at the PLs to see how much confidence there is in the variant call. Is this site filtered out after VQSR?

    -Sheila

  • galmorangalmoran Member

    Hi Sheila,
    i am using 3.1.1. version of GATK, and i have not used VQSR filtering in this project.
    where did you see the AD values? there were not omitted in the VCF records. .
    without getting too much into the details of GATK algorithm, can you say that variants that are omitted this way (with apparent DP=0 and no AD info) are less reliable? would you trust them?
    i am not interested in these specific variants, but asked about them because i am calculating a general statsistics about the DP and AD (which reflect the ratio of teh mutant allele) of the VarCall in this project. where can i retrieve the DP and AD for them?

    Thanks again for your help
    Moran

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    @galmoran If you are planning to do anything with AD you should redo the analysis with the latest version because that bug produced unreliable AD values. The calls themselves were good but you can't trust the AD annotation.

  • galmorangalmoran Member

    I see, but what it their DP? it cannt be DP=0 (As written in the info field) ...?

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @galmoran
    Hi Moran,

    If you upgrade to the latest version, the DP will be fixed, along with the AD.

    -Sheila

  • jlrfloresjlrflores ✭✭ Member ✭✭

    I just observed this same issue in 3.4, is it fixed in 3.5?
    I noticed the problem because I was running a script that pulls the depth for each variant (typically 3rd column in the formatted genotype string). Any reason why these varants only have 3 columns (GT:GQ:PL instead of GT:AD:DP:GQ:PL).

    20 1592383 rs139926854 A ATT 53.70 . AC=2;AF=1.00;AN=2;DB;DP=0;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=0.00;SOR=0.693 GT:GQ:PL 1/1:6:90,6,0

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    No idea -- what was the command line? Are all calls affected or just some?

    Fastest way to know if the problem will remain in 3.5 is to try it on a snippet...

Sign In or Register to comment.