Estimating and comparing microbial diversity in the presence of sequencing errors

Introduction

Advances in high-throughput DNA sequencing have opened a novel way to assess highly diverse microbial communities (Sogin et al., 2006; Roesch et al., 2007; Fierer et al., 2008; Turnbaugh & Gordon, 2009). However, the measurement and comparison of microbial diversity are challenging issues due to sampling limitations (Bohannan & Hughes, 2003; Schloss & Handelsman, 2006; Schloss & Handelsman, 2008; Øvreås & Curtis, 2011). These issues become more challenging when sequencing errors generate spurious low frequency counts especially singletons (Quince et al., 2009; Dickie, 2010; Kunin et al., 2010; Quince et al., 2011; Bunge et al., 2012a; Bunge et al., 2012b). In this paper, we use “species” to refer to taxa or operational taxonomic units (OTUs) under a pre-specified percentage of identity of sequences (Schloss & Handelsman, 2005; Schloss & Handelsman, 2008). We also use “individuals” to refer to sequences or any sampling unit.

In macro-ecology, Hill numbers have been increasingly used to quantify species diversity. An Ecology Forum led by Ellison (2010) (and papers that followed it) surprisingly achieved a consensus in the use of Hill numbers as the proper choice of diversity measure, despite intense debates existing in earlier literature regarding this issue. Hill numbers (or the effective number of species) are a mathematically unified family of diversity indices differing among themselves only by an exponent q that determines the measure’s sensitivity to species relative abundances. This family includes the three most important diversity measures: species richness (q = 0), Shannon diversity (q = 1, the exponential of Shannon entropy), and Simpson diversity (q = 2, the inverse of Simpson index). See below for its mathematical formula and interpretation. Hill numbers were first used in ecology by MacArthur (1965), developed by Hill (1973), and reintroduced to ecologists by Jost (2006) and Jost (2007). Hill numbers have been extended to incorporate evolutionary history and species traits; see Chao, Chiu & Jost (2014) for a recent review.

Various ecological measures have been applied to quantify the diversity of microbial communities (Hughes et al., 2001; Curtis, Sloan & Scannell, 2002). Hill et al. (2003) reviewed and discussed the suitability of a wide range of ecological diversity measures for use with highly diverse bacterial communities. Members of Hill numbers are also proposed as promising measures for quantifying microbial diversity. For example, Haegeman et al. (2008), Haegeman et al. (2013) and Haegeman et al. (2014) recommended the use of Shannon diversity and Simpson diversity to measure and compare microbial diversity; Doll et al. (2013) suggested using a continuous diversity profile, a plot of Hill numbers as a continuous function of q ≥ 0. In this paper, we adopt the general framework of Hill numbers and use continuous profiles to quantify microbial diversity. The diversity profile for q ≥ 0 conveys all information contained in a species relative abundance distribution if community parameters (species richness and relative abundances) are known. In practice, however, community parameters are unknown and thus the true diversity must be estimated from sampling data, meaning that statistical methods are required.

In this paper, we propose two statistical approaches (non-asymptotic and asymptotic) to make fair comparisons of microbial diversity across multiple communities. A non-asymptotic approach refers to the comparison of estimated diversities of standardized samples with a common finite sample size or sample completeness (as measured by sample coverage; see below). This approach aims to compare diversity estimates for equally-large or equally-complete samples; it is based on the seamless rarefaction and extrapolation sampling curves of Hill numbers, specifically for q = 0, 1 and 2. Traditional sample-size-based rarefaction for species richness has been widely applied in ecology as a standardization method and also suggested by Dickie (2010) for molecular surveys. For species richness, Colwell et al. (2012) proposed an integrated rarefaction and extrapolation sampling curve for standardizing sample size; Chao & Jost (2012) proposed the corresponding curve for standardizing sample completeness. Hill numbers calculated from a sample, like species richness, are an increasing function of sampling effort and thus tend to increase with sample completeness. Chao et al. (2014) generalized previous papers (Chao & Jost, 2012; Colwell et al., 2012) on species richness to the family of Hill numbers and developed two types of standardization methods (sample-size- and sample-coverage-based). The sample-size- and sample-coverage-based integration of rarefaction and extrapolation together represent a unified non-asymptotic and non-parametric framework for estimating diversity and for making statistical inferences based on these estimates. The rarefaction and extrapolation curves for measures of small value of q (say, 0 ≤ q < 2) heavily depend on the low frequency counts, especially singletons (Chao et al., 2014).

