Question: How To Calculated Coverage Of The Genome/Contig Slice?
5
gravatar for Leszek
8.4 years ago by
Leszek4.0k
IIMCB, Poland
Leszek4.0k wrote:

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?

ADD COMMENTlink modified 8.4 years ago by lh331k • written 8.4 years ago by Leszek4.0k

possible duplicate? http://biostar.stackexchange.com/questions/5181/tools-to-calculate-average-coverage-for-a-bam-file

ADD REPLYlink written 8.4 years ago by Casey Bergman18k

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

ADD REPLYlink written 8.4 years ago by Leszek4.0k
11
gravatar for lh3
8.4 years ago by
lh331k
United States
lh331k wrote:

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.

ADD COMMENTlink written 8.4 years ago by lh331k
7
gravatar for Alex Stoddard
8.4 years ago by
Alex Stoddard190
Wisconsin, USA
Alex Stoddard190 wrote:

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

ADD COMMENTlink written 8.4 years ago by Alex Stoddard190
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1121 users visited in the last hour