Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible — TEI-tagged view

Paul J. McMurdiepersName, Susan HolmespersName (2014date). “Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible.” PLoS Computational Biologytitle·j.

Object of study — normalization and statistical inference for overdispersed microbiome sequencing-count data

Research question — whether rarefyingterm is statistically defensible for microbiome clustering and differential-abundance analysis, and whether mixture-model methods perform better

Contribution — a formal critique plus reproducible empirical-data simulationsterm·mth showing the costs of rarefyingterm and supporting negative-binomial alternatives

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Waste Not, Want Not: Why Rarefyingterm Microbiome Data Is Inadmissible Abstract Current practice in the normalization of microbiome count dataterm is inefficient in the statistical sense. For apparently historical reasons, the common approach is either to use simple proportionsterm (which does not address heteroscedasticityterm) or to use rarefyingterm of counts, even though both of these approaches are inappropriate for detection of differentially abundant species. Well-established statistical theory is available that simultaneously accounts for library sizeterm differences and biological variability using an appropriate mixture model. Moreover, specific implementations for DNA sequencing read count data (based on a Negative Binomialterm model for instance) are already available in RNA-Seq focused R packages such as edgeRterm and DESeqterm. Here we summarize the supporting statistical theory and use simulationsterm·mth and empirical data to demonstrate substantial improvements provided by a relevant mixture model framework over simple proportionsterm or rarefyingterm. We show how both proportionsterm and rarefied counts result in a high rate of false positivesterm in tests for species that are differentially abundant across sample classes. Regarding microbiome sample-wise clustering, we also show that the rarefyingterm procedure often discards samples that can be accurately clustered by alternative methods. We further compare different Negative Binomialterm methods with a recently-described zero-inflated Gaussianterm mixture, implemented in a package called metagenomeSeqterm . We find that metagenomeSeqterm performs well when there is an adequate number of biological replicatesterm, but it nevertheless tends toward a higher false positive rateterm. Based on these results and well-established statistical theory, we advocate that investigators avoid rarefyingterm altogether. We have provided microbiome-specific extensions to these tools in the R package, phyloseqterm.

This is a PLOS Computational Biology Methods article.

Introduction Modern, massively parallel DNA sequencing technologies have changed the scope and technique of investigations across many fields of biology [1] , [2] . In gene expression studies the standard measurement technique has shifted away from microarray hybridization to direct sequencing of cDNA, a technique often referred to as RNA-Seq [3] . Analogously, culture independent [4] microbiome research has migrated away from detection of species through microarray hybridization of small subunit rRNA gene PCR amplicons [5] to direct sequencing of highly-variable regions of these amplicons [6] , or even direct shotgun sequencing of microbiome metagenomic DNA [7] . Even though the statistical methods available for analyzing microarray data have matured to a high level of sophistication [8] , these methods are not directly applicable because DNA sequencing data consists of discrete counts of sequence reads rather than continuous values derived from the fluorescence intensity of hybridized probes. In recent generation DNA sequencing the total reads per sample ( library sizeterm ; sometimes referred to as depths of coverageterm ) can vary by orders of magnitude within a single sequencing run. Comparison across samples with different library sizesterm requires more than a simple linear or logarithmic scaling adjustment because it also implies different levels of uncertainty, as measured by the sampling varianceterm of the proportionterm estimate for each feature (a feature is a gene in the RNA-Seq context, and is a species or O perational T axonomic U nit, OTUterm, in the context of microbiome sequencing). In this article we are primarily concerned with optimal methods for addressing differences in library sizesterm from microbiome sequencing data.

Variation in the read counts of features between technical replicatesterm have been adequately modeled by Poisson random variables [9] . However, we are usually interested in understanding the variation of features among biological replicatesterm in order to make inferences that are relevant to the corresponding population; in which case a mixture model is necessary to account for the added uncertainty [10] . Taking a hierarchical model approach with the Gamma-Poissonterm has provided a satisfactory fit to RNA-Seq data [11] , as well as a valid regression framework that leverages the power of generalized linear models [12] . A Gamma mixture of Poisson variables gives the negative binomialterm (NB) distribution [10] , [11] and several RNA-Seq analysis packages now model the counts, K , for gene i , in sample j according to: (1) where s j is a linear scaling factor for sample j that accounts for its library sizeterm, μ i is the mean proportionterm·var for gene i , and φ i is the dispersion parameterterm for gene i . The variance is , with the NB distribution becoming Poisson when φ = 0. Recognizing that φ >0 and estimating its value is necessary in gene-level tests in order to control the rate of false positiveterm genes. Many false positiveterm genes appear significantly differentially expressed between experimental conditions under the assumption of a Poisson distributionterm, but are nevertheless not-significant in tests that account for the larger variance that results from non-zero dispersion.

