High-throughput (HT) DNA sequencing [1] is allowing major advances in microbial ecology studies [2], where our understanding of the presence and abundance of microbial species relies heavily on the observation of their nucleic acids in a “culture independent” manner [3]. This nucleic acid sequencing based census of the inhabitants of microbiome samples is very often now accompanied with other experimental observations (e.g. clinical, environmental, metabolomic, etc.), in addition to phylogenetic tree reconstruction and/or taxonomic classification of the sequences. Here we refer to this as “phylogenetic sequencing” data if it can be usefully represented as a contingency table of taxonomic units and samples, and integrated with the other aforementioned data types. Importantly, this term – also the namesake of the software here described – is defined so as to not be specific to the method by which the phylogenetically relevant microbial census data was obtained, reflecting the intended level of data abstraction in the software. The following are two examples of common methods for producing phylogenetic sequencing data.

Barcoded [2] amplicon sequencing of dozens to hundreds of samples [4] is a method of phylogenetic sequencing of microbiomes, often targeting the small subunit ribosomal RNA (16S rRNA) gene [3], for which there are also convenient tools [5] and large reference databases [6]–[8]. The task of decoding the sample source of each sequence read by its barcode, followed by similarity clustering to define operational taxonomic units (OTUs, sometimes referred to as taxa) [9], [10] can be performed by publicly available packages/pipelines, including QIIME [11], mothur [12], and PANGEA [13]; as well as virtual machine (VM) and cloud-based solutions such as the RDP pipeline [7], Pyrotagger [14], CLoVR-16S [15], Genboree [16], QIIME EC2 image [17], n3phele [18], and MG-RAST [19].

An alternative experimental method is random “shotgun” sequencing [20], [21] of un-amplified metagenomic DNA [22], in which case OTU clustering and counting is based upon one or more detectable phylogenetic markers in the metagenomic sequence fragments, using tools such as phylOTU [23]. It is worth noting that bias from PCR amplification is avoided in this latter approach – at the expense of per-sequence efficiency [23] – and both methods are now commonly used for phylogenetic sequencing (Figure 1).

Many of the previously mentioned OTU-clustering applications also perform additional downstream analyses (File S1). However, typically an investigator must port the human-unreadable output data files to other software for additional processing and statistical analysis specific to the goals of the investigation. The powerful statistical, ecological, and graphics tools available in R [24] make it an attractive option for this post-clustering stage of analysis. While the computational efficiency of compiled languages like [25] make them appropriate for the expensive but well-defined requirements of the initial sequence-processing, the subsequent analysis is vaguely-defined and project specific; requiring instead a broad set of interactive calculations that is often less computationally expensive and for which R is well-suited [26]. The public repositories of open-source R extensions (“packages” or “libraries”) include many dedicated ecology and phylogenetic packages. For instance, there are several dozen packages listed in the CRAN Ecology Task View [27], as well as distory [28], phangorn [29], picante [30], and now phyloseq [31]. Furthermore, R includes infrastructure for documenting an analysis in such a way that it can be easily reproduced and modified by peers [32], [33].

In spite of all of these highly relevant tools, we recently described the lack of a satisfactory standard within Bioconductor [34] (or R generally) for importing the data files from the most popular OTU-clustering applications, or representing this data in a complete, integrated class [31]. One Bioconductor package, OTUbase [35], pursues some of these goals, but has no support for phylogenetic trees in its data class, nor support for importing data from popular/recent OTU-clustering output formats [35], [36] (File S1). We have proposed a new Bioconductor package, phyloseq (from “phylogenetic sequencing”), dedicated to the object-oriented representation and analysis of phylogenetic sequencing data in R [31], and supporting common OTU-clustering output formats like QIIME [11], mothur [12], the RDP-pipeline [7], Pyrotagger [14], and the biom-format [37].

In this article we describe the conceptual framework and toolbox of a substantially enhanced phyloseq codebase, including especially some advanced ordination and graphics capabilities. We further note that data imported by phyloseq is also accessible to analyses encoded by a large number of freely available R packages, in addition to the capabilities directly supported by phyloseq itself. We will end by discussing the notion of “reproducible research” in the context of phylogenetic sequencing data, and how phyloseq and R can be used in analyses that are more open and reproducible than those found in recent common practice.

The phyloseq package provides an object-oriented programming infrastructure that simplifies many of the common data management and preprocessing tasks required during analysis of phylogenetic sequencing data. This simplified syntax helps mitigate inconsistency errors and encourages interaction with the data during preprocessing. The phyloseq package also provides a set of powerful analysis and graphics functions, building upon related packages available in R and Bioconductor. It includes or supports some of the most commonly-needed ecology and phylogenetic tools, including a consistent interface for calculating ecological distances and performing dimensional reduction (ordination). The graphics functions allow users to interactively produce annotated publication-quality graphics in just one or two lines of code. The phyloseq package includes extensive documentation in the form of function- and package-level manuals embedded in the package's documentation interface and in a PDF version on Bioconductor [38], as well as extended reproducible examples on the phyloseq homepage [39], and open collaborative development on GitHub [40].

