CoverM: read alignment statistics for metagenomics

1Microbial communities significantly impact the health of our planet and ourselves personally. Through the advent of genome-centric metagenomics, it has become possible to recover complete or near-complete genomes from these populations through deep sequencing of community DNA. Piecing together sequencing reads into metagenome assembled genomes (MAGs) has become increasingly common and high-throughput, with recovery of tens to hundreds of thousands of genomes from a diverse array of environments.

2A key step in MAG recovery is binning, where assembled contigs are partitioned into groups representing putative genomes. Commonly, contig abundances across samples are calculated based on alignments of reads to contigs. Contigs with similar abundances across many samples are likely to be from the same genome, enabling improved, sequence-agnostic genome recovery (differential coverage binning). However, the specific algorithms metagenomic binning tools use to derive coverage statistics from read alignments [BAM files] vary substantially. For instance, some use the number of aligned reads divided by contig length. Others calculate a coverage based on the number of reads aligned to each reference position, taking the mean or trimmed mean across all bases in a contig. Still others find approximate coverage using a k-mer counting approach eschewing the alignment step altogether for faster calculation.

3The backend data structure used for coverage calculation can also vary. Notably, “Mosdepth arrays” were shown to be two times faster than naive methods for exact calculation of mean coverage. Mosdepth arrays record the change in the number of mapped segments at each position, such that coverage at a given position can be calculated by summing all preceding values.

4Recovered MAGs can be used to estimate the community composition by aligning reads to their contigs and calculating coverage statistics across the genome. For example, calculating the mean coverage for a genome is achieved by summing the number of aligning bases across all contigs in a genome and dividing by their total length. Researchers may also opt to apply filters to prevent genomes being deemed present (i.e. having non-zero coverage) on the basis of a small number of positions having aligned reads.

5Despite the differences between per-contig and per-genome coverage statistics, they are conceptually similar enough to motivate a single tool that calculates both. Here, we present CoverM, which uses Mosdepth arrays to efficiently calculate a number of per-contig and per-genome coverage statistics.

6CoverM is a software tool for calculating coverage statistics of contigs or genomes with reads derived from metagenomes or microbial isolates (Fig. 1). It is written in the Rust programming language, relying on independent tools for genome dereplication and read alignment. It makes substantial use of rust-htslib, which provides Rust bindings to HTSlib for parsing intermediate alignment (BAM) files.

7CoverM can be compiled locally using the standard cargo method in the Rust ecosystem, by downloading the source code from the repository directly or by downloading a version statically compiled with musl libc from https://github.com/wwood/coverm/releases. Apart from the standard command line version, Python and Julia interfaces to CoverM are also available.

8The inputs to CoverM are metagenomic reads (or precalculated BAM files) and reference sequences: either contigs or genomes. Since aligning against multiple near-identical reference genomes can distort read alignments, CoverM uses Galah to dereplicate input genomes, interpreting ANI calculated by skani, FastANI, or Dashing.

9Next, the reads are aligned to the reference sequences using Strobealign, Minimap2, or BWA. Spurious alignments can be filtered on the basis of identity, alignment length, or percent of the read aligned. Paired reads can be filtered independently or together with their mate. The filtered alignments are then processed using Mosdepth arrays to calculate a range of coverage statistics (Table 1). The coverage of a genome is aggregated by considering all bases in the genome together. As an additional filtering step, to account for off-target alignment in conserved regions, by default CoverM requires that at least 10% of a genome’s length is covered by at least one read before assigning non-zero coverage.

10The default per-genome coverage statistic is relative abundance of each genome [Equation (1)]. The relative abundance of a genome is an estimate of the proportion of the cells in a given sample belonging to the species with that genome. Assuming no bias in DNA extraction, sequencing or read alignment, the ratio of genome copies in a sample is equal to the ratio of their mean coverages. One complication is that in many metagenomes the reference database is incomplete, so many reads remain unaligned. To account for this, the fraction of reads that are unaligned is used to scale the relative abundances of each genome.

11In Equation (1), is the estimated relative abundance of the lineage as a percentage, is the calculated mean coverage of that lineage, is the sum of the mean coverage of all genomes, is the number of reads aligning to any genome, and is the total number of reads in the sequencing library. Implicitly, this formula assumes that (i) the average genome size of community members is the same between genomes in the reference set and those that are missing, and (ii) species with reference genomes are roughly complete and uncontaminated.

12CoverM provides an efficient calculation method for several coverage metrics applicable to contig-based (metagenomic binning) and genome-based (community profiling) tasks. The implementation is flexible, with >50 arguments and options for controlling input/output formats, alignment, filtering, dereplication, and coverage metric.

13Conflict of interest: None declared.