Calculate depth of every base position

Lv Ray · 10-16-2014, 01:00 AM

I have the sorted bam file, and I want to calculate the depth of every position on genome including some positions uncovered by reads(marked as: chr pos 0). I have tried samtools depth, but it seems don't count some position(because I count the row number, and not equal to the total number of genome bases).
Any advice about this would be much appreciated, thank you.

blakeoft · 10-16-2014, 01:21 PM

I use bedtools genomecov. There are several output options, but I prefer to run

Code:

bedtools genomecov -bg -ibam your_bam.bam > your_bam.bedgraph

The fields in the output are chr, start, end and then depth. For example, you might see something like

Code:

chr1  554304  554309  5
chr1  554309  554313  6

This means that the depth is 5 from 554304 to 554309 and then the depth is 6 from 554309 to 554313 (I'm pretty sure it's not inclusive on the right endpoint) on chromosome 1.

There are some other options that I haven't actually explored myself, including adding "-d" to the code above that seems to output the depth at every single position, but this would make the file much larger than it would need to be and still contain exactly the same information that you'd have without "-d".

Read more about genomecov here.

Edit: samtools depth might do the same thing as I described above. I don't know since I haven't actually used it. You might want examine its output and see if it's presented in (what I like to call) the RLE way, like genomecov does with just the "-bg" parameter.

WhatsOEver · 10-17-2014, 02:56 AM

Just a small addition:

You need to call bedtools genomcov with -bga to get regions with 0 coverage. The rest blakeoft explained is correct.

blakeoft · 10-17-2014, 05:35 AM

Oops. I misread what Lv Ray posted. I thought that "marked as: chr pos 0" meant that you wanted to start the indexing at zero rather than one. Yes, WhatsOEver is correct. Thanks for clarifying that by the way.

For the record, this is the line you'll want to execute

Code:

bedtools genomecov -bga -ibam your_bam.bam > your_bam.bedgraph

Brian Bushnell · 10-17-2014, 09:53 AM

BBTools contains a tool called pileup.sh, which can also do this rapidly from an unsorted sam or bam file.

pileup.sh -Xmx16g in=mapped.sam basecov=coverage.txt

It can also generate the coverage while mapping, if desired.

Lv Ray · 10-17-2014, 11:04 PM

Thank you all. I will try the bedtools, and i will let u know what my result. Thanks again.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
how to calculate base-by-base coverage?	shuoguo	Bioinformatics	4	02-21-2014 08:06 PM
how to calculate the mean sequencing depth of RNA-seq	wangli	RNA Sequencing	3	02-19-2014 05:44 PM
calculate the base distribution of a SAM/BAM file	weberma	Bioinformatics	3	12-05-2012 12:22 AM
calculate a base's distance to the read's 3' end	lyz1030	Bioinformatics	2	07-12-2012 06:36 AM
How to calculate proportion of reads with each base at every reference position	lindseyjane	Bioinformatics	7	01-19-2010 08:43 AM

10-16-2014, 01:00 AM	#1
Lv Ray Member Location: GZ,China Join Date: Jun 2014 Posts: 42	Calculate depth of every base position I have the sorted bam file, and I want to calculate the depth of every position on genome including some positions uncovered by reads(marked as: chr pos 0). I have tried samtools depth, but it seems don't count some position(because I count the row number, and not equal to the total number of genome bases). Any advice about this would be much appreciated, thank you.

10-16-2014, 01:21 PM	#2
blakeoft Member Location: Connecticut Join Date: Oct 2013 Posts: 79	I use bedtools genomecov. There are several output options, but I prefer to run Code: bedtools genomecov -bg -ibam your_bam.bam > your_bam.bedgraph The fields in the output are chr, start, end and then depth. For example, you might see something like Code: chr1 554304 554309 5 chr1 554309 554313 6 This means that the depth is 5 from 554304 to 554309 and then the depth is 6 from 554309 to 554313 (I'm pretty sure it's not inclusive on the right endpoint) on chromosome 1. There are some other options that I haven't actually explored myself, including adding "-d" to the code above that seems to output the depth at every single position, but this would make the file much larger than it would need to be and still contain exactly the same information that you'd have without "-d". Read more about genomecov here. Edit: samtools depth might do the same thing as I described above. I don't know since I haven't actually used it. You might want examine its output and see if it's presented in (what I like to call) the RLE way, like genomecov does with just the "-bg" parameter. Last edited by blakeoft; 10-16-2014 at 01:25 PM.

10-17-2014, 05:35 AM	#4
blakeoft Member Location: Connecticut Join Date: Oct 2013 Posts: 79	Oops. I misread what Lv Ray posted. I thought that "marked as: chr pos 0" meant that you wanted to start the indexing at zero rather than one. Yes, WhatsOEver is correct. Thanks for clarifying that by the way. For the record, this is the line you'll want to execute Code: bedtools genomecov -bga -ibam your_bam.bam > your_bam.bedgraph

10-17-2014, 02:56 AM	#3
WhatsOEver Senior Member Location: Germany Join Date: Apr 2012 Posts: 215	Just a small addition: You need to call bedtools genomcov with -bga to get regions with 0 coverage. The rest blakeoft explained is correct.

10-17-2014, 09:53 AM	#5
Brian Bushnell Super Moderator Location: Walnut Creek, CA Join Date: Jan 2014 Posts: 2,707	BBTools contains a tool called pileup.sh, which can also do this rapidly from an unsorted sam or bam file. pileup.sh -Xmx16g in=mapped.sam basecov=coverage.txt It can also generate the coverage while mapping, if desired.

10-17-2014, 11:04 PM	#6
Lv Ray Member Location: GZ,China Join Date: Jun 2014 Posts: 42	Thank you all. I will try the bedtools, and i will let u know what my result. Thanks again.