An asymptotic approach refers to the comparison of the estimated asymptotic diversity profiles. That is, this approach compares the estimated profiles for samples with size tending to be sufficiently large or samples with sample completeness tending to unity. It is based on statistical estimation of the true Hill number of any order q ≥ 0. This profile is typically generated by substituting species sample proportions into the diversity formula. However, this empirical approach generally underestimates the true profile, because samples usually miss some of the community’s species due to under-sampling. Finding an analytic reduced-bias continuous diversity profile has been a long-standing challenge. Chao & Jost (2015) recently proposed a resolution to obtain a diversity profile estimator, which infers the asymptotic or true diversities. By applying their diversity profile estimator, the negative bias associated with the empirical diversity curve due to undetected species can be greatly reduced. The authors also used real data sets to demonstrate that the empirical and their estimated diversity profiles may give qualitatively different answers when comparing biodiversity surveys. Chao & Jost’s (2015) diversity profile estimator for low values of q (0 ≤ q < 2) is strongly affected by the low frequency counts. This is mainly because the observed rare species that produce low frequencies carry nearly all the information about the undetected species and play an important role in almost all statistical inferences in diversity estimation.

However, unlike macro-community ecological data, the low frequency counts, especially singletons from high-throughput DNA sequencing, are subject to various types of sequencing errors at different stages of processing (Quince et al., 2009; Huse et al., 2010; Quince et al., 2011). Some sequences may be misclassified as new taxa, and, accordingly, are also misclassified as singletons. Consequently, the observed singletons are greatly inflated and can comprise more than 60% of taxa in a sample (Buee et al., 2009). Since singletons play crucial roles in both asymptotic and non-asymptotic analyses as described above, our suggested approaches will be seriously affected if the inflated singleton count is not adjusted. A wide range of methods have been developed to reduce or correct sequencing errors (Buee et al., 2009; Quince et al., 2011) at the bioinformatics-processing stage. Without knowledge of the sources of measurement errors, statistical sampling-based methods have also recently been proposed to correct the number of spurious singletons and estimate diversity. Bunge et al. (2012a) and Bunge, Willis & Walsh (2014) proposed a parametric mixture model and a method using “left-censored” data; Willis & Bunge (2015) proposed an approach using the ratio of two successive frequency counts. These statistical approaches generally require different parametric assumptions about the sampling models or about the functional form of the ratio of frequency counts. Some of these parametric assumptions may not be reliably tested, and the different parametric assumptions may disallow comparisons across communities.

In this paper, we propose a novel nonparametric approach to estimate the true number of singletons in the presence of sequencing errors. Here we derive a relationship between the expected frequency of singletons and the expected frequencies of doubletons, tripletons and quadrupletons, based on a modified Good–Turing frequency formula originally developed by the founder of modern computer science Alan Turing, and his colleague Good (1953) and Good (2000). Our estimator of singleton count is thus expressed in terms of the observed frequency counts of doubletons, tripletons and quadrupletons, provided these three frequency counts are reliable. Simulation results are reported to demonstrate an important finding about our proposed singleton count estimator. That is, when there are no sequencing errors and sample sizes are reasonably large, our estimator differs from the true singleton count only to a limited extent; when there are sequencing errors, our estimator is substantially lower than the observed singleton count. Therefore, the discrepancy between the estimated and the observed singleton counts can also be used to assess whether or not sequencing errors were present in the observed data.

