One of the most common question after your Next-Generation Sequencing (NGS) run is: “What is the coverage I achieved?” Unfortunately there a different definitions available for the term “coverage”, but in general I would suggest to use “How many different reads are aligned at any position of the target genome on average?”
To answer this question the straightforward calculation is called Lander-Waterman equation and states:
C = L * N / G
meaning the coverage C is the length of the reads (L) multiplied by the number of reads (N) divided by the length of the target genome (G).
This should be done after duplicates and non-aligning reads have been removed. From a BAM file you can use samtools and Unix awk magic to get an average with standard deviation:
samtools depth -a aligned_reads.bam | awk '{sum+=$3; sumsq+=$3*$3} END { print "Average = ",sum/NR; print "Stdev = ",sqrt(sumsq/NR - (sum/NR)**2) }'
It is however important to remember that this will still only be an estimate of the average coverage of each base of the genome. There are regions in most genomes that are harder to generate sequencing reads for and regions that are harder to get a reliable alignment for (Mappability). This means that even though you might have the average coverage of 30x required for many applications for the analysis of the (human) genome, there will be regions where there are no reads at all! To assess this kind of coverage the term “breadth of coverage” has been established to describe how many regions of the genome are covered by reads.
Further details, in particular for Illumina NGS runs and how to calculate the target depth can be found at the Illumina site.
A tools discussion at Biostars.
