Calculating coverage based on full insert size from paired-end data

[rnaeye]

Hi! all,

I would like to calculate coverage from paired-end data. I can do this with bedtools but I want to make sure of one thing. That is to include full length insert in the coverage. For instance, I have 2 x 100 bp paired-end sequencing data and my average insert size is 500 bp. When I calculate the coverage I want to make sure that I include unsequenced part between two reads.

Read1       unsequenced part                Read2
--------->|-------------------------------|<-----------
100 nt    |         300 nt                |   100 nt

How can I make sure that coverage calculation includes unsequenced 300 nt part in a proper pair. Also, Is there a way of calculating coverage for both properly paired reads and singletons at the same time. Any tool is ok, Perl script, bedtools, samtools, bamtools etc.

Thank you for suggestions in advance.

[pli]

Prolly way to late, but it might help others, here's a script I wrote recently to do just that. Note that it produces a SAM file with all sequences just N's, and a CIGAR string of just M's.

https://gist.github.com/plijnzaad/592145e2a57dfbd6ac07

Enjoy,

Philip

[rnaeye]

Thanks Philip. It's not too late. It's good to check the same thing with different approaches. I wil give a try. Best,

[Brian Bushnell]

BBMap has a "physcov" flag that allows it to report physical rather than sequenced coverage. It can be used directly in BBMap, or with pileup, if you already have a sam file. For example:

pileup.sh in=mapped.sam covstats=coverage.txt