The phyloseq package adheres to the requirements for standard R packages set forth in the official “Writing R Extensions” manual [40]. It also satisfies additional requirements of the Bioconductor Repository [34], and uses a literate-programming framework based on structured in-source comments, called roxygen2 [41], for (re)building the R documentation (.Rd) files and the namespace specifications. The phyloseq package can be installed on any system on which R is supported, including Mac OS X, Windows, and most Linux distributions.

R packages can include example data that is documented with the same help system as other package objects [58]. This data becomes available in the R session by invoking the data function after the package has been loaded. Unless otherwise noted, the examples provided in this manuscript use example data that is included in the phyloseq package.

The phyloseq project includes an object-oriented class that integrates the heterogeneous components of OTU-clustered phylogenetic sequencing data. Although Bioconductor provides many utilities for efficient manipulation of DNA sequences, phyloseq does not currently re-implement any methods for DNA sequence decoding, processing, or OTU-clustering (Figure 1, File S1). Instead, phyloseq provides tools to read the output files of the most common OTU-clustering applications [7], [11], [12], [14], and represents this data in R as an instance of the main data class. This multi-component “experiment-level” class — named “phyloseq”, and referred to here as “the phyloseq-class” — is a key design feature of the phyloseq project, with subsequent user-accessible functions expecting to operate on an instance of this class as their sole or primary input data. These functions are described in detail in the phyloseq manual [7], and are part of a modular workflow summarized in Figure 2.

Figure 3 summarizes the structure of the phyloseq-class and its components. Each of the slots are empty (NULL) by default, although an instance missing an otu_table component is invalid. Tools in phyloseq that truncate dimensions of one component (that is, remove samples or OTUs) automatically propagate the change across all relevant components. In general, researchers only need to manipulate their “experiment-level” object, making data (pre)processing less prone to mistakes, and often simplifying analysis commands to just one data argument.

Complementing the data infrastructure, the phyloseq package provides a set of functions that take a phyloseq object as the primary data, and performs an analysis and/or graphics task. Figure 2 summarizes the general workflow within phyloseq, and lists some of the main functions/tools.

Comparisons of the type and quantity of OTUs observed between microbiome samples (“beta diversity”) is often approached through the calculation of pairwise ecological distances [42], [43], and through dimensional reduction (ordination) methods. The phyloseq package provides a consistent interface for the most common approaches to distance calculations and ordination. This interface is also the foundation for the custom ordination and heatmap graphics functions described in the next subsection.

In phyloseq the interface for ecological distance calculations is a single function, distance, that takes a phyloseq object as its data argument as well as a character string indicating the distance method, with explicit support for more than 40 ecological distance methods. This includes a R -native, optionally-parallel implementation of Fast UniFrac [44] (both weighted [45] and unweighted [46]). The output is a “dist” class distance matrix (lower-triangle) appropriate for standard clustering analysis in core R (e.g. hclust), as well as certain dimensional reduction (ordination) methods.

The interface for performing ordination methods is also a single function, called ordinate, that takes a phyloseq object as its primary data argument and a character string indicating the desired ordination method. For example, the following would perform (unconstrained) correspondence analysis on the included “Global Patterns” dataset [47].

The ordinate function currently supports correspondence analysis (CA) [48], constrained correspondence analysis (CCA) [49], detrended correspondence analysis (DCA) [50], redundancy analysis (RDA) [51], principal components analysis (PCA) [52], double principle coordinates analysis (DPCoA) [53], multidimensional scaling (MDS, PCoA) [54], and non-metric multidimensional scaling (NMDS) [55]. For CA, CCA, DCA, RDA, and DPCoA, the ordination is based upon an evaluation of abundance values (in the case of DPCoA, the patristic distances between OTUs on the phylogenetic tree is also used), but not an ecological distance. For MDS and NMDS, the ordinate function requires a pre-calculated distance matrix (“dist” object) or the name of a supported ecological distance method. For example, PCoA/MDS can be calculated on an unweighted UniFrac distance matrix [46], using the following command:

There are many combinations of approaches possible (even extending into time-series of table pairs), and the optimal approach depends on the goals of the experiment and characteristics of the data [56]. The phyloseq package also includes a specialized function for displaying ordination results in different ways, described in the following section.

