Calculating average coverage for .bam files (sequence data)

Question

(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...

Answer 1

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

1. 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.

2. 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 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.