samtools depth -a while removing duplicates
2
0
Entering edit mode
8.4 years ago

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!

samtools depth • 3.6k views
ADD COMMENT
1
Entering edit mode
8.4 years ago

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.

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

depth() is essentially a wrapper around mpileup, so just use the latter.

ADD REPLY
0
Entering edit mode
8.4 years ago
winter_li ▴ 60

I think if you wanna remove duplicates reads ,you must mark duplicates first. I recommend you to use Picard to do this.

ADD COMMENT

Login before adding your answer.

Traffic: 2710 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6