Samtools mean depth and max depth results are discordant

Question

[rachel.shadoff]

I have used the samtools coverage function to compute mean depth of coverage for a bam file, and the samtools depth function in conjunction with cat to determine the maximum and minimum depth of coverage values for this same bam. Oddly, the maximum depth of coverage value from the txt created using samtools depth is a smaller value than the mean depth of coverage given by samtools coverage.

This is the code I am applying to the BAM file:

samtools coverage 3_sorted.bam

This returns a value of 15,177 as the mean depth of coverage.

samtools depth 3_sorted.bam -o 3depth.txt    
cat 3depth.txt | sort -nk3,3 | head -1

This returns a value of 8040 as the maximum depth of coverage (maximum value in the coverage column of the txt file)

I am hoping someone can point out the mistake in my code, or can explain why the max depth value is smaller than the average (mathematically impossible). Thanks!

[raphael.B]

Isn't mean depth 15.177 instead of 15177? That number seems absolutely enormous...

[rachel.shadoff]

Hi Raphael, I belive 15177 is correct. When I run

samtools coverage 3_sorted.bam

the result is:

#rname    startpos  endpos  numreads    covbases    coverage    meandepth   meanbaseq   meanmapq
NC_012920.1 1   16569   1947653         16569           100          15177.4    37.5             60

[raphael.B]

This is huge. Just to know, what interest do you have in getting such a sequencing depth?

Regarding your first question, my guess is that samtools depth does not provide you the depth of all positions. Try the -aa option.

Answer 1

[jkbonfield]

Older versions of samtools depth had a very dumb depth sub-command with a maximum depth limit (approx 8000). This was rather unexpected behaviour! Some versions support a -d option to adjust this maximum depth, although in the very earliest releases that wasn't adjustable.

The depth command was completely rewritten in 1.13 and no longer has this odd limitation. It's also now significantly faster. On an excessive ~130K depth covid-19 sample the new code was 17x quicker!

I'd advise upgrading to the latest release.

[rachel.shadoff]

Using the 1.13 version fixed the issue, and I'm now seeing depths well over 8000. Thanks so much!!!