<p>With the current versions of samtools, you can also use:</p>
<pre><code>samtools coverage sample.bam
</code></pre>
<p>The documents are <a href="https://www.htslib.org/doc/samtools-coverage.html" rel="nofollow noreferrer">here</a>,</p>
<blockquote>
<p>samtools coverage – produces a histogram or table of coverage per
chromosome</p>
</blockquote>
<p>A further example from the site,</p>
<pre><code>samtools coverage -r chr1:1M-12M input.bam
</code></pre>
<p><strong>Output</strong></p>
<pre><code>#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8
</code></pre>
<p>The histogram option is</p>
<pre><code>samtools coverage -A -w 32 -r chr1:1M-12M input.bam
</code></pre>
<p>... and output available on the website link.</p>

CC BY-SA 4.0
2022-10-12T11:38:51.013Z
2022-10-12T15:30:16.433Z
2022-10-12T15:30:16.433Z
pdp10
<p>For a quick estimate you’re making it more complicated than necessary.</p>
<ol>
<li><p>The theoretical average coverage is <span class="math-container">$\frac{n \cdot \hat l}{N}$</span> where <span class="math-container">$n$</span> is the number of reads, <span class="math-container">$\hat l$</span> is the average read length and <span class="math-container">$N$</span> is the genome size.</p>
</li>
<li><p><code>samtools idxstats</code> gives you the chromosome lengths and number of mapped reads in one convenient list.</p>
</li>
</ol>
<p>Putting this together, we get</p>
<pre><code>samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'
</code></pre>
<p>(In this case I’ve used 120 as the read length.)</p>
<p>However, note that this over-estimates <em>empirical</em> sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average <em>mapped</em> length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.</p>
<hr />
<p>If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the <code>-a</code> argument to <code>samtools depth</code>: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:</p>
<pre><code>sum=$(samtools depth &quot;$1&quot; | awk '{sum+=$3} END {print sum}')
</code></pre>
<p>Some more remarks on your code:</p>
<ul>
<li>The <code>-P</code> flag for <code>grep</code> is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: <code>grep '^@SQ'</code> works just as well.</li>
<li>While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:
<ul>
<li>You can often combine <code>cut</code> and <code>awk</code>, in this case by setting the field separator:
<pre><code>awk -vFS=: '{sum+=$3} END {print sum}'
</code></pre>
</li>
<li>You can often combine <code>grep</code> and <code>awk</code>:
<pre><code>awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'
</code></pre>
</li>
</ul>
</li>
<li><em>Always</em> double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters <em>are valid characters in paths!</em></li>
<li>Your code is clean and readable. Well done.</li>
</ul>

CC BY-SA 4.0
2020-07-04T13:26:23.603Z
2020-07-06T08:30:06.137Z
2020-07-06T08:30:06.137Z
Konrad Rudolph
<p><em>(Full discolosure that this is my first time working with sequence data, and with the bash scripting.)</em></p>
<p>I need to calculate the average coverage for any .bam file.</p>
<p>After some searching I wrote the following script:</p>
<pre class="lang-bsh prettyprint-override"><code># Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=<span class="math-container">$(samtools view -H $</span>1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=<span class="math-container">$1} END {print sum}')
echo $</span>tot

#compute the coverage at each point
samtools depth -ao cov_<span class="math-container">$1 $</span>1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=<span class="math-container">$(awk '{sum+=$3} END {print sum}' cov_$</span>1)
echo $sum
echo

avg=<span class="math-container">$(echo "$</span>sum/<span class="math-container">$tot" | bc -l)
echo ''The average coverage of the sample $</span>1 is $avg x.''
</code></pre>
<p>The logic behind the calculation of tot is as follows:</p>
<p><code>samtools view -H </code> gets the header of the bam file</p>
<p><code>grep -P</code> gets all lines in this header with SQ</p>
<p><code>cut -f 3 -d ':'</code> removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes</p>
<p><code>awk '{sum+=$1} END {tot=sum}'</code> calculates the length/total bases in the chromosome</p>
<p>However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.</p>
<p><em>One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...</em></p>

<sequence-analysis><bash><coverage>
Calculating average coverage for .bam files (sequence data)
# Calculating average coverage for .bam files (sequence data)

(Asked by: username_0 on 2020-07-04)

(Full discolosure that this is my first time working with sequence data, and with the bash scripting.)
I need to calculate the average coverage for any .bam file.
After some searching I wrote the following script:
# Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=$(samtools view -H $1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}')
echo $tot

#compute the coverage at each point
samtools depth -ao cov_$1 $1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=$(awk '{sum+=$3} END {print sum}' cov_$1)
echo $sum
echo

avg=$(echo "$sum/$tot" | bc -l)
echo ''The average coverage of the sample $1 is $avg x.''

The logic behind the calculation of tot is as follows:
samtools view -H  gets the header of the bam file
grep -P gets all lines in this header with SQ
cut -f 3 -d ':' removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes
awk '{sum+=$1} END {tot=sum}' calculates the length/total bases in the chromosome
However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.
One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...

**Tags:** sequence-analysis, bash, coverage

---

## 2 answers

### Answer by username_2 on 2020-07-04. Score: 9 | Accepted answer

For a quick estimate you’re making it more complicated than necessary.

The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.

samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list.


Putting this together, we get
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

(In this case I’ve used 120 as the read length.)
However, note that this over-estimates empirical sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average mapped length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the -a argument to samtools depth: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

Some more remarks on your code:

The -P flag for grep is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: grep '^@SQ' works just as well.
While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:

You can often combine cut and awk, in this case by setting the field separator:
awk -vFS=: '{sum+=$3} END {print sum}'


You can often combine grep and awk:
awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'




Always double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters are valid characters in paths!
Your code is clean and readable. Well done.


---

### Answer by username_1 on 2022-10-12. Score: 2

With the current versions of samtools, you can also use:
samtools coverage sample.bam

The documents are here,

samtools coverage – produces a histogram or table of coverage per
chromosome

A further example from the site,
samtools coverage -r chr1:1M-12M input.bam

Output
#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

The histogram option is
samtools coverage -A -w 32 -r chr1:1M-12M input.bam

... and output available on the website link.


---