Throughout the paper, “adjusted data/estimators” refer to those with the observed singleton count being replaced by the estimated count (the observed singleton count is discarded), whereas “original or observed data” refer to the observed data with spurious singletons possibly present. To quantify and compare microbial diversity, here we propose applying both non-asymptotic and asymptotic analyses to the adjusted data whenever the singleton count is uncertain in measurement. That is, for adjusted data, we present seamless sample-size- and coverage-based rarefaction and extrapolation sampling curves of Hill numbers (focusing on measures of q = 0, 1, and 2) and a continuous diversity profile estimator. Simulation results based on various taxa abundance distributions are reported to examine the performance of our method and to compare our results with those obtained from a previous ratio-based method (Bunge, Willis & Walsh, 2014; Willis & Bunge, 2015). Sequencing data from viral metagenomes (Allen et al., 2011; Allen et al., 2013) are used for illustration and comparison. The generalization of our methods to phylogenetic diversity is discussed.

Methods

Simulation Results

Since both non-asymptotic and asymptotic analyses depend on the quality of the estimated singleton count, it is essential to investigate the performance of the proposed estimator in Eq. (5). We conducted a simulation by generating data from six species abundance distributions with various degrees of heterogeneity in species relative abundances (details are provided in Supplemental Text S2). In each model, we fixed the number of species at S = 2,000 to mimic microbial communities. Then for each given model, we considered a range of sample sizes (n = 2,000 to 10,000 in an increment of 2,000). The degree of heterogeneity in species relative abundances is quantified by the CV (which is the ratio of the standard deviation over the mean) of species relative abundances. When all species relative abundances are equal, CV = 0. A larger value of CV indicates a higher degree of heterogeneity among species abundances.

For each combination of abundance model and sample size, we generated two types of data: (i) data without sequencing errors, and (ii) spurious data with a sequencing error rate of 10%, i.e., there was a 10% chance that a sampled individual was misclassified to a new species and thus became a spurious singleton. In Fig. 1, we show the plots of the average values (over 1,000 simulation trials) of four singleton counts as a function of sample size that was used in data generation. The four singleton counts include the true singleton count generated from the data without sequencing error, the spurious singleton count generated from the data with sequencing error, the adjusted singleton count based on Eq. (5), and the count obtained from the ratio-based method of Bunge, Willis & Walsh (2014) and Willis & Bunge (2015) through the R package “breakaway,” available from CRAN (Comprehensive R Archive Network). In Fig. 2, the corresponding root mean squared errors (RMSEs) for the ratio-based and the proposed methods are shown. The patterns revealed by these plots are summarized below.

Figure 1 reveals that the number of singletons for the data without sequencing error (dotted curve in each panel) generally declines with sample size when sample size becomes sufficiently large, whereas the number of singletons for data with sequencing error (dash-dotted curve in each panel) always increases with sample size, revealing a drastically different pattern; see Dickie (2010) for a similar finding. This pattern can be used to detect whether sequencing error exists in the original data when an empirical accumulation curve for the singleton count can be recorded in the data-collecting procedures.

Simulation results also show that when heterogeneity is low as reflected by low CV values (Model 1 to Model 4) the ratio-based method (dashed curve) and our proposed method (solid curve) yield similar singleton counts that are close to the true data (i.e., data without sequencing error, dotted curve). The RMSEs of the two methods are thus generally comparable (Fig. 2). However, in the highly heterogeneous cases as reflected by relatively high CV values (Model 5 and Model 6), the ratio-method produces much higher singleton counts compared to the true data and thus much larger root mean squared errors than the proposed method, as shown in Fig. 2. In these high-CV cases, our estimator of singleton count still closely matches the true number of singletons, although it exhibits negative bias when sample size is relatively small especially when species abundances are highly heterogeneous.

