Sequence coverage (or sequencing depth) refers to the number of reads that include a specific nucleotide of a reference genome. In the screenshot below a small region of the Human genome is shown in a genome browser, where the alignment of a re-sequencing experiment was loaded. Each gray arrow in the main area of the window represents a single read, while the highlighted area is a coverage plot, summarizing for each nucleotide how many times was sequenced in this experiment.
but will stop rendering if the region is >50Kbp.
An extended version of this article was published on Medium.
Visualizing coverage from a BAM file
As shown in the introduction, IGV is a simple and effective tool to inspect an alignment file. You need the reference genome used to perform the alignment (load it via Genome → Load from file…) and the alignment file in BAM format (sorted and indexed).
The major limitation is that you can view up to 50 kbp with IGV, if you zoom out you’ll stop seeing alignments. This makes IGV perfect for detailed inspection of small regions (you’ll see the alignment read by read!) but less useful for broader inspection.
Typical case: which regions have a very low (or no) coverage? Which ones have a very high coverage?
Extracting coverage information with samtools
There is a samtools subprogram, called depth, that calculates the sequence coverage at each position³:
samtools depth -a FILE.bam > FILE.txt
The output is a tabular three columns table: chromosome name, position, coverage. It is worth noting that samtools checks the CIGAR field and if your sample has a deletion will report zero coverage on the deletion. Most of the times this is not wanted (it is required for the variant calling process, but that’s another story).
Bedtools genome coverage
Bedtools has an even better option to calculate coverage, as it produces a standard BED file.
bedtools genomecov -bga FILE.bam > FILE.coverage.bed
The output is a standard BED file: a tabular file with these columns: chromosome, start, end, coverage. The difference with the previous output is that adjacent bases with the same coverage will be collapsed in the same interval.
Covtobed: from conda with love
We recently published covtobed, a simple tool to extract coverage information from sorted BAM files.
To install the tool you can use Miniconda:
conda install -y -c bioconda covtobed
And then, the usage to extract regions covered bedtween 0X (included) and 1X (excluded) is as simple as:
covtobed -m 0 -M 1 input.bam > not_covered.bed