Estimating DNA coverage and abundance in metagenomes using a gamma approximation

Abstract

Motivation: Shotgun sequencing generates large numbers of short DNA reads from either an isolated organism or, in the case of metagenomics projects, from the aggregate genome of a microbial community. These reads are then assembled based on overlapping sequences into larger, contiguous sequences (contigs). The feasibility of assembly and the coverage achieved (reads per nucleotide or distinct sequence of nucleotides) depend on several factors: the number of reads sequenced, the read length and the relative abundances of their source genomes in the microbial community. A low coverage suggests that most of the genomic DNA in the sample has not been sequenced, but it is often difficult to estimate either the extent of the uncaptured diversity or the amount of additional sequencing that would be most efficacious. In this work, we regard a metagenome as a population of DNA fragments (bins), each of which may be covered by one or more reads. We employ a gamma distribution to model this bin population due to its flexibility and ease of use. When a gamma approximation can be found that adequately fits the data, we may estimate the number of bins that were not sequenced and that could potentially be revealed by additional sequencing. We evaluated the performance of this model using simulated metagenomes and demonstrate its applicability on three recent metagenomic datasets.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

The effect of community complexity on the expected number of reads per bin (λ_i) for each bin i. In a simple community with a single genome, λ is approximately equal for all bins. Here, the observed bin spectrum (number of reads per bin) follows a Poisson distribution. However, in metagenomic samples from complex communities, bins will be drawn from different genomes that are present in varying abundances. Therefore, the value of λ is not the same for all bins. If λ follows a gamma distribution, then the bin spectrum will follow a negative binomial distribution and can be modeled.

Fig. 2.

(a) Blue curve: estimated log bin spectrum for the Lake Washington formate dataset. Red stars: the log number of observed bins. Note that the observed value at zero reads per contig is zero. The χ²-score for this fit is 3.7; we cannot reject the assumption that the bin abundance is gamma-like. (b) Blue curve: observed and estimated log bin abundance distribution for the termite hindgut dataset. The χ²-value is 1.0.