This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to get sequencing depths from VCF with Rsamtools

I got from our collaborator some pindel output and the BAM files that is used as pindel input. Unfortunately, they used an older of version of pindel (before approximately 0.2.4u) that doesn't record total depth at breakpoint (i.e. the DP field in VCF). While our collaborator's standard practice is to obtain depths through Rsamtools, I haven't found a good way to do this. Specifically, I know how to read sequencing depths to an Rle objects, and on the other end convert pindel outputs to VCF and then to GenomicRanges, but how do I integrate the two?

r pindel rsamtools
• 3,755 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by John Ma John Ma

"By design, the call annotations in the VCF reflect the data that the caller/genotyper actually used in their calculation. " how could you known the number of reads used by pindel ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

The docs for the current versions of pindel states depth is the same I'd get from samtools depth.

• 0 views • link
written by John Ma John Ma

1 answer

I've written a tools to update the DP fields from a set of BAM files: http://lindenb.github.io/jvarkit/FixVcfMissingGenotypes.html

$ find ~/src/gatk-ui/testdata/ -name "*.bam" > input.list

$ tail -2 input.vcf
rotavirus   1064    .   G   A   21.5606 .   DP=250;VDB=2.70971e-16;SGB=8.40135;RPB=0.935144;MQB=1;BQB=0.683886;MQ0F=0;AF1=0.25;G3=0.75,2.37734e-17,0.25;HWE=0.033921;AC1=2;DP4=0,219,0,31;MQ=60;FQ=22.8019;PV4=1,1.22605e-06,1,1    GT:PL   0/0:0,244,70    0/0:0,199,65    0/0:0,217,68    1/1:69,84,0
rotavirus   1064    .   G   A   21.5606 .   DP=250;VDB=2.70971e-16;SGB=8.40135;RPB=0.935144;MQB=1;BQB=0.683886;MQ0F=0;AF1=0.25;G3=0.75,2.37734e-17,0.25;HWE=0.033921;AC1=2;DP4=0,219,0,31;MQ=60;FQ=22.8019;PV4=1,1.22605e-06,1,1    GT:PL   ./. ./. ./. ./.

$ java -jar dist/fixvcfmissinggenotypes.jar -d 50 --fixDP --filtered zz -B input.list input.vcf | tail -2
rotavirus   1064    .   G   A   21.56   .   AC1=2;AF1=0.25;BQB=0.683886;DP=188;DP4=0,219,0,31;FQ=22.8019;G3=0.75,2.37734e-17,0.25;HWE=0.033921;MQ=60;MQ0F=0;MQB=1;PV4=1,1.22605e-06,1,1;RPB=0.935144;SGB=8.40135;VDB=2.70971e-16    GT:DP:PL    0/0:48:0,244,70 0/0:63:0,199,65 0/0:53:0,217,68 1/1:24:69,84,0
rotavirus   1064    .   G   A   21.56   .   AC1=2;AF1=0.25;BQB=0.683886;DP=72;DP4=0,219,0,31;FQ=22.8019;G3=0.75,2.37734e-17,0.25;HWE=0.033921;MQ=60;MQ0F=0;MQB=1;PV4=1,1.22605e-06,1,1;RPB=0.935144;SGB=8.40135;VDB=2.70971e-16 GT:DP:FT:FXG    ./.:48:PASS 0/0:63:zz:1 0/0:53:zz:1 ./.:24:PASS
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this