The uncertainty in estimating φ i for every gene when there is a small number of samples — or a small number of biological replicatesterm — can be mitigated by sharing information across the thousands of genes in an experiment, leveraging a systematic trend in the mean-dispersion relationshipterm [11] . This approach substantially increases the power to detect differences in proportionsterm (differential expression) while still adequately controlling for false positivesterm [13] . Many R packages implementing this model of RNA-Seq data are now available, differing mainly in their approach to modeling dispersion across genes [14] . Although DNA sequencing-based microbiome investigations use the same sequencing machines and represent the processed sequence data in the same manner — a feature-by-sample contingency table where the features are OTUsterm instead of genes — to our knowledge the modeling and normalization methods currently used in RNA-Seq analysis have not been transferred to microbiome research [15] – [17] .

Instead, microbiome analysis workflows often begin with an ad hoc library sizeterm normalization by random subsampling without replacementterm, or so-called rarefyingterm [17] – [19] . There is confusion in the literature regarding terminology, and sometimes this normalization approach is conflated with a non-parametric resampling technique — called rarefactionterm·mth [20] , or individual-based taxon re-sampling curves [21] — that can be justified for coverage analysis or species richness estimation in some settings [21] , though in other settings it can perform worse than parametric methods [22] . Here we emphasize the distinction between taxon re-sampling curves and normalization by strictly adhering to the terms rarefyingterm or rarefied counts when referring to the normalization procedure, respecting the original definition for rarefactionterm·mth . Rarefyingterm is most often defined by the following steps [18] .

Select a minimum library sizeterm, . This has also been called the rarefactionterm·mth level [17] , though we will not use the term here.

Discard libraries (microbiome samples) that have fewer reads than .

Subsample the remaining libraries without replacement such that they all have size .

Often is chosen to be equal to the size of the smallest library that is not considered defective , and the process of identifying defective samples comes with a risk of subjectivity and bias. In many cases researchers have also failed to repeat the random subsampling step (3) or record the pseudorandom number generation seed/process — both of which are essential for reproducibility. To our knowledge, rarefyingterm was first recommended for microbiome counts in order to moderate the sensitivityterm of the UniFrac distance [23] to library sizeterm, especially differences in the presence of rare OTUsterm [24] . In these and similar studies the principal objective is an exploratory/descriptive comparison of microbiome samples, often from different environmental/biological sources; a research task that is becoming increasingly accessible with declining sequencing costs and the ability to sequence many samples in parallel using barcoded primers [25] , [26] . Rarefyingterm is now an exceedingly common precursor to microbiome multivariate workflows that seek to relate sample covariates to sample-wise distanceterm matrices [19] , [27] , [28] ; for example, integrated as a recommended option in QIIME's [29] beta_diversity_through_plots.py workflow, in Sub.sample in the mothur software library [30] , in daisychopper.pl [31] , and is even supported in phyloseqterm's rarefy_even_depth function [32] (though not recommended in its documentation). The perception in the microbiome literature of “rarefyingterm to even sampling depth” as a standard normalization procedure appears to explain why rarefied counts are also used in studies that attempt to detect differential abundanceterm of OTUsterm between predefined classes of samples [33] – [37] , in addition to studies that use proportionsterm directly [38] . It should be noted that we have adopted the recently coined term differential abundanceterm [39] , [40] as a direct analogy to differential expression from RNA-Seq. Like differentially expressed genes, a species/OTUterm is considered differentially abundant if its mean proportionterm·var is significantly different between two or more sample classes in the experimental design.

Statistical motivation Despite its current popularity in microbiome analyses rarefyingterm biological count data is statistically inadmissible because it requires the omission of available valid data. This holds even if repeated rarefyingterm trials are compared for stability as previously suggested [17] . In this article we demonstrate the applicability of a variance stabilizationterm technique based on a mixture model of microbiome count dataterm. This approach simultaneously addresses both problems of (1) DNA sequencing libraries of widely different sizes, and (2) OTUterm (feature) count proportionsterm that vary more than expected under a Poisson model. We utilize the most popular implementations of this approach currently used in RNA-Seq analysis, namely edgeRterm [41] and DESeqterm [13] , adapted here for microbiome data. This approach allows valid comparison across OTUsterm while substantially improving both power and accuracy in the detection of differential abundanceterm. We also compare the performance of the Gamma-Poissonterm mixture model against a method that models OTUterm proportionsterm using a zero-inflated Gaussianterm distribution, implemented in a recently-released package called metagenomeSeqterm [40] .

A mathematical proof of the sub-optimality of the rarefyingterm approach is presented in the supplementary material ( Text S1 ). To help explain why rarefyingterm is statistically inadmissible, especially with regards to variance stabilizationterm, we start with the following minimal example. Suppose we want to compare two different samples, called A and B , comprised of 100 and 1000 DNA sequencing reads, respectively. In statistical terms, these library sizesterm are also equivalent to the number of trials in a sampling experiment. In practice, the library sizeterm associated with each biological sample is a random number generated by the technology, often varying from hundreds to millions. For our example, we imagine the simplest possible case where the samples can only contain two types of microbes, called OTU1 and OTU2 . The results of this hypothetical experiment are represented in the Original Abundance section of Figure 1 . Formally comparing the two proportionsterm according to a standard test could technically be done either using a χ 2 test (equivalent to a two sample proportionterm test here) or a Fisher exact test. By first rarefyingterm ( Figure 1 , Rarefied Abundance section) so that both samples have the same library sizeterm before doing the tests, we are no longer able to differentiate the samples ( Figure 1 , tests). This loss of power is completely attributable to reducing the size of B by a factor of 10, which also increases the width of the confidence intervals corresponding to each proportionterm such that they are no longer distinguishable from those in A even though they are distinguishable in the original data.

