  The interpretation of count data originating from the current generation of
DNA sequencing platforms requires special attention. In particular, the
per-sample library sizes often vary by orders of magnitude from the same
sequencing run, and the counts are overdispersed relative to a simple Poisson
model These challenges can be addressed using an appropriate mixture model that
simultaneously accounts for library size differences and biological
variability. This approach is already well-characterized and implemented for
RNA-Seq data in R packages such as edgeR and DESeq.
  We use statistical theory, extensive simulations, and empirical data to show
that variance stabilizing normalization using a mixture model like the negative
binomial is appropriate for microbiome count data. In simulations detecting
differential abundance, normalization procedures based on a Gamma-Poisson
mixture model provided systematic improvement in performance over crude
proportions or rarefied counts -- both of which led to a high rate of false
positives. In simulations evaluating clustering accuracy, we found that the
rarefying procedure discarded samples that were nevertheless accurately
clustered by alternative methods, and that the choice of minimum library size
threshold was critical in some settings, but with an optimum that is unknown in
practice. Techniques that use variance stabilizing transformations by modeling
microbiome count data with a mixture distribution, such as those implemented in
edgeR and DESeq, substantially improved upon techniques that attempt to
normalize by rarefying or crude proportions. Based on these results and
well-established statistical theory, we advocate that investigators avoid
rarefying altogether. We have provided microbiome-specific extensions to these
tools in the R package, phyloseq.
