I want to run samtools depth with the -a flag to capture positions with 0 coverage.
However due to the circumstances I do not have bam files with the duplicates removed. I realize I can run rmdups but we chose not to due to financial and time constraints (analyzing 1000s of whole genomes on cloud service)
So I would like to ensure I'm not counting duplicate reads with the samtools depth command. Is it possible to use a flag for this?
I would also like to include the -a flag because positions with 0 coverage are informative for my analysis.
Thanks!
•
written
by
QVINTVS_FABIVS_MAXIMVS
There's no way to do this directly with samtools depth. If you really need this, then you can do it in either pysam (this is easier) or htslib (assuming you know C). In either case you'll need to write the duplicate removal code.
I'm actually doing this in pysam but there is no depth function for Aligned_File class (I think). So I'm using pysam.depth()
I suppose you are correct. I could fetch the reads within a region, remove dups and hard code some per base pair depth function.
depth() is essentially a wrapper around mpileup, so just use the latter.
I think if you wanna remove duplicates reads ,you must mark duplicates first. I recommend you to use Picard to do this.
Log in to answer this question.