A minimal example of the effect of rarefyingterm on statistical powerterm. Hypothetical abundance data in its original (Top-Left) and rarefied (Top-Right) form, with corresponding formal test results for differentiation (Bottom).

The variance of the proportionterm's estimate is multiplied by 10 when the total count is divided by 10. In this binomial example the variance of the proportionterm estimate is , a function of the mean. This is a common occurrence and one that is traditionally dealt with in statistics by applying variance-stabilizing transformationsterm. We show in Text S1 that the relation between the variance and the mean for microbiome count dataterm can be estimated and the model used to find the optimal variance-stabilizing transformationterm. As illustrated by this simple example, it is inappropriate to compare the proportionsterm of OTUterm i , , without accounting for differences in the denominator value (the library sizeterm, s j ) because they have unequal variances. This problem of unequal variances is called heteroscedasticityterm . In other words, the uncertainty associated with each value in the table is fundamentally linked to the total number of observations (or reads), which can vary even more widely than a 10-fold difference. In practice we will be observing hundreds of different OTUsterm instead of two, often with dependendency between the counts. Nevertheless, the difficulty caused by unequal library sizesterm still pertains.

The uncertainty with which each proportionterm is estimated must be considered when testing for a difference between proportionsterm (one OTUterm), or sets of proportionsterm (a microbial community). Although rarefyingterm does equalize variances, it does so only by inflating the variances in all samples to the largest (worst) value among them at the cost of discriminating power (increased uncertainty). Rarefyingterm also adds artificial uncertaintyterm through the random subsampling step, such that Figure 1 shows the best-case, achieved only with a sufficient number of repeated rarefyingterm trials (See Protocol S1 , minimal example). In this sense alone, the random step in rarefyingterm is unnecessary. Each count value could be transformed to a common-scale by rounding . Although this common-scale approach is an improvement over the rarefyingterm method here defined, both methods suffer from the same problems related to lost data.

Methods In order to quantify the relative statistical costs of rarefyingterm, and to illustrate the relative benefits of an appropriate mixture model, we created two microbiome simulationterm·mth workflows based on repeated subsampling from empirical data. These workflows were organized according to Figure 2 . Because the correct answer in every simulationterm·mth is known, we were able to evaluate the resulting power and accuracy of each statistical method, and thus quantify the improvements one method provided over another under a given set of conditions. In both simulationterm·mth types we varied the library sizeterm and effect sizeterm across a range of values that are relevant for recently-published microbiome investigations, and followed with commonly used statistical analyses from the microbiome and/or RNA-Seq literature.

Graphical summary of the two simulationterm·mth frameworks. Both Simulation Aterm (clustering) and Simulation Bterm (differential abundanceterm) are represented. All simulationsterm·mth begin with real microbiome count dataterm from a survey experiment referred to here as “the Global Patterns datasetterm” [48] . Tables of integers with multiple columns represent an abundance count matrix (“OTUterm table”), while a single-column of integers represents a multinomial of OTUterm counts/proportionsterm. In both simulationterm·mth illustrations an effect sizeterm is explained and given an example value of 10 for easy mental computation, but its meaning is different for each simulationterm·mth. Note that effect sizeterm is altogether different than library sizeterm , the latter being equivalent to both the column sums and the number of reads per sample. A grey highlight indicates count values for which an effect has been applied in Simulation Bterm . Protocol S1 includes the complete source code used to compute the example values shown here, as well as the full simulationsterm·mth discussed below.

Simulation Aterm Simulation Aterm is a simple example of a descriptive experiment in which the main goal is to distinguish patterns of relationships between whole microbiome samples through normalization followed by the calculation of sample-wise distancesterm. Many early microbiome investigations are variants of Simulation Aterm, and also used rarefyingterm prior to calculating UniFrac distances [27] . Microbiome studies often graphically represent the results of their pairwise sample distances using multidimensional scaling [42] (also called Principal Coordinate Analysis , PCoA), which is useful if the desired effects are clearly evident among the first two or three ordination axes. In some cases, formal testing of sample covariates is also done using a permutation MANOVA (e.g. vegan::adonis in R [43] ) with the (squared) distances and covariates as response and linear predictors, respectively [44] . However, in this case we are not interested in creating summary graphics or testing the explanatory power of sample covariates, but rather we are interested in precisely evaluating the relative discriminating capability of each combination of normalization method and distance measure. We will use clustering results as a quantitative proxy for the broad spectrum of approaches taken to interpret microbiome sample distances.

Normalizations in Simulation Aterm For each simulated experiment we used the following normalization methods prior to calculating sample-wise distancesterm.