These simulation results thus imply (i) when there are no sequencing errors (so that the dotted curves represent the singleton counts for data), our estimator differs only to a limited extent from the true data, yielding almost the same diversity inference; (ii) when there are sequencing errors (so that the dash-dotted curves represent the spurious singleton counts for data), our estimator can greatly reduce the raw singleton count and make proper corrections. Therefore, the discrepancy between our proposed estimator of singleton count and the singleton count from the observed data can be used to assess whether sequencing errors were present in data processing. Moreover, this implies that whenever the singletons are uncertain or in doubt, it is worth applying our proposed estimator of singleton count. More simulation results on the effect of spurious singletons on the estimation of asymptotic diversities are provided in Supplemental Text S2; see Discussion.

Application Results

A number of data sets on frequency counts of contig (contiguous groups of sequences) spectra of viral phage metagenomes from similar or different environments were analyzed in Allen et al. (2013). We select two samples with different environments to illustrate the use of our methods: one sample includes the pooled contig spectra from seven non-medicated swine feces, and the other sample includes the pooled contig spectra from four reclaimed fresh water samples. For simplicity, these two samples/viromes are respectively referred to as “swine feces” sample/virome and “reclaimed water” sample/virome in the following analysis. The frequency counts for the two samples originally provided in the additional file of Allen et al. (2013) are reproduced in Table 1. The empirical diversities and asymptotic analyses are shown in Table 2.

In the swine feces original data, there were 8,833 taxa among 9,988 individuals (sequences); the number of singletons was f₁ = 8,025, and the number of doubletons was f₂ = 605. In the reclaimed water data, there were 8,739 taxa among 9,973 individuals, and the first two frequency counts are f₁ = 7,986 and f₂ = 518. In these two original samples, most of the frequencies are concentrated on singletons. Consequently, based on the original data, the Chao1 lower bounds, 62,057 and 70,299 respectively for swine feces and reclaimed water viromes, are greatly inflated due to the presence of spurious singletons. Using Eq. (5), we obtain an estimated singleton count of 2,831 for the swine feces sample, and 2,105 for the reclaimed water sample (Table 1). For each sample, the estimated singleton count is substantially less than the observed singleton count, indicating that sequencing errors were present. The empirical and estimated diversities for the original and adjusted data are shown in Table 2. We also compare in Table 2 our estimates with those based on a ratio-based method (Bunge, Willis & Walsh, 2014; Willis & Bunge, 2015), and with those proposed in Allen et al. (2013) based on the CatchAll software.

From Table 2, as expected, the estimated diversity (species richness, Shannon diversity and Simpson diversity) based on the adjusted data for each sample is much lower than that based on the original data. For species richness, the CatchAll software yields excessively low estimates, even lower than the observed richness of the adjusted data. The ratio-based method, however, yields extremely large estimates for the number of species. In our simulations on species richness as described in Supplemental Text S2, we show that the ratio-based method might severely overestimate the true species richness when the heterogeneity among species abundance is relatively high. The empirical CV values for the swine feces and reclaimed water samples for adjusted date are respectively 0.62 and 0.79. As there are many undetected rare species, the true CV should be much higher than the empirical CV, leading to extremely large species richness estimates for the ratio-based method. All the following analyses are based on our adjusted data, unless otherwise stated.

Before we present the non-asymptotic analyses, we plot in Fig. 3 the sample completeness curve as a function of sample size. The sample completeness of the adjusted swine feces sample is 41%, which is lower than that for the adjusted reclaimed water sample, 48.6%. When the sample size is extrapolated to a size of 10,000 (approximately double the adjusted sample size for swine feces), the coverage of the swine feces sample is increased from 41.0% to 62.9%, whereas the coverage of the reclaimed water sample is increased from 48.6% to 74.7%. For any standardized sample size, Fig. 3 shows that the sample completeness of the swine feces sample is lower than that for the reclaimed water sample of the same size.

For non-asymptotic analysis, we present in Fig. 4 the sample-size- and coverage-based rarefaction and extrapolation curves along with 95% confidence intervals for three measures: q = 0, 1 and 2. The sample-size-based sampling curve is extrapolated up to a maximum size of 10,000, whereas the coverage-based sampling curve is extended up to the coverage of the size 10,000, i.e., the maximum coverage is up to 62.9% for the swine feces sample and 74.7% for the reclaimed water sample.