One of the key features of the phyloseq package is a set of graphics functions custom-tailored for phylogenetic sequencing analysis, built using the ggplot2 package [57]. The ggplot2 package is an implementation of Wilkinson's The Grammar of Graphics, which provides an object-oriented description of analytical graphics that emphasizes the separation of data and its mapping to aesthetic attributes [58]. In the phyloseq package, functions having names beginning with “plot_” require a phyloseq object as input data, and return a ggplot2 graphics object. These plot_ functions support optional mapping of color, size, and shape aesthetics to sample or OTU variables — usually by providing the name of the variable or taxonomic rank as a character string (E.g. color = “SampleType”). Legends are automatically generated based on the data and aesthetic mappings (not true of the base R graphics), and all features of these graphics can be further modified in R via functions/options in the ggplot2 package.

The following list summarizes the key graphics-producing functions in phyloseq, which are also demonstrated in Figure 4, and in phyloseq's online tutorials [39]. File S2 provides the complete R code for creating Figures 4 and 5. We have also included some additional examples of graphics created by plot_ordination (Figure 5). They emphasize different aspects of ordination results, and the best choice depends heavily on characteristics of the data and research questions. The provided code also demonstrates a custom modification to the ggplot2 graphic, in this case the addition of a two-dimensional density estimate to the “OTUs-only” plot (File S2).

In multivariate analyses such as PCA, large differences in variances between columns are corrected by standardizing each column; i.e. dividing each column by its standard deviation. Thus each column will have the same weight in the multivariate analysis. For OTU abundance tables, such a procedure is inappropriate as the disparities in column sums can be 100-fold. Methods based on chi-squared distances rather than variances deal with this by comparing weighted column profiles [62], computed as relative abundances for each OTU within a column, with the overall column sum retained as a weighting factor. However, chi-square distances are sums of squares and can be overly sensitive to outliers and sequencing “jackpot” effects such as those occurring in pyrosequencing data [63]. Bray-Curtis distances can be a useful alternative, as it is based on the distance between profiles, as long as the differences in actual column sums are also accounted for in the final study. The other approach to the problem of disparities between column sums has been to subsample the over-abundant columns down to the same number as the smaller ones. However this results in a loss of information, rarely an optimal procedure in statistical contexts. This subsampling procedure is inspired by the popular idea of rarefaction in coverage studies first invented by Sanders [64], but has yet to be proved beneficial for all microbial community structures. The parallels between gene expression microarray analyses and microbial abundance analyses was mentioned in [65], which proposed several expression-inspired strategies for robustifying abundance measurements. The main points were that rankings and thresholding are important in the presence of noise and high variability in sequence depths. As in gene expression analysis filtering the OTUs is beneficial, especially in the latter multiple testing adjustments. The phyloseq package enables easy filtering and rank transformations in the same vein as robust multi-array averaging (rma) [66]. We provide further details in (McMurdie and Holmes, [67]).

Although useful for exploring and summarizing microbiome data, many of the graphics and ordination methods discussed here are not formal tests of any particular hypothesis. The most common framework for testing in microbiome studies is the comparisons of samples from different categories (e.g. healthy and obese; control and treated; different environments). Standard test statistics include the t-test, the paired permutation t-test, and ANOVA type tests based on F or pseudo-F statistics. However, microbiome data have two particularities. First, the raw abundance counts are never normally distributed, so the preferred methods are nonparametric. Second, there is contiguous information available about the relationships between OTUs, as well as for variables measured on the samples, so testing is sometimes more elaborate than a two-sample test. The hypergeometric test, also known as Fisher's exact test, is used in cases when we have a test statistic for each of the different OTUs. The goal is to confirm that a certain property of these significant OTUs is overrepresented compared to the general population of OTUs, often called “the universe”. For instance in Holmes et al [65] and Nelson et al [68] several phyla were shown to be significantly over-abundant in IBS rats as compared to healthy controls using this hypergeometric test.

An organizing principle in many nonparametric testing protocols is that the repetition of an analysis multiple times enables the user to control for multiple testing, or to evaluate the quality of estimators or the optimal values of tuning parameters. Modern confirmatory analyses currently depend on these repeated analyses under various data perturbation schemes, of which resampling, permutations, and Monte Carlo simulations are the most common. For instance the bootstrap uses many thousands of analyses of resampled data to address problems such as statistical stability or bias estimation [69], and can even provide confidence regions [69] for nonstandard parameters, such as phylogenetic trees [70]. Repeating analyses on permuted data can allow for control of the probability of encountering 1 or more false positives (falsely rejected nulls) among your group of simultaneous hypotheses, also called the Family Wise Error Rate (FWER). For instance, Westfall and Young's permutation-based minP procedure controls the FWER [71] and is implemented within the multtest package [72]. The phyloseq package interfaces with minP in multtest through a wrapper function, called mt. In the following example code we use the mt wrapper to control the FWER while simultaneously testing whether each OTU correlates with the “Enterotypes” classification of the samples. Note that we first remove samples that were not assigned an enterotype by the original authors (Table 1).

