Nonpareil 3: Fast Estimation of Metagenomic Coverage and Sequence Diversity

Nonpareil 3 can use parallelization across nodes and threads by using Message Passing Interface (MPI) and pthreads, respectively. We executed Nonpareil with the original alignment kernel by using both parallelization methods, independently and in combination, on a 2.3-Gbp test data set in order to evaluate the speedup. These optimizations resulted in up to 500 times faster computation of coverage. Compared to Amdahl’s law, Nonpareil 3 speedups corresponded to around 99.5% parallelization with MPI alone and around 99.8% parallelization with MPI and four threads, while the observed speedup for multithreads in one node was essentially linear (Fig. 1). In addition, the read redundancy can now be estimated by using perfect matches of one k-mer per query read corrected by sequencing error estimation, instead of the complete ungapped alignment. This implementation results in similar coverage and diversity estimates, as well as highly correlated projections of sequencing effort, but reduces the computational time by about 300 times (Table 1; see Fig. S1 in the supplemental material). We recommend using k = 24 (the default value in Nonpareil), the smaller value producing stable estimates (Fig. S2). The runtime of Nonpareil with alignment kernel is approximately proportional to N × Q × L², where N is the total number of sequencing reads, Q is the number of query sequencing reads searched against the complete metagenomic data set, and L is the length of a read. In contrast, the complexity of Nonpareil with k-mer kernel is simply N × L, i.e., directly proportional to the total data set size (Fig. S3). As the runtime for the k-mer kernel is independent of the read length, Nonpareil 3 is directly compatible with long-read data, although current technologies (e.g., PacBio or MinION) exhibit error rates that are higher than optimal for the current Nonpareil implementation (i.e., <5% expected sequencing error). Thus, additional optimizations will be necessary to allow reads with higher sequencing error to be analyzed by Nonpareil. In addition, since the runtime is not proportional to Q, it is practical to use a larger value for Q (default of 10,000 instead of 1,000 in the alignment kernel), increasing the precision of Nonpareil over that of the original version.