All plots in Fig. 4 exhibit a consistent pattern, with the diversity curve for the swine feces samples lying above the curve of the reclaimed water sample. In all plots, the 95% confidence intervals for the two samples in any rarefaction/extrapolation curve are disjoint, implying a significant difference. As stated earlier, the extrapolation for Shannon and Simpson diversity, unlike that of species richness, can often be reliably extended to infinity size or complete coverage to reach the asymptotic diversity estimate. Therefore, for Shannon diversity (common taxa richness) and Simpson diversity (dominant taxa richness), the data indicate that the swine feces virome is significantly more diverse than the reclaimed water virome. This is valid not only for the standardized sample size and sample coverage values plotted in Fig. 4, but also for entire viromes. (This is also supported by the asymptotic analysis below). For species richness, the data support this conclusion up to a standardized 62.9% fraction of each virome (Fig. 4B). Beyond that, the data do not provide sufficient information for comparison. This is because the asymptotic species richness estimator is only a lower bound (as opposed to point estimates for the other two asymptotic diversities).

For the asymptotic analysis, we plot the empirical and estimated asymptotic diversity profiles along with 95% confidence intervals in Fig. 5 when q is between 0 and 3. (The empirical and estimated asymptotes of diversities for the special cases of q = 0, 1 and 2 are shown in Table 2, and the asymptotic diversity estimates are also shown next to an arrow at the right-hand end of each rarefaction/extrapolation plot in Fig. 5). The empirical diversities (Table 2 and Fig. 5) imply that the two viromes have limited difference in each of the three measures. In contrast, the plots in Fig. 5 reveal that for the asymptotic Shannon diversity, the swine feces virome is substantially more diverse than the reclaimed water virome. A similar conclusion is also valid for the Simpson diversity, confirming our earlier statement in the preceding paragraph.

Table 2 and Fig. 5 show that the adjusted Chao1 estimator in Eq. (6a) gives an estimate of 10,261 taxa for swine feces and 7,134 taxa for reclaimed water virome. Each is five times that obtained from CatchAll (Allen et al., 2013). Since the Chao1 estimate represents only minimum richness, it cannot be used to rank the taxa richness of the two entire viromes. Nevertheless, taxa richness can be compared through the coverage-based non-asymptotic approach, as discussed earlier; see Discussion. By contrast, for diversity of order q ≥ 1, we can compare not only the estimated diversities for standardized sample size/completeness but also the estimated asymptotic diversities across communities. In Supplemental Table S1, we also give all the estimated asymptotes of diversities for other data sets provided in Allen et al. (2013).

Conclusion and Discussion

Whenever the singletons are uncertain or in doubt in sequencing data, it is worth applying our proposed estimator to estimate the singleton count; see Eq. (5). The discrepancy between our estimated singleton count and the observed count can be used to infer whether sequencing errors were present in data processing. Using the estimated number of singleton count and the original non-singleton frequency counts, we can quantify and compare microbial diversity for data sets with different sequencing error rates through non-asymptotic analysis (based on the plots of the sample-size- and coverage-based rarefaction and extrapolation sampling curves) and asymptotic analysis (based on the plot of a continuous asymptotic diversity profile estimator). Illustrative plots for sequencing data from viral metagenomes are shown in Fig. 4 (the non-asymptotic analysis) and Fig. 5 (the asymptotic analysis). Although we have focused on microbial data with spurious singleton counts, both our asymptotic and non-asymptotic approaches are also recommended for analyzing data with reliable singleton counts.

