Calculate depth of every base position

[Lv Ray]

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]

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

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

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]

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]

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

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

[Brian Bushnell]

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]

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