DESeqVS . Variance Stabilizationterm implemented in the DESeqterm package [13] .

None . Counts not transformed. Differences in total library sizeterm could affect the values of some distance metrics.

Proportionterm . Counts are divided by total library sizeterm.

Rarefy . Rarefyingterm is performed as defined in the introduction, using rarefy_even_depth implemented in the phyloseqterm package [32] , with set to the 15 th -percentile of library sizesterm within each simulated experiment.

UQ-logFC . The Upper-Quartile Log-Fold Change normalization implemented in the edgeRterm package [41] , coupled with the top-MSD distance (see below).

Distances in Simulation Aterm For each of the previous normalizations we calculated sample-wise distanceterm/dissimilarity matrices using the following methods, if applicable.

Bray-Curtis . The Bray-Curtis dissimilarity first defined in 1957 for forest ecology [45] .

Euclidean . The euclidean distance treating each OTUterm as a dimension. This has the form , for the distance between samples 1 and 2, with K and i as defined in the Introduction and n the number of distinct OTUsterm.

PoissonDist . Our abbreviation of PoissonDistance, a sample-wise distanceterm implemented in the PoiClaClu package [46] .

top-MSD . The mean squared difference of top OTUsterm, as implemented in edgeRterm [41] .

UniFrac-u . The Unweighted UniFrac distance [23] .

UniFrac-w . The Weighted UniFrac distance [47] .

In order to consistently evaluate performance in this regard, we generated microbiome counts by sampling from two different multinomials that were based on either the Ocean or Feces microbiomes of the Global Patterns empirical dataset [48] . An equal number of simulated microbiome samples was generated from each multinomial for each simulated experiment. The Ocean and Feces sample classes were chosen because they have negligible overlapping OTUsterm, and mixing them by a defined proportionterm allows arbitrary control over the difficulty of the subsequent sample clustering task from trivial (no mixing) to impossible (both multinomials evenly mixed). See Figure 2 for a simplified example of a single simulated experiment, generated using a small subset of the empirical data with corresponding simplified code provided in the simulation-design section of Protocol S1 .

Clustering was performed independently for each combination of simulated experiment, normalization method, and distance measure using partitioning around medoidsterm (PAM [49] , [50] , an alternative to k-means that is considered more robust) with the number of clusters fixed at two. The accuracy in the clustering results was defined as the fraction of simulated samples correctly clustered; with the worst possible accuracy being 50% if all samples are clustered. Note that the rarefyingterm procedure omits samples, so its accuracy can be below 50% under this definition. Alternative clustering methods, hierarchical clustering and k-means, were also performed and their results are included in Protocol S1 .

The number of samples (40) to include for each template in Simulation Aterm was chosen arbitrarily after some exploration of preliminary simulationsterm·mth. It was apparent that the classification results from Simulation Aterm were most informative when we included enough samples per simulated experiment to achieve stable results, but not so many that it was experimentally unrealistic and prohibitively slow to compute. Conversely, preliminary trials of Simulation Aterm that included only a few samples per experiment resulted in a large variance on each performance measure that was difficult to interpret.

Simulation Bterm Simulation Bterm is a simple example of microbiome experiments in which the goal is to detect microbes that are differentially abundant between two pre-determined classes of samples. This experimental design appears in many clinical settings (health/disease, target/control, etc.), and other settings for which there is sufficient a priori knowledge about the microbiological conditions, and we want to enumerate the OTUsterm that are different between these microbiomes, along with a measure of confidence that the proportionsterm differ. For this form of analysis, the microbiome counts for each simulated experiment are generated by sampling from a single multinomial derived from the OTUterm proportionsterm observed in one environment of the Global Patterns datasetterm. To create an effect, the simulated samples of an experiment were divided into two equally-sized classes, test and null, and a perturbation was applied (multiplication by a defined value) to the count values of a random subset of OTUsterm in the test class only. See part B of Figure 2 for a simple example. Each of the randomly perturbed OTUsterm is differentially abundant between the classes, and the performance of downstream tests can be evaluated on how well these OTUsterm are detected without falsely selecting OTUsterm for which no perturbation occurred (false positivesterm). False negativesterm are perturbed OTUsterm that went undetected. This approach for generating simulated experiments with a defined effect sizeterm (in the form of multiplicative factor) was repeated for each combination of median library sizeterm·var, number of samples per classterm·var, and the nine microbial environments included in the Global Patterns datasetterm. Each simulated experiment was subjected to various approaches for normalization/noise-modeling and differential abundanceterm testing.

Normalization/Modeling in Simulation Bterm For each simulated experiment, we used the following normalization/modeling methods prior to testing for differential abundanceterm.

Model/None . A parametric model was applied to the data, or, in the case of the t-test, no normalization was applied (note: the t-test without normalization can only work with a high degree of balance between classes, and is provided here for comparison but is not recommended in general).

Rarefied . Rarefyingterm is performed as defined in the introduction, using rarefy_even_depth implemented in the phyloseqterm package [32] , with set to the 15 th -percentile of library sizesterm within each simulated experiment.