In highly diverse microbial communities, unless strong assumptions or parametric models are made, sampling data often do not provide sufficient information to accurately infer the number of undetected taxa in the sample. Thus, it is statistically infeasible to provide reliable estimates of taxa richness for the entire community. Our estimated species richness (q = 0 measure in our asymptotic analysis) theoretically is a lower bound. This implies that fair comparison of asymptotic species richness among multiple communities is not statistically feasible. In this case, fair comparison of taxa richness across multiple assemblages can be made by standardizing sample completeness (i.e., comparing taxa richness for a standardized fraction of population) based on coverage-based rarefaction and extrapolation sampling curves, as illustrated in the real data analysis. By contrast, when the diversity order q is away from 0 (say, q ≥ 1), rare species have less impact on these diversities, and we generally can infer these diversities up to asymptotes and compare them across communities; see our illustrative example for interpretations. We recommend the use of an estimated diversity profile such as Fig. 5 for asymptotic analysis. If only one or two measures are desired in the inferences of highly diverse microbial diversity, then a perspective from Shannon diversity and Simpson diversity, instead of taxa richness, is more promising and more practical because we can accurately estimate Shannon and Simpson diversity not only for standardized samples but also their asymptotes. Besides, as shown in our simulation results (Figure S1 in Supplemental Text S2), the taxa richness estimator is seriously inflated or affected by spurious singleton counts, whereas the effect on Shannon diversity and Simpson diversity is less serious.

Our proposed estimator of singleton count is in terms of f₂, f₃ and f₄, provided these counts are reliable. A slight generalization of our method can be applied to estimate any frequency count. For example, supposing that singletons and doubletons are both uncertain, we can similarly derive an estimator of doubleton count based on f₃, f₄ and f₅ following exactly the same approach proposed in this paper. Subsequently, Eq. (5) then gives an estimate of singleton count based on the estimated doubleton count, f₃ and f₄. Consequently, our proposed non-asymptotic and asymptotic analyses can be similarly applied to data with the first two frequency counts being replaced by the estimated values. However, the sampling variance of the estimated diversity would be unavoidably increased.

In our approach, the original singleton count is discarded and replaced by our estimated count. In 16S rRNA sequencing or metagenomic sequencing, it is often standard practice to compare sequencing reads against a reference database, such as Greengenes, e.g., see Turnbaugh et al. (2009) or the software MOTHUR (Schloss et al., 2009). The Greengenes alignment tool helps adjust the original singleton count and alleviate the problem of sequencing error. Also, when there are multiple samples, singletons in a given sample have different probabilities of being spurious depending on their total number of reads across samples. Because many local singletons are not global singletons, the cross-sample information may also help adjust the original singleton count. Further investigation examining how to extend our framework to incorporate related covariates (such as cross-sample and database information) is merited.

Finally, we briefly discuss the phylogenetic diversity (PD) because of its broad interest and applications (Martin, 2002; Lozupone & Knight, 2005) in microbial studies. In this paper, all taxa are treated as if they were equally distinct and thus differences among sequences are not considered. Faith’s (1992) PD is the most widely used PD metric to take into account phylogenetic differences among taxa. Faith’s (1992) PD is defined as the total sum of branch lengths of a phylogenetic tree connecting all focal species. Based on sampling data, Chao et al. (2015) recently proposed a non-parametric estimator of the true PD (PD of the entire community, i.e., the observed PD in the sample plus the undetected PD). In the presence of sequencing errors, the inflated singleton count will also affect the estimation of the true PD. Since error-induced singletons will mostly likely fall in a closely related taxon, the effect may not be as pronounced as that in species richness estimation. More investigation is needed to tackle sequencing error and to adjust the Chao et al. (2015) PD estimator. Since Faith’s (1992) PD does not incorporate taxa abundances, Chao, Chiu & Jost (2010) developed a class of abundance-sensitive PD measures which generalize Faith’s (1992) PD to incorporate taxa abundances, and also extend Hill numbers to take into account phylogenetic relationships among taxa. How to extend the proposed analyses presented in this paper (the asymptotic and non-asymptotic analyses) to the class of abundance-sensitive PD is a worthwhile topic of future research.

Supplemental Information