In exploratory statistical work it is easy to produce biased results [73] through poorly chosen metrics or tests, a failure to properly control for multiple inferences, undisclosed data “pruning”, and probably many other means. Although not commonly required [74]–[76], an important defense against biased (or poorly-supported) findings is a higher standard for reproducibility in published research [77], in which journal articles are accompanied by sufficient data and software such that all presented analyses, tables, and figures can be reproduced exactly and with minimal effort [75]. In this context of highly-parallel phylogenetic-sequencing experiments, reproducible research can be partially facilitated by emerging standards for experimental design [78] and file format [37]. Virtual machine image and cloud-deployed “pipeline” analyses [11], [15], [19] can further increase accessibility of analyses by mitigating the need for expensive computing hardware while also avoiding complicated installation procedures. However, the use of publicly available “pipeline” tools does not fully meet the reproducibility standard unless accompanied with the complete code and data used in the analysis being published [75]. This is especially important when considering the many choices that are involved in decoding, OTU-clustering, and preprocessing; as well as the many varied approaches to incorporating sample covariates and performing multivariate analyses on complex data. The recent release of the HMP data and multiple articles on the results from their analyses underscore this fact. Thresholding and noise filtering were done independently by each team, but no overall robustness study was performed [79]. Changes early in the analysis pipeline could have downstream effects that are now prohibitively difficult or impossible to evaluate. Generally speaking, the preprocessing of OTU abundance data through filtering, normalizing, centering, shrinking, and other transformations is a common practice and necessary for analysis [66], but varies widely among researchers and is often difficult to reproduce. This is particularly true when the preprocessing transformations are the result of “manual” adjustments in a spreadsheet, custom code/script that is not included in the publication, or random subsampling (“rarefying” to even sequencing effort) with no reported seed. A related example is the (often not-so) reproducible choice of tuning parameters and perturbation-based statistical validation procedures, allowing for the easy testing of alternatives and robustness of results. To a large extent this revisits many of the same issues of reproducible research [80]–[82] that have been addressed over the last decade for the analysis of microarray data [66], and for which there are many proven tools already available in Bioconductor/R. The emphasis of preprocessing tools in phyloseq is intended to decrease the extent to which these steps constitute opaque and idiosyncratic efforts by investigators, while making the results of different studies more comparable.

One of the goals of the phyloseq project is to help close the gap in reproducible research that presently exists between pipeline results and the additional analyses required by investigators. This can be achieved when phyloseq is used (possibly with other R packages) in conjunction with documentation tools such as Sweave [32], knitr [33], iPython [83] Notebook invoking the rmagic extension, or “ R flavored markdown” (RFM) [84]. The Sweave-format approach is part of the reproducible research standards strongly encouraged by the journal Biostatistics [81], as well as many disciplines related to statistics and bioinformatics [77], [85]. The recently-described RFM format and iPython Notebook can also work very well for cases where a web-browser is a satisfactory documentation delivery medium, with RFM being our preferred source format for publishing reproducible online tutorials with embedded code and figures (HTML5) [39], [86]. We emphasize that the benefits of reproducibility are not contingent on “pretty” code [87], and we encourage researchers in the field to make their code available even if they feel insecure about its programmatic elegance. As an illustrative example, we have made available the Sweave (.Rnw) and supporting files required to completely reproduce this article, including especially the complete source as an RFM file (.Rmd) with its associated output HTML file, both of which provide the preprocessing steps and graphics commands needed to exactly reproduce each figure (File S2). We have also published a GitHub repository dedicated to reproducible demonstrations of analyses with phyloseq [86].

It is important to note that the new phyloseq-class is a significant departure from the originally-proposed phyloseq-class structure [31], which used nested multiple inheritance and a naming convention. It was a valid approach in principle, but was an overly complex approach for the goal of representing a phylogenetic sequencing experiment as a single object. The updated phyloseq-class is simple to extend for developers and easy to explain to users (Figure 3). In general, the downstream analysis and plotting functions that might operate on an instance of the phyloseq-class do not need to (re)perform common validity checks because these checks are consolidated as part of the phyloseq-constructor method.

Analysis tools available in R but not explicitly wrapped in phyloseq are nevertheless available to users and developers via accessors and other data infrastructure tools. This leverages the fact that phyloseq data components are based on standard R data classes and easily used in other package settings in R. For example, we have included example code that illustrates the use of the bioenv function from the vegan package, starting with data represented by the phyloseq-class (See File S2 for code, and the phyloseq demo [86]). Similarly, as an open-source package in an open language/framework (R), phyloseq can be easily included at the relevant steps in pipelines, workbenches, and GUIs now under active development (E.g. ClovR [15], MG-RAST [19], QIIME [11], mcaGUI [88]). This represents a means for investigators with limited programming literacy to still benefit from some of the tools included in, or facilitated by, phyloseq.