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.