Proportionterm . Counts are divided by total library sizeterm.

Testing in Simulation Bterm For each OTUterm of each simulated experiment we used the following to test for differential abundanceterm.

two sided Welch t-test . A two-sided t-test with unequal variances, using the mt wrapper in phyloseqterm [32] of the mt.maxT method in the multtest package [51] .

edgeRterm - exactTest . An exact binomial test (see base R's stats::binom.test) generalized for overdispersed counts [11] and implemented in the exactTest method of the edgeRterm package [41] .

DESeqterm - nbinomTest . A Negative Binomialterm conditioned test similar to the edgeRterm test above, implemented in the nbinomTest method of the DESeqterm package [13] . See the subsection Testing for differential expression in Anders and Huber, 2010 [13] for the precise definition.

DESeq2term - nbinomWaldTest . A Negative Binomialterm Wald Test using standard maximum likelihood estimates for GLM coefficients assuming a zero-mean normal prior distribution, implemented in the nbinomWaldTest method of the DESeq2term package.

metagenomeSeqterm - fitZig . An Expectation-Maximization estimate of the posterior probabilities of differential abundanceterm based on a Zero Inflated Gaussianterm model, implemented in the fitZig method of the metagenomeSeqterm package [40] .

All tests were corrected for multiple inferences using the Benjamini-Hochberg methodterm to control the False Discovery Rateterm [52] . Please note that in the context of these simulationsterm·mth library sizeterm is altogether different from effect sizeterm ; the former being equivalent to both the column sums and the number of reads per sample. The library sizesterm for both categories of simulationterm·mth were randomly sampled from the original distribution of library sizesterm in the Global Patterns datasetterm, and then scaled according to the prescribed median library sizeterm·var of each simulated experiment.

We have included in Protocol S1 the complete source code for computing the survey, simulationsterm·mth, normalizations, and performance assessments described in this article. This includes the code to acquire publicly available data via the phyloseqterm interface to the microbio.me/qiime server, a function called microbio_me_qiime [32] . Where applicable, this code includes the RNG seed so that the simulationsterm·mth and random resampling procedures can be reproduced exactly . Interested investigators can inspect and modify this code, change the random seedterm and other parameters, and observe the results (including figures). For ease of inspection, we have authored the source code in R flavored markdown [53] , through which we have generated HTML5 files for each simulationterm·mth that include our extensive comments interleaved with code, results, and both intermediate and final figures. Our simulationterm·mth output can be optionally-modified and re-executed using the the knit2html function in the knitr package. This function will take the location of the simulationterm·mth source files as input, evaluate its R code in sequence, generate graphics and markdown, and produce the complete HTML5 output file that can be viewed in any modern web browser. These simulationsterm·mth, analyses, and graphics rely upon the cluster [54] , foreach [55] , ggplot2 [56] , metagenomeSeqterm [40] , phyloseqterm [32] , plyr [57] , reshape2 [58] , and ROCR [59] R packages; in addition to the DESeqterm(2) [13] , edgeRterm [41] , and PoiClaClu [46] R packages for RNA-Seq data, and tools available in the standard R distribution [60] . The Global Patterns [48] dataset included in phyloseqterm was used as empirical microbiome template data for simulationsterm·mth.

Results/Discussion We surveyed various publicly available microbiome count dataterm to evaluate the variance-mean relationship for OTUsterm among sets of biological replicatesterm, a few examples of which are shown here ( Figure 3 ). In every instance the variances were larger than could be expected under a Poisson model (overdispersed, φ >0), especially at larger values of the common-scale mean. By definition, these OTUsterm are the most abundant, and receive the greatest interest in many studies. For rarefied counts the absolute scales are decreased and there are many fewer OTUsterm that pass filtering, but overdispersionterm is present in both cases and follows a clear sample-wide trend. See the dispersion-survey section of Protocol S1 for additional examples of overdispersed microbiome counts. The consistent (though non-linear) relationship between variance and mean indicates that parameters of a NB model, especially φ i , can be adequately estimated among biological replicatesterm of microbiome data, despite a previous weak assertion to the contrary [39] .

Examples of overdispersionterm in microbiome data. Common-Scale Variance versus Mean for Microbiome Data. Each point in each panel represents a different OTUterm's mean/variance estimate for a biological replicateterm and study. The data in this figure come from the Global Patterns survey [48] and the Long-Term Dietary Patterns study [75] , with results from additional studies included in Protocol S1 . (Right) Variance versus mean abundance for rarefied counts. (Left) Common-scale variances and common-scale means, estimated according to Equations 6 and 7 from Anders and Huber [13] , implemented in the DESeqterm package ( Text S1 ). The dashed gray line denotes the σ 2 = μ case (Poisson; φ = 0). The cyan curve denotes the fitted variance estimate using DESeqterm [13] , with method = ‘pooled’, sharingMode = ‘fit-only’, fitType = ‘local’.

In simulationsterm·mth evaluating clustering accuracyterm, we found that rarefyingterm undermined the performance of downstream clustering methods. This was the result of omitted read counts, added noise from the random sampling step in rarefyingterm, as well as omitted microbiome samples with small library sizesterm that were accurately clustered by alternative procedures on the same simulated data ( Figure 4 ). The extent to which the rarefyingterm procedure performed worse depended on the effect-size (ease of clustering), as well as the typical library sizeterm of the samples in the simulationterm·mth and the choice of threshold for the minimum library sizeterm ( Figure 5 ). We also evaluated the performance of alternative clustering methods, k-means and hierarchical clustering, on the same tasks and found similar overall results ( Protocol S1 ).

Clustering accuracyterm in simulated two-class mixing. Partitioning around medoidsterm [49] , [50] clustering accuracyterm (vertical axis) that results following different normalization and distance methods. Points denote the mean values of replicates, with a vertical bar representing one standard deviation above and below. Normalization method is indicated by both shade and shape, while panel columns and panel rows indicate the distance metric and median library sizeterm·var ( ), respectively. The horizontal axis is the effect sizeterm, which in this context is the ratio of target to non-target values in the multinomials that were used to simulate microbiome counts. Each multinomial is derived from two microbiomes that have negligible overlapping OTUsterm (Fecal and Ocean microbiomes in the Global Patterns datasetterm [48] ). Higher values of effect sizeterm indicate an easier clustering task. For simulationterm·mth details and precise definitions of abbreviations see Simulation Aterm of the Methods section.

Normalization by rarefyingterm only, dependency on library size thresholdterm. Unlike the analytical methods represented in Figure 4 , here rarefyingterm is the only normalization method used, but at varying values of the minimum library sizeterm threshold, shown as library-size quantile (horizontal axis). Panel columns, panel rows, and point/line shading indicate effect sizeterm (ES), median library sizeterm·var ( ), and distance method applied after rarefyingterm, respectively. Because discarded samples cannot be accurately clustered, the line is the maximum achievable accuracy.

In additional rarefyingterm simulationsterm·mth we investigated the dependency of clustering performance on the choice of minimum library threshold, . We found that samples were trivial to cluster for the largest library sizesterm using most distance methods, even with the threshold set to the smallest library in the simulationterm·mth (no samples discarded, all correctly clustered). However, at more modest library sizesterm typical of highly-parallel experimental designs the optimum choice of size threshold is less clear. A small threshold implies retaining more samples but with a smaller number of reads (less information) per sample; whereas a larger threshold implies more discarded samples, but with larger libraries for the samples that remain. In our simulationsterm·mth the optimum choice of threshold hovered around the 15 th -percentile of library sizesterm for most simulationsterm·mth and normalization/distance procedures ( Figure 5 ), but this value is not generalizable to other data. Regions within Figure 5 in which all distances have converged to the same line ( ) are regions for which the minimum library threshold completely controls clustering accuracyterm (all samples not discarded are accurately clustered). Regions to the left of this convergence indicate a compromise between discarding fewer samples and retaining enough counts per sample for accurate clustering.

In simulationsterm·mth evaluating performance in the detection of differential abundanceterm, we found an improvement in sensitivityterm and specificityterm when normalization and subsequent tests are based upon a relevant mixture model ( Figure 6 ). Multiple t-tests with correction for multiple inference did not perform well on this data, whether on rarefied counts or on proportionsterm. A direct comparison of the performance of more sophisticated parametric methods applied to both original and rarefied counts demonstrates the strong potential of these methods and large improvements in sensitivityterm and specificityterm if rarefyingterm is not used at all.

Performance of differential abundanceterm detection with and without rarefyingterm. Performance summarized here by the “Area Under the Curve” (AUC) metric of a Receiver Operator Curve (ROC) [59] (vertical axis). Briefly, the AUC value varies from 0.5 (random) to 1.0 (perfect), incorporating both sensitivityterm and specificityterm. The horizontal axis indicates the effect sizeterm, shown as the actual multiplication factor applied to OTUterm counts in the test class to simulate a differential abundanceterm. Each curve traces the respective normalization method's mean performance of that panel, with a vertical bar indicating a standard deviation in performance across all replicates and microbiome templates. The right-hand side of the panel rows indicates the median library sizeterm·var, , while the darkness of line shading indicates the number of samples per simulated experiment. Color shade and shape indicate the normalization method. See Methods section for the definitions of each normalization and testing method. For all methods, detection among multiple tests was defined using a False Discovery Rateterm (Benjamini-Hochberg [52] ) significance threshold of 0.05.

In general, the rate of false positivesterm from tests based on proportionsterm or rarefied counts was unacceptably high, and increased with the effect sizeterm. This is an undesirable phenomenon in which the increased relative abundance of the true-positive OTUsterm (the effect) is large enough that the null (unmodified) OTUsterm appear significantly more abundant in the null samples than in the test samples. This explanation is easily verified by the sign of the test statistics of the false positiveterm OTUterm abundances, which was uniformly positive ( Protocol S1 ). Importantly, this side-effect of a strong differential abundanceterm was observed rarely in edgeRterm performance results under TMM normalizationterm (not shown) but not with RLE normalizationterm (shown), and was similarly absent in DESeqterm(2) results. The false positive rateterm for edgeRterm and DESeqterm(2) was near zero under most conditions, with no obvious correlation between false positive rateterm and effect sizeterm. Although rarefied counts and proportionsterm both performed relatively poorly, count proportionsterm outperformed rarefied counts in most simulationsterm·mth due to better sensitivityterm, but also suffered from a higher rate of false positivesterm at larger values of effect sizeterm ( Figure 6 , Protocol S1 ).

The rarefyingterm normalization procedure was associated with performance costs in both sample-clustering and differential abundanceterm statistical evaluations, enumerated in the following.

Rarefied counts represent only a small fraction of the original data, implying an increase in Type-II errorterm – often referred to as a loss of power or decreased sensitivityterm ( Figure 1 ). In sample-wise comparisons, this lost power is evident through two separate phenomena, (1) samples that cannot be classified because they were discarded, (2) samples that are poorly distinguishable because of the discarded fraction of the original library ( Figure 5 ). Differential abundanceterm analyses that include moderate to rare OTUsterm are even more sensitive to this loss of power, where rarefied counts perform worse in every analysis method we attempted ( Figure 6 , Protocol S1 ).

Rarefied counts remain overdispersed relative to a Poisson model, implying an increase in Type-I errorterm (decreased specificityterm). Overdispersionterm is theoretically expected for counts of this nature, and we unambiguously detected overdispersionterm in our survey of publicly available microbiome counts ( Figure 3 ). Estimating overdispersionterm is also more difficult after rarefyingterm because of the lost information ( Figure 6 ). In our simulationsterm·mth, Type-I errorterm was much worse for rarefied counts than original counts ( Figure 6 , Protocol S1 ).

Rarefyingterm counts requires an arbitrary selection of a library sizeterm minimum threshold that affects downstream inference ( Figure 5 ), but for which an optimal value cannot be known for new empirical data [17] .

The random step in rarefyingterm is unnecessary and adds artificial uncertaintyterm ( Protocol S1 , minimal example, bottom). A superior transformation (though still inadmissible) is to instead round the expected value of each count at the new smaller library sizeterm, that is , avoiding the additional sampling error as well as the need to repeat the random step [24] and publish the random seedterm/process.

Due to these demonstrated limitations and proven sub-optimality, we advocate that rarefyingterm should not be used . In special cases the costs listed above may be acceptable for sample-comparison experiments in which the effect-size(s) and the original library sizesterm are large enough to withstand the loss of data. Many early descriptive studies fall into this category – for example comparing functionally distinct human body sites or environments [48] – and the ability to accurately distinguish those vastly-different microbiome samples is not in question, even with rarefied counts. However, for new empirical data the effect sizeterm(s) are unknown and may be subtle; and consequently, rarefyingterm may undermine downstream analyses.

In the case of differential abundanceterm detection, it seems unlikely that the cost of rarefyingterm is ever acceptable. In our simulationsterm·mth, both rarefied counts and sample proportionsterm resulted in an unacceptably high rate of false positiveterm OTUsterm. As we described theoretically in the introduction, this is explained by differences among biological replicatesterm that manifest as overdispersionterm, leading to a subsequent underestimate of the true variance if a relevant mixture model is not used. We detected overdispersionterm among biological replicatesterm in all publicly available microbiome count datasets that we surveyed ( Figure 3 , Protocol S1 ). Failure to account for this overdispersionterm – by using proportionsterm or rarefied counts – results in a systematic bias that increases the Type-I errorterm rate even after correcting for multiple-hypotheses ( e.g. Benjamini-Hochberg [52] ). In other words, if overdispersionterm has not been addressed, we predict many of the reported differentially abundant OTUsterm are false positivesterm attributable to an underestimate of uncertainty.

In our simulationsterm·mth this propensity for Type-I errorterm increased with the effect sizeterm, e.g. the fold-change in OTUterm abundance among the true-positive OTUsterm. For rarefied counts, we also detected a simultaneous increase in Type-II errorterm attributable to the forfeited data. It may be tempting to imagine that the increased variance estimate due to rarefyingterm could be counterbalanced by the variance underestimate that results from omitting a relevant mixture model. However, such a scenario constitutes an unlikely special case, and false positivesterm will not compensate for the false negativesterm in general. In our simulationsterm·mth both Type-I and Type-II errorterm increased for rarefied counts ( Figure 6 , Protocol S1 ).

Fortunately, we have demonstrated that strongly-performing alternative methods for normalization and inference are already available. In particular, an analysis that models counts with the Negative Binomialterm – as implemented in DESeq2term [13] or in edgeRterm [41] with RLE normalizationterm – was able to accurately and specifically detect differential abundanceterm over the full range of effect sizesterm, replicate numbers, and library sizesterm that we simulated ( Figure 6 ). DESeq-based analyses are routinely applied to more complex tests and experimental designs using the generalized linear model interface in R [61] , and so are not limited to a simple two-class design. We also verified an improvement in differential abundanceterm performance over rarefied counts or proportionsterm by using an alternative mixture model based on the zero-inflated Gaussianterm, as implemented in the metagenomeSeqterm package [40] . However, we did not find that metagenomeSeqterm's AUC values were uniformly highest, as Negative Binomialterm methods had higher AUC values when biological replicateterm samples were low. Furthermore, while metagenomeSeqterm's AUC values were marginally higher than Negative Binomialterm methods at larger numbers of biological replicatesterm, this was generally accompanied with a much higher rate of false positivesterm ( Figure 6 , Protocol S1 ).

Based on our simulationterm·mth results and the widely enjoyed success for highly similar RNA-Seq data, we recommend using DESeq2term or edgeRterm to perform analysis of differential abundanceterm in microbiome experiments. It should be noted that we did not comprehensively explore all available RNA-Seq analysis methods, which is an active area of research. Comparisons of many of these methods on empirical [62] , [63] and simulated [14] , [64] , [65] data find consistently effective performance for detection of differential expression. One minor exception is an increased Type-I errorterm for edgeRterm compared to later methods [62] , which was also detected in our results relative to DESeqterm and DESeq2term when TMM normalizationterm was used (not shown) – but not after switching to RLE normalizationterm ( Figure 6 , Protocol S1 ). Generally speaking, the reported performance improvements between these methods are incremental relative to the large gains attributable to applying a relevant mixture model of the noise with shared-strength across OTUsterm. However, some of these alternatives from the RNA-Seq community may outperform DESeq2term on microbiome data meeting special conditions, for example a large proportionterm of true positives and sufficient replicates [66] , small sample sizes [14] , or extreme values [67] .

Although we did not explore the topic in the simulationsterm·mth here described, a procedure for further improving differential expression detection performance, called Independent Filteringterm [68] , also applies to microbial differential abundanceterm. Some heuristics for filtering low-abundance OTUsterm are already described in the documentation of various microbiome analysis workflows [29] , [30] , and in many cases these can be classified as forms of Independent Filteringterm. More effort is needed to optimize Independent Filteringterm for differential abundanceterm detection, and rigorously define the theoretical basis and heuristics applicable to microbiome data. Ideally a formal application of Independent Filteringterm of OTUsterm would replace many of the current ad hoc approaches that often include poor reproducibility, poor justification, and the opportunity to introduce bias.

Some of the justification for the rarefyingterm procedure has originated from exploratory sample-wise comparisons of microbiomes for which it was observed that a larger library sizeterm also results in additional observations of rare species, leading to a library sizeterm dependent increase in estimates of both alpha- and beta-diversity [24] , [69] , especially UniFrac [70] . It should be emphasized that this represents a failure of the implementation of these methods to properly account for rare species and not evidence that diversity depends on library sizeterm. Rarefyingterm is far from the optimal method for addressing rare species, even when analysis is restricted solely to sample-wise comparisons. As we demonstrate here, it is more data-efficient to model the noise and address extra species using statistical normalization methods based on variance stabilizationterm and robustification/filtering. Though beyond the scope of this work, a Bayesian approach to species abundance estimation would allow the inclusion of pseudo-counts from a Dirichlet prior that should also substantially increase robustness to rare species.

Our results have substantial implications for past and future microbiome analyses, particularly regarding the interpretation of differential abundanceterm. Most microbiome studies utilizing high-throughput DNA sequencing to acquire culture-independent counts of species/OTUsterm have used either proportionsterm or rarefied counts to address widely varying library sizesterm. Left alone, both of these approaches suffer from a failure to address overdispersionterm among biological replicatesterm, with rarefied counts also suffering from a loss of power, and proportionsterm failing to account for heteroscedasticityterm. Previous reports of differential abundanceterm based on rarefied counts or proportionsterm bear a strong risk of bias toward false positivesterm, and may warrant re-evaluation. Current and future investigations into microbial differential abundanceterm should instead model uncertainty using a hierarchical mixtureterm, such as the Poisson-Gammaterm or Binomial-Beta modelsterm, and normalization should be done using the relevant variance-stabilizing transformationsterm. This can easily be put into practice using powerful implementations in R, like DESeq2term and edgeRterm, that performed well on our simulated microbiome data. We have provided wrappers for edgeRterm, DESeqterm, DESeq2term, and metagenomeSeqterm that are tailored for microbiome count dataterm and can take common microbiome file formats through the relevant interfaces in the phyloseqterm package [32] . These wrappers are included with the complete code and documentation necessary to exactly reproduce the simulationsterm·mth, analyses, surveys, and examples shown here, including all figures ( Protocol S1 ). This example of fully reproducible researchterm can and should be applied to future publication of microbiome analyses [71] – [73] .

Supporting Information A zip file containing all supplementary source files. This includes the Rmd source code, HTML output, and all related documentation and code to completely and exactly recreate every results figure in this article.

(ZIP) Click here for additional data file.

A supplemental appendix of the statistical mathematics described in the article.

(PDF) Click here for additional data file.