In order to reduce the impact of sequencing error on Nonpareil with k-mer kernel, we implemented a redundancy correction based on the sequencing read quality scores. Briefly, we estimated the fraction of selected k-mers that are expected to include at least one sequencing error from the quality scores and removed that fraction from the list of k-mers without matches, assuming that the introduced errors would most frequently result in unobserved variants (see Materials and Methods for details). To test this correction, low- and high-coverage simulated metagenome-like data sets of 101-bp reads (507,813 and 7,093,697 reads, respectively) were generated in silico from 30 bacterial and archaeal genomes by a previously described method (3) (Table S1; data sets are available at http://enve-omics.ce.gatech.edu/data/nonpareil). For each test, 10,000 reads were randomly selected from the simulated data sets and modified with substitution errors at a constant error probability per base. The 10,000 reads were then used in Nonpareil with k-mer kernel as query sequences with error correction enabled and disabled (Table S2; Fig. S4). In general, error rates affected the estimation of both coverage and required sequencing effort when error correction was disabled but not when it was enabled, indicating that the correction effectively removes the effect of sequencing errors when the error probability estimation is accurate.

The Nonpareil index of sequence diversity (N_d), described in Materials and Methods, is expressed in units of the natural logarithm of base pairs and summarizes the community diversity in sequence space, i.e., how redundant (or conversely unique) the sequences of a data set are among themselves. This metric depends on the joint distribution of genome size and abundance, as well as intragenome gene duplication. Therefore, given a small variation in genome size and a small impact of genomic duplications, e.g., for prokaryote-only communities, N_d can be used as a database-independent metric of alpha diversity. Since the shapes of the Nonpareil curves from replicates and subsamples closely resemble each other regardless of coverage (3), we propose N_d as a coverage-independent measurement of the diversity of the sampled community (Fig. S5).

We compared N_d and Shannon diversity index values in natural units (H′) from 16S rRNA gene OTU tables in different collections of metagenomic data sets, including: set I, metagenomes from different biomes including six mock data sets; set II, metagenomes and 16S rRNA gene amplicons from human microbiomes of different body sites; and set III, a collection of marine metagenomes from different sampling sites. First, we observed a high correlation (Pearson’s r = 0.804 [P < 10⁻¹⁵]; analysis of variance [ANOVA], 64.6% variance explained [P < 10⁻¹⁵]) between N_d and H′ values from diverse environments (set I), as well as a monotonic trend for interquartile ranges (IQRs) between different biomes (Fig. 2). This correlation was maintained in the subset of mock samples alone (r = 0.87 [P = 0.023]), and the residuals were only slightly affected by sequencing technology (ANOVA after including N_d, 3.7% variance [P = 0.008]). However, most outliers had low Turing-Good coverage estimates for the 16S rRNA gene OTU count profiles (labeled points in Fig. 2). Indeed, we observed a significant effect of the 16S rRNA gene Turing-Good coverage estimates on the residuals (ANOVA after N_d, excluding mock samples, 5.7% variance [P = 10⁻⁴]; ANOVA after N_d, including mock samples assuming complete coverage, 14.7% variance [P = 10⁻¹¹]) (Fig. 2, inset). We did not observe a significant effect of the Nonpareil-estimated coverage on the model residuals (ANOVA after N_d, 0.008% variance [P = 0.9]), indicating that a significant component of the difference between the two estimates of diversity was due to insufficient data to robustly estimate H′ but not N_d. We were able to differentiate between the effect of coverage on both estimates independently because the correlation between 16S rRNA gene-derived Turing-Good coverage estimates and metagenome-derived Nonpareil coverage estimates was partially lost by subsampling of most metagenomic data sets in the estimation of N_d but not in the estimation of H′ (r = 0.64). (Note that H′ was estimated on the basis of previously constructed OTU tables on the complete data sets available [11], not subsamples as in the case of N_d; see Materials and Methods for details.)

Next, we evaluated paired metagenomic and 16S rRNA gene amplicon samples from the Human Microbiome Project (HMP) (set II) and observed a similarly high correlation between median N_d and median H′ values per body site (r = 0.93), as well as between N_d and H′ values per data set (r = 0.84; Fig. S6).

Finally, we explored a collection of marine data sets (set III) and compared N_d and H′ diversity indices against available metadata by using ANOVA with alpha = 0.01. After controlling for the variance introduced by the source project (as the first independent variable in ANOVAs), we evaluated the effects of different variables on diversity variation (see Materials and Methods). For both measurements of diversity (N_d and H′), we observed significant effects of size fraction (H′ variance explained, 16%; N_d variance explained, 25% [P < 10⁻¹¹]) and latitude (H′ variance explained, 6.7%; N_d variance explained, 1.5% [P < 0.01]), similar to previously reported results for estimated richness as the measure of alpha diversity (12). Geographic location was also found to significantly affect the variance of N_d (12% variance, P = 10⁻⁸) and much more weakly so the variance of H′ (4.3% variance, P = 0.02). Moreover, N_d captured a pattern of increased diversity toward the winter (wintriness; explained variance, 2.3% [P = 0.0011]) not observed with H′ (explained variance, 0.14% [P = 0.46]).

The processing of a sample in Nonpareil is divided into two main steps. The first is redundancy estimation, where sequences are compared and the estimated redundancy is subsampled at different values of sequencing effort. The large number of pairwise alignments makes this step the most resource intensive. This task is now distributed across multiple nodes using MPI and processors using C++ pthreads. Further, Nonpareil 3 offers a k-mer-based method for redundancy estimation (see below) that accelerates this step by using short fragments of the sequencing reads with no errors allowed. Finally, redundancies are subsampled in multiple threads. In Nonpareil 1, subsamples were estimated linearly, resulting in sparse values toward the left side of the Nonpareil curve. Although this strategy is still available, the default in Nonpareil 3 is logarithmic subsampling; sample size is iteratively multiplied by a density factor (default 0.7) until only two reads remain.

The second step is estimation of the abundance-weighted average coverage at different sequencing efforts (Nonpareil curves), fitting to a sigmoidal model (projection), and graphical representation (e.g., Fig. S1). Note that this step relies on the assumption of independence of events between sequencing reads; therefore, Nonpareil should be applied on single reads, one of the sister reads in paired-end reads, or merged paired-end reads (3). This step has modest resource requirements and is implemented in the Nonpareil R package. In Nonpareil 3, we have streamlined this analysis by using headers in the redundancy output files and included an estimation of the sequence diversity derived from the fitted model (see below). All of the experiments in this report were executed with Nonpareil v3.3.

To accelerate read-to-read comparisons, Nonpareil 3 now implements a k-mer-based comparison kernel. Briefly, query k-mers are derived from a randomly selected subset of the metagenomic data set reads, with one target k-mer selected from the 5′ end of each read. Determination of the number of times each target k-mer (or its reverse complement) is found at any position of any sequence read in the complete data set can be performed in time proportional to the size of the metagenomic data set and is independent of the number of target k-mers. Note that this test is unable to detect if any of the last k − 1 bases in a sequencing read cover the target k-mer, so for each metagenomic read of length L, only L − k + 1 positions can be tested for matches, reducing the effective size of the metagenomic data set scanned with respect to the complete alignment kernel. Nonpareil curves (R package) account for the effect of this difference between sequencing effort and effective size, and its effect is corrected.

Because the target k-mers are derived from the metagenomic sequence reads, some of the k-mers will contain sequence errors, but if k is large enough, these k-mers will likely have zero coverage. Although it is not possible to determine which of the zero-coverage target k-mers contain errors, we utilized the sequencing quality (Q) scores from each target k-mer’s parent sequencing read to estimate E, the expected number of k-mers with errors. To correct for these errors, we reduced the target set by removing E zero-coverage target k-mers. The coverage counts for the remaining target k-mers, along with the reduced effective metagenomic data set size, are passed through the remaining original Nonpareil steps to estimate coverage, required sequencing effort, and sequence diversity.

To compare the results obtained with the traditional alignment kernel and the novel k-mer kernel, we executed both analyses in seven metagenomic data sets with different degrees of diversity (Table 1; Fig. S1). Metagenomic data sets were processed by using SolexaQA (21) with a maximum expected error of 1% and a minimum length of 50 bp, and adapter contamination was clipped by using Scythe (https://github.com/vsbuffalo/scythe). For paired-end samples, only the forward reads were used. Short Read Archive (SRA) identifiers are provided in Table 1 for all of the data sets except Iowa continuous cornfield soil. For the latter, seven lanes from one run of Illumina HiSeq were retrieved from the JGI Genome Portal (http://genome.jgi.doe.gov) on 21 July 2016 from project 402461.

Nonpareil results with k-mer kernel were obtained by using a 27-in. iMac with 8 gigabytes of random-access memory and an Intel Core i5 3.2-GHz processor using k = 24, 10,000 queries, and two threads (only for subsampling, the k-mer matching portion of Nonpareil is single threaded). Nonpareil alignment kernel results for Iowa soil were obtained by using the Michigan State University High-Performance Computing Center nodes with 20 cores and the following options: 20 threads, 1,000 queries, 50% overlap, and a redundancy to coverage transformation factor of 1.0. All other data sets were processed with a 27-in. iMac as described above and with the same Nonpareil alignment options at 95% identity and two threads. Iowa soil central processing unit (CPU) time for alignment kernel was estimated by using the linear regression of the other data sets.

Nonpareil curves are plots of abundance-weighted average coverage (Ĉ) per sequencing effort (LR) fitted to the cumulative probability function of the gamma distribution (3) with parameters α and β as follows: Ĉ = γ[α, β × log(LR + 1)]/Γ(α), where Γ is the gamma function and γ is the lower incomplete gamma function. Hence, we can use the mode of the corresponding gamma distribution to identify the value of log(LR + 1) corresponding to the inflection point of the curve, which we propose as a measurement of sequence diversity as follows: N_d = (α − 1)/β.

To evaluate the correlation between N_d and traditional measurements of alpha diversity derived from taxonomic affiliation, we compiled three collections of metagenomic data sets. Set I consisted of 90 samples from multiple environments, set II consisted of 54 human-associated microbiome metagenomic samples, and set III consisted of 292 marine samples. Set I was used to evaluate the general agreement between traditional taxonomy-based alpha diversity and N_d. It consisted of 84 metagenomic data sets from multiple environments divided into six distinct biomes plus six mock data sets (Table S3). Subsamples of processed nucleotide reads (between one and seven file chunks of 0.5 gigabyte zipped, as necessary to surpass 60% coverage or as many files as available) and OTU tables based on metagenome-derived 16S rRNA gene-containing reads were obtained from EBI Metagenomics (11). The Shannon index (H′) of each OTU table was estimated on the basis of Bayesian estimates of frequencies by using the Dirichlet multinomial pseudocount model with Laplace prior, as implemented in the R package entropy (22), with the exception of mock samples, for which maximum-likelihood entropy of input concentrations was used. Note that the 16S rRNA gene-containing reads derived from metagenomes do not necessarily overlap; hence, the EBI Metagenomics Pipeline uses closed-reference OTU picking (23), potentially biasing the results by database completeness. Therefore, we extended this analysis by using set II derived from the HMP (24). This collection consisted of samples from 13 body sites including 54 metagenomic data sets and 3,613 16S rRNA gene amplicon data sets. Fifty-three samples included both types of data sets. An OTU table including all of the processed samples was obtained from HMP Qiime Community Profiling (23, 24), from which H′ values per sample were estimated. Values were compared against N_d values derived from the metagenomic data sets by body site (medians per site) and by sample (intersection samples). Finally, we evaluated the resolution of N_d compared to H′ by using set III, a collection of marine samples derived from two global sampling projects, 228 samples from the Tara Oceans expedition between September 2009 and March 2012 (12) and 64 samples from the Global Ocean Sampling expedition between March 2009 and December 2010 funded by the Beyster Family Fund and the Life Technology Foundation (SRA BioProject accession no. PRJEB10418). Sample metadata, as well as preprocessed reads and 16S rRNA gene OTU tables, were obtained from EBI Metagenomics (11). ANOVA was performed to evaluate the effects of different metadata variables on both N_d and Bayesian analysis-corrected H′ of extracted 16S rRNA genes. The variables considered were size fraction (categorical), absolute value of latitude (i.e., degrees from the equator), latitude (degrees), geographic location (Mediterranean Sea, North Atlantic Ocean, North Pacific Ocean, Indian Ocean, Red Sea, South Atlantic Ocean, South Pacific Ocean, or Southern Sea), sampling date, and two decomposed seasonal components of sampling date. The decomposed seasonal components were estimated as the sine (vernality) and cosine (wintriness) of the date in radians (1 year = 2π) for samples in the Northern Hemisphere or the negative sine and negative cosine, respectively, for samples in the Southern Hemisphere.

Nonpareil 3 is available for online analyses at http://enve-omics.ce.gatech.edu/nonpareil. The Nonpareil code is freely distributed under the artistic license 2.0 and is available at https://github.com/lmrodriguezr/nonpareil.