How To Calculated Coverage Of The Genome/Contig Slice?

Question

[Leszek]

Hi,

I have mapped reads (SAM/BAM) from resequencing and I need an elegant way of calculating depth of coverage for a slices of assembly (let's say 100bp). It's possible to use samtools:

samtools view -c reads.bam contig1.142:2000-1000

but this just count number of reads that are completely aligned to given fragment and a strand (FLAG: 0 - forward 16 - reverse) if I'm right. Do you know any better way of doing it?

[Casey Bergman]

possible duplicate? Tools To Calculate Average Coverage For A Bam File?

[Leszek]

I need coverage of a small fragments, not the depth of coverage for the genome as whole

Answer 1

[lh3]

If you only want the coverage in a couple of regions:

samtools mpileup -ABr chr1:100-200 -d 1000000 in.bam

If you want to get the depth in multiple regions, but not so many (less than 100,000):

awk '{print $1":"($2+1)"-"$3}' in.bed | xargs -i samtools mpileup -ABr {} -d 1000000 in.bam > output

If you want to get the depth in many regions, use BEDTools:

coverageBed -abam in.bam -b in.bed

Which method to use also depends on the total length of regions. If the total length is above 1% of the genome and everywhere, BEDTools is probably more efficient.

As to your question, samtools view retrieves alignment overlapping the specified region, not only contained. It outputs alignments on both strands, unless you use -f/-F to control the output.

Answer 2

[Alex Stoddard]

I have used BEDTools for this purpose.

Despite the name BEDTools works with a whole range of common input formats including BAM.

See the manual for coverageBed. You can use a bed file to specify your region(s) of interest and get coverage details and summaries based on reads in a BAM format file.

e.g.

$ coverageBed -abam yourBAMfile.bam -b regionsOfInterest.bed