Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. 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.

A snip and deletion without DP in g.vcf file

I have noticed that some times g.vcf files will have calls where the DP field is completely missing from the FORMAT string, even though a variant is called with good GQ. When this happens, the INFO string does have a DP=0, as can be seen below.

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  102637-001-093
gi|194447306|ref|NC_011083.1|   415537  .   A   <NON_REF>   .   .   END=415547  GT:DP:GQ:MIN_DP:PL  0:2:86:2:0,86
gi|194447306|ref|NC_011083.1|   415548  .   G   A,<NON_REF> 13.22   .   DP=0;MLEAC=1,0;MLEAF=1.00,0.00;RAW_MQ=0.00  GT:GQ:PL:SB 1:43:43,0,43:0,0,0,0
gi|194447306|ref|NC_011083.1|   415549  .   A   <NON_REF>   .   .   END=415549  GT:DP:GQ:MIN_DP:PL  0:2:86:2:0,86

gi|194447306|ref|NC_011083.1|   4683672 .   A   <NON_REF>   .   .   END=4683672 GT:DP:GQ:MIN_DP:PL  0:2:0:2:0,0
gi|194447306|ref|NC_011083.1|   4683673 .   AATC    A,<NON_REF> 6.95    .   DP=0;MLEAC=1,0;MLEAF=1.00,0.00;RAW_MQ=0.00  GT:GQ:PL:SB 1:45:45,0,45:0,0,0,0
gi|194447306|ref|NC_011083.1|   4683677 .   A   <NON_REF>   .   .   END=4683677 GT:DP:GQ:MIN_DP:PL  0:1:0:1:0,0

I'm using gatk version 3.5-0-g36282e4, and I've attached the files I used. The exact command I used is as folows:

bamfile=093_interval.bam
output=093_interval.g.vcf
bamout=093_interval_bamout.bam
reference=NC_011083.fasta

gatk -T HaplotypeCaller \
            --sample_ploidy 1 \
            -R $reference \
            -I $bamfile \
            -o $output \
            -ERC GVCF \
           -bamout $bamout 

How can I make sense of this, to my mind DP=0 (from the INFO field at least) should mean there are no reads, so therefore no call can be made, right?

Tagged:

Issue · Github
by Sheila

Issue Number
561
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
chandrans

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Redmar_van_den_Berg
    Hi,

    I have been trying to figure out what is going on, but I can't! There are a few weird things going on in the bamout file and graph out file. Let me check with the team about what to do next. I will get back to you.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Redmar_van_den_Berg
    Hi again,

    I should note, after GenotypeGVCFs, those two odd sites are not output in the final VCF.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Looks like a bug to me.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Redmar_van_den_Berg
    Hi,

    Can you tell us how many sites you see this issue with? Are there a lot of sites or just a handful? Will most of them fail filters downstream?

    Thanks,
    Sheila

  • With the latest gatk, this happens to a low number of sites. I've just re-ran a sample with the latest gatk and checked for missing DP. 093 is a salmonella sample (genome size ~ 4.5M) with 50K snips, and I only found 2 occurrences of missing DP in the .g.vcf file.

    I don't recall exactly how many instances I found when I originally reported this, I can re-run this analysis with gatk 3.5 if you need to know that figure as well.

    $ gvcf2fasta.py 093_sorted_dedup_realn.g.vcf --reference /usr/share/ngs/Reference/NC_011083.fasta
    Parsing 093_sorted_dedup_realn.g.vcf
    Warning, no DP for sample '093' in Record(CHROM=gi|194447306|ref|NC_011083.1|, POS=415548, REF=G, ALT=[A, <NON_REF>])
    Warning, no DP for sample '093' in Record(CHROM=gi|194447306|ref|NC_011083.1|, POS=4683673, REF=AATC, ALT=[A, <NON_REF>])
    
    $ grep "415548\|4683673" 093_sorted_dedup_realn.g.vcf
    gi|194447306|ref|NC_011083.1|   415548  .   G   A,<NON_REF> 13.22   .   DP=0;MLEAC=1,0;MLEAF=1.00,0.00;RAW_MQ=0.00  GT:GQ:PL:SB 1:43:43,0,43:0,0,0,0
    gi|194447306|ref|NC_011083.1|   4683673 .   AATC    A,<NON_REF> 6.95    .   DP=0;MLEAC=1,0;MLEAF=1.00,0.00;RAW_MQ=0.00  GT:GQ:PL:SB 1:45:45,0,45:0,0,0,0
    
    $ gatk --version
    3.6-0-g89b7209
    
    # Final vcf file
    $ wc -l ../snps/joint_snips.vcf 
    49772 ../snps/joint_snips.vcf
    
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Redmar_van_den_Berg
    Hi,

    No need to re-run anything. We are just trying to get a sense of how many sites this issue occurs at. A few other users have reported the same issue, but all have agreed it occurs only at a few sites.

    Thanks,
    Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Redmar_van_den_Berg
    Hi,

    This has been fixed and is available in the nightly builds.

    -Sheila

Sign In or Register to comment.