phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data — TEI-tagged view

Paul J. McMurdiepersName, Susan HolmespersName (2013date). “phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data.” PLoS ONEtitle·j.

Object of study — the representation, analysis, visualization and reproducibility of microbiome census dataterm in Rterm

Research question — how to provide an integrated, extensible and reproducible post-clustering workflow for heterogeneous phylogenetic sequencing dataterm·def

Contribution — the open-source phyloseqterm package, its experiment-level data class, analysis interfaces, graphics and reproducible examples

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.

phyloseqterm: An Rterm Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Dataterm Abstract The analysis of microbial communities through DNA sequencing brings many challenges: the integration of different types of data with methods from ecology, genetics, phylogenetics, multivariate statistics, visualization and testing. With the increased breadth of experimental designs now being pursued, project-specific statistical analyses are often needed, and these analyses are often difficult (or impossible) for peer researchers to independently reproduce. The vast majority of the requisite tools for performing these analyses reproducibly are already implemented in Rterm and its extensions (packages), but with limited support for high throughput microbiome census dataterm.

Here we describe a software project, phyloseqterm, dedicated to the object-oriented representationterm and analysis of microbiome census dataterm in Rterm. It supports importing data from a variety of common formats, as well as many analysis techniques. These include calibration, filtering, subsetting, agglomeration, multi-table comparisons, diversity analysis, parallelized Fast UniFracterm, ordination methodsterm·var, and production of publication-quality graphics; all in a manner that is easy to document, share, and modify. We show how to apply functions from other Rterm packages to phyloseq-represented data, illustrating the availability of a large number of open source analysis techniques. We discuss the use of phyloseqterm with tools for reproducible researchterm, a practice common in other fields but still rare in the analysis of highly parallel microbiome census dataterm. We have made available all of the materials necessary to completely reproduce the analysis and figures included in this article, an example of best practices for reproducible researchterm.

The phyloseqterm project for Rterm is a new open-source software packageterm, freely available on the web from both GitHub and Bioconductorterm.

Introduction Phylogenetic Sequencingterm 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 observationsterm (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 sequencingterm” data if it can be usefully represented as a contingency tableterm of taxonomic unitsterm 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 dataterm·def.

Barcoded [2] amplicon sequencingterm of dozens to hundreds of samples [4] is a method of phylogenetic sequencingterm of microbiomes, often targeting the small subunit ribosomal RNA (16S rRNAterm) 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 unitsterm (OTUsterm, sometimes referred to as taxaterm ) [9] , [10] can be performed by publicly available packages/pipelines, including QIIMEterm [11] , mothurterm [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] , QIIMEterm EC2 image [17] , n3phele [18] , and MG-RAST [19] .

An alternative experimental method is random “shotgunterm” sequencing [20] , [21] of un-amplified metagenomic DNA [22] , in which case OTUterm 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 amplificationterm is avoided in this latter approach – at the expense of per-sequence efficiency [23] – and both methods are now commonly used for phylogenetic sequencingterm ( Figure 1 ).

Example of a phylogenetic sequencingterm workflow. A diagram of an experimental and analysis workflow for amplicon or shotgunterm phylogenetic sequencingterm. The intended role for phyloseqterm is indicated.

The phyloseqterm Project Many of the previously mentioned OTU-clustering applicationsterm 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 Rterm [24] make it an attractive option for this post-clustering stageterm 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 Rterm is well-suited [26] . The public repositories of open-source Rterm 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 phyloseqterm [31] . Furthermore, Rterm 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 Bioconductorterm [34] (or Rterm generally) for importing the data files from the most popular OTU-clustering applicationsterm, or representing this data in a complete, integrated class [31] . One Bioconductorterm package, OTUbaseterm [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 Bioconductorterm package, phyloseqterm (from “ phylo genetic seq uencing”), dedicated to the object-oriented representationterm and analysis of phylogenetic sequencing dataterm·def in Rterm [31] , and supporting common OTU-clustering output formats like QIIMEterm [11] , mothurterm [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 phyloseqterm codebase, including especially some advanced ordinationterm and graphics capabilities. We further note that data imported by phyloseqterm is also accessible to analyses encoded by a large number of freely available Rterm packages, in addition to the capabilities directly supported by phyloseqterm itself. We will end by discussing the notion of “reproducible researchterm” in the context of phylogenetic sequencing dataterm·def, and how phyloseqterm and Rterm can be used in analyses that are more open and reproducible than those found in recent common practice.

Methods phyloseqterm Project Key Features The phyloseqterm package provides an object-oriented programming infrastructureterm that simplifies many of the common data management and preprocessing tasks required during analysis of phylogenetic sequencing dataterm·def. This simplified syntax helps mitigate inconsistency errors and encourages interaction with the data during preprocessing. The phyloseqterm package also provides a set of powerful analysis and graphics functions, building upon related packages available in Rterm and Bioconductorterm. It includes or supports some of the most commonly-needed ecology and phylogenetic tools, including a consistent interface for calculating ecological distancesterm and performing dimensional reductionterm (ordinationterm). The graphics functions allow users to interactively produce annotated publication-quality graphics in just one or two lines of code. The phyloseqterm 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 Bioconductorterm [38] , as well as extended reproducible examples on the phyloseqterm homepage [39] , and open collaborative development on GitHub [40] .

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

Data Availability Rterm packages can include example data that is documented with the same help system as other package objects [58] . This data becomes available in the Rterm 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 phyloseqterm package.

Data Infrastructure and Design The phyloseqterm project includes an object-oriented class that integrates the heterogeneous components of OTU-clustered phylogenetic sequencing dataterm·def. Although Bioconductorterm provides many utilities for efficient manipulation of DNA sequences, phyloseqterm does not currently re-implement any methods for DNA sequence decoding, processing, or OTU-clustering ( Figure 1 , File S1). Instead, phyloseqterm provides tools to read the output files of the most common OTU-clustering applicationsterm [7] , [11] , [12] , [14] , and represents this data in Rterm as an instance of the main data class. This multi-component “experiment-level” class — named “phyloseqterm”, and referred to here as “the phyloseq-classterm·def” — is a key design feature of the phyloseqterm 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 phyloseqterm manual [7] , and are part of a modular workflow summarized in Figure 2 .

Analysis workflow using phyloseqterm. The workflow starts with the results of OTUterm clustering and independently-measured sample data (Input, top left), and ends at various analytic procedures available in Rterm for inference and validation. In between are key functions for preprocessing and graphics. Rounded rectangles and diamond shapes represent functions and data objects, respectively, further described in Figure 3 .

Figure 3 summarizes the structure of the phyloseq-classterm·def and its components. Each of the slots are empty (NULL) by default, although an instance missing an otu_tableterm component is invalid. Tools in phyloseqterm that truncate dimensions of one component (that is, remove samples or OTUsterm) 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.

The “phyloseqterm” class. The phyloseq classterm is an experiment-level data storage classterm defined by the phyloseqterm package for representing phylogenetic sequencing dataterm·def. Most functions in the phyloseqterm package expect an instance of this class as their primary argument. See the phyloseqterm manual [38] for a complete list of functions.

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

Comparisons of the type and quantity of OTUsterm observed between microbiome samples (“beta diversityterm”) is often approached through the calculation of pairwise ecological distancesterm [42] , [43] , and through dimensional reductionterm (ordinationterm) methods. The phyloseqterm package provides a consistent interface for the most common approaches to distanceterm calculations and ordinationterm. This interface is also the foundation for the custom ordinationterm and heatmap graphics functions described in the next subsection.

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

The interface for performing ordination methodsterm·var is also a single function, called ordinateterm, that takes a phyloseqterm object as its primary data argument and a character string indicating the desired ordination methodterm·var. For example, the following would perform (unconstrained) correspondence analysisterm on the included “Global Patterns” dataset [47] .

The ordinateterm function currently supports correspondence analysisterm (CA) [48] , constrained correspondence analysisterm (CCA) [49] , detrended correspondence analysisterm (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 scalingterm (NMDS) [55] . For CA, CCA, DCA, RDA, and DPCoA, the ordinationterm is based upon an evaluation of abundance values (in the case of DPCoA, the patristic distancesterm between OTUsterm on the phylogenetic tree is also used), but not an ecological distanceterm. For MDS and NMDS, the ordinateterm function requires a pre-calculated distance matrixterm (“dist” object) or the name of a supported ecological distanceterm method. For example, PCoA/MDS can be calculated on an unweighted UniFrac distance matrixterm [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 phyloseqterm package also includes a specialized function for displaying ordinationterm results in different ways, described in the following section.

Specialized Graphics One of the key features of the phyloseqterm package is a set of graphics functions custom-tailored for phylogenetic sequencingterm analysis, built using the ggplot2term package [57] . The ggplot2term 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 attributesterm [58] . In the phyloseqterm package, functions having names beginning with “plot_” require a phyloseqterm object as input data, and return a ggplot2term graphics object. These plot_ functions support optional mapping of color, size, and shape aesthetics to sample or OTUterm 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 Rterm graphics), and all features of these graphics can be further modified in Rterm via functions/options in the ggplot2term package.

The following list summarizes the key graphics-producing functions in phyloseqterm, which are also demonstrated in Figure 4 , and in phyloseqterm's online tutorials [39] . File S2 provides the complete Rterm code for creating Figures 4 and 5 . We have also included some additional examples of graphics created by plot_ordinationterm ( Figure 5 ). They emphasize different aspects of ordinationterm 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 ggplot2term graphic, in this case the addition of a two-dimensional density estimate to the “OTUs-only” plot (File S2).

Graphic functions of the phyloseqterm package. The phyloseq classterm is an experiment-level data storage classterm defined by the phyloseqterm package for representing phylogenetic sequencing dataterm·def. Most functions in the phyloseqterm package expect an instance of this class as their primary argument. See the phyloseqterm manual The Global Patterns [47] and Enterotypes [91] datasets are included with the phyloseqterm package. The Global Patterns data was preprocessed such that each sample was transformed to the same total read depth, and OTUsterm were trimmed that were not observed at least 3 times in 20% of samples or had a coefficient of variation ≤ 3.0 across all samples. For the plot_treeterm and plot_bar subplots, only the Bacteroidetes phylum is shown. Each subplot title indicates the plot function that produced it. Complete details for reproducing this figure are provided in File S2. All of these functions return a ggplot object that can be further customized/modified by tools in the ggplot2term package [57] . See additional descriptions of each function in the body text, and at the phyloseqterm homepage [39] .

plot_ordinationterm display methods included in phyloseqterm. Each panel uses a “Bacteroidetes-only” subset of the preprocessed “Global Patterns” dataset that was also used in Figure 4 . The coordinates are derived from an unconstrained correspondence analysisterm [62] . Different panels illustrate different displays of the ordinationterm results using the type argument to the plot_ordinationterm function. (Top Left) Example of a samples-only display, with the “SampleType” mapped to the color aesthetic, and a filled-polygon layer to emphasize plot regions where sample types co-occur. (Top Left Insert) A “scree” plot of the eigenvalues associated with each axis, which indicates the proportion of total variability represented in each axis. (Top Right) Biplot representation in which samples and OTUsterm ordinationterm results are overlaid. Clumps of OTUsterm appear to co-occur with different sample types, and some correlation with taxonomic phylum is also evident. (Middle) An OTUs-only plot that has been faceted (separated into panels) by class, with a two-dimensional density estimate overlain in blue. This view shows clearly a lack of association between the Sphingobacteria and Flavobacteria classes with fecal samples, which appear to be enriched in a subset of the Bacteroidia (relative to other OTUsterm in this Bacteroidetes-only dataset). Meanwhile, subsets of Bacteroidia appear to be enriched within multiple sample types. (Bottom) The “split” type for this graphic, in which both samples-only and OTUs-only plots are created, and shown side-by-side with one legend and shared vertical axis. Both the “biplot” and “split” options allow dual projections of both OTU- and sample-space.

plot_ordinationterm. This is the main function for plotting the results of an ordinationterm. It currently supports four different representations of the ordinationterm results: samples-only, OTUs-only, “biplot” (combined) representation, and “split”. A demonstration of these different options is provided in Figure 5 . As can be seen in these examples, the “biplot” and “split” options support dual projections of both OTU- and sample-space. Additional parameters easily map the respective sample variable or taxonomic rank to color, size, or shape aesthetics.

plot_heatmapterm. This is a special implementation of the ordination-organized heat map similar to the NeatMap package [59] . Briefly, the abundance matrix is represented as a grid of colored tiles, with the color of the tiles mapped to the (usually transformed) abundance value. The ordering of the OTUsterm and sample indices in this representation is critical for discriminating any patterns. Traditionally, hierarchical clusteringterm methods have been used for this organization; but, as Rajaram and Oono recently pointed out [59] , this has the potential to misrepresent the data when deeply-branching elements are placed next to one another arbitrarily. Instead, the samples (and optionally, OTUsterm) are reordered based on their radial coordinate angle in the first two axes of an ordinationterm. For the plot_heatmapterm function, any of the distancesterm/ordinationsterm supported by the distanceterm and ordinateterm functions can be used, with the default being non-metric multidimensional scalingterm. Any arbitrary color scale can be selected, as well as any choice of numerical transformation for scaling the mapping of color shades to abundance.

plot_networkterm. This function plots an igraph-class network [60] representing binary relationships between samples or OTUsterm. The network is calculated using the make_network function with phyloseqterm data as input and a desired ecological distanceterm and threshold value. Unlike ordinationterm, where most of the data structure is summarized by the relative position in two or more axes, the data is instead summarized by connections between samples (or OTUsterm) drawn with straight lines. Two samples are considered “connected” if the distanceterm between them is less than a user-defined threshold. The relative position of points is optimized for the visual display of network properties, but is otherwise arbitrary. Any of the ecological distancesterm supported by the distanceterm function can be selected, and this can be a powerful representation of major clusters among samples or OTUsterm, provided the value of the distance thresholdterm has been chosen carefully.

plot_treeterm. This function facilitates easy graphical rendering/investigation of the phylogenetic tree, with sample data overlaid. In some cases an annotated tree can be a powerful representation of an underlying evolutionary structure. The plot_treeterm function optionally places successive points next to the tips of the tree, indicating the samples in which each OTUterm was observed. These points can have their color, shape, and size aesthetics mapped to sample variables, revealing the correspondence of environmental variables on specific regions of the evolutionary tree. Standard ggplot2term customizations are supported, and this is, to our knowledge, the only function for ggplot2-based phylogenetic trees currently available in the CRAN/Bioconductorterm repositories. For phylogenetic sequencingterm of samples with large richness, some of the options in this function will be prohibitively slow to render or too dense to be interpretable, a drawback to summarizing phylogenetic sequencing dataterm·def using trees. One suggestion is to either agglomerate or subset the data such that there are not more than 200 or so OTUsterm (tree tips) on a given plot, sometimes less depending on the complexity of the additional annotations being mapped to the tree. In many modern datasets 200 OTUsterm (or less) will be insufficient to summarize the entire dataset, in which case one or more of the other plot methods is suggested.

plot_bar. Although sometimes very complicated, a well-organized bar plot can be an effective graphical means for direct quantitative comparison of abundance values, and we note that statisticians generally discourage the use of pie-charts [61] . The plot_bar function takes as input a phyloseqterm dataset and a collection of arbitrary expressions for grouping the data based upon taxonomic rank and sample variables. The returned graphic represents each abundance value as the height of a rectangular block that is outlined by a thin black line and filled with the corresponding color of the user-specified sample or taxonomic variable, grey by default. Each of these OTU abundanceterm·var rectangles corresponding to the same horizontal position (usually sample, or sample group) are stacked in order of abundance, such that the aggregate height of the stacked bar is also quantitatively informative.

plot_richness . This function creates plots of richness estimates of each sample in a phyloseqterm data object, allowing for horizontal grouping and color shading according to additional sample variables. Differences in richness (alpha diversity) between samples is often one of the first questions asked of phylogenetic sequencing dataterm·def.

Normalization and Standardization 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 abundanceterm·var tables, such a procedure is inappropriate as the disparities in column sumsterm·var can be 100-fold. Methods based on chi-squared distancesterm rather than variances deal with this by comparing weighted column profiles [62] , computed as relative abundances for each OTUterm within a column, with the overall column sumterm·var retained as a weighting factor. However, chi-square distancesterm 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 distancesterm can be a useful alternative, as it is based on the distanceterm between profiles, as long as the differences in actual column sumsterm·var are also accounted for in the final study. The other approach to the problem of disparities between column sumsterm·var has been to subsample the over-abundant columns down to the same number as the smaller ones. However this results in a loss of informationterm, rarely an optimal procedure in statistical contexts. This subsampling procedure is inspired by the popular idea of rarefactionterm 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 OTUsterm is beneficial, especially in the latter multiple testing adjustments. The phyloseqterm package enables easy filtering and rank transformations in the same vein as robust multi-array averaging (rma) [66] . We provide further details in (McMurdiepersName and HolmespersName, [67] ).

Confirmatory Analyses Although useful for exploring and summarizing microbiome data, many of the graphics and ordination methodsterm·var 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 countsterm are never normally distributed, so the preferred methods are nonparametricterm. Second, there is contiguous information available about the relationships between OTUsterm, 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 OTUsterm. The goal is to confirm that a certain property of these significant OTUsterm is overrepresented compared to the general population of OTUsterm, often called “the universe”. For instance in HolmespersName 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 nonparametricterm 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 parametersterm·var. 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 Rateterm (FWERterm). For instance, Westfall and Young's permutation-based minPterm procedure controls the FWERterm [71] and is implemented within the multtest package [72] . The phyloseqterm package interfaces with minPterm in multtest through a wrapper function, called mt. In the following example code we use the mt wrapper to control the FWERterm while simultaneously testing whether each OTUterm 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 ).

Results from the mt function on the “Enterotypes” dataset.

Results and Discussion As the complexity and sophistication of phylogenetic sequencingterm experiments continues to increase, it is clear that a “one-analysis fits all” approach is not sufficient. While it is often useful and convenient to have common analyses coupled within the application that decodes the sequences and clusters OTUsterm, we posit that a separate set of flexible open-source analytical tools is also needed that can be reproduced consistently by peers, and easily applied to new datasets and data sources. It should include a large library of statistical functions, and be independent of the choice of OTU-clustering method or sequencing technology. The phyloseqterm package helps satisfy this need by reducing the effort necessary to analyze OTU-clustered phylogenetic sequencing dataterm·def via the Rterm language and interactive computing environment.

Reproducible Researchterm and Sharing 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 researchterm 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 dataterm 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 abundanceterm·var 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 parametersterm·var 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 researchterm [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 Bioconductorterm/Rterm. The emphasis of preprocessing tools in phyloseqterm 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 phyloseqterm project is to help close the gap in reproducible researchterm that presently exists between pipeline results and the additional analyses required by investigators. This can be achieved when phyloseqterm is used (possibly with other Rterm packages) in conjunction with documentation tools such as Sweaveterm [32] , knitrterm [33] , iPython [83] Notebook invoking the rmagic extension, or “ R flavored markdownterm” (RFM) [84] . The Sweave-format approach is part of the reproducible researchterm 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 Sweaveterm (.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 phyloseqterm [86] .

Extending phyloseqterm It is important to note that the new phyloseq-classterm·def is a significant departure from the originally-proposed phyloseq-class structureterm [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 sequencingterm experiment as a single object. The updated phyloseq-classterm 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-classterm·def do not need to (re)perform common validity checks because these checks are consolidated as part of the phyloseq-constructor methodterm.

Analysis tools available in Rterm but not explicitly wrapped in phyloseqterm are nevertheless available to users and developers via accessorsterm and other data infrastructure tools. This leverages the fact that phyloseqterm data components are based on standard Rterm data classes and easily used in other package settings in Rterm. 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-classterm·def (See File S2 for code, and the phyloseqterm demo [86] ). Similarly, as an open-source package in an open language/framework (Rterm), phyloseqterm can be easily included at the relevant steps in pipelines, workbenches, and GUIs now under active development (E.g. ClovR [15] , MG-RAST [19] , QIIMEterm [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, phyloseqterm.

Conclusions The phyloseqterm project is a new open-source software tool for statistical analysis of phylogenetic sequencing dataterm·def within the Rterm programming language and environment. The tools in phyloseqterm make it easy to read the data output of several of the most common OTUterm clustering pipelines, and also represents this data in a unified, integrated form amenable to many modern analysis methods. With this integrated representation of the data it is easy to use supervised methods — such as canonical correspondence analysisterm, discriminant correspondence analysisterm, sparse linear discriminant analysis, etc. — to explain clinical or environmental response variables. We hope that this will provide a gateway for users to take their analyses towards more robust nonparametricterm alternatives to classical least squares methods, and allow them to interact graphically with their data more easily and efficiently. By leveraging existing Rterm infrastructure for reproducible researchterm, the phyloseqterm project also enables reproducible preprocessing, analysis, and publication-quality graphics production — such that it is easy to document, share, and modify analyses of phylogenetic sequencing dataterm·def. The phyloseqterm package is released on Bioconductorterm [34] and developed collaboratively on GitHub [39] .

Availability and Requirements Project name: phyloseqterm Project Stable Release: http://www.bioconductorterm.org/packages/release/bioc/html/phyloseqterm.html Project Home Page: http://joey711.github.com/phyloseqterm/ Project Issue Tracker: https://github.com/joey711/phyloseqterm/issues Project Demo Page: http://joey711.github.com/phyloseq-demo/ Operating System(s): Platform Independent Programming Language(s): Rterm Other Requirements: Rterm, Rterm packages (ade4, ape, Biostrings, foreach, ggplot2term, igraph0, multtest, picante, plyr, reshape, RJSONIO, scales, vegan) License: AGPL-3 Supporting Information Summary of comparison between phyloseqterm and currently available software. This PDF file contains a table summarizing a comparison of supported capabilities between phyloseqterm and QIIMEterm [11] , mothurterm [12] , and the pair of packages OTUbaseterm [35] and mcaGUI [88] . A “+” or “–” indicates that the capability is not directly supported, respectively. A symbol or word instead of “+” implies that the capability is supported, but with an extra caveat or detail, further defined below the table, if necessary. This is not a comprehensive summary of the capabilities of each packages, but rather the capabilities of relevance to this article. The abbreviations CA, DCA, RDA, and DPCoA stand for the ordination methodsterm·var correspondence analysisterm, detrended correspondence analysisterm, redundancy analysis, and double principal coordinates analysis, respectively. Note that in some cases the capabilities deemed “+” in this table are only supported for amplicon sequencingterm based data, sometimes from a specific sequencing platform and with the 16S rRNAterm gene as target. However, the phyloseqterm package is implemented at a stage in the analysis process that can be more generally applied to any phylogenetic sequencingterm, including non-standard amplicon targets, shotgunterm metagenome sequencing, etc.

(PDF) Click here for additional data file.

Source materials for reproducing this manuscript. This is a compressed .zip directory containing the main source file in Sweaveterm .Rnw format [32] , as well as the additional files necessary to completely recreate the original manuscript submitted to PLoS ONE. For the uninitiated, Sweaveterm is a Rterm/LaTeX2e interleaved hybrid language format [32] that allows advanced typesetting description to accompany Rterm code and its output (including graphics). Also included is the RFM source file that was used to create Figures 4 and 5 , and its accompanying HTML output that includes additional documentation details, links, and intermediate graphics. This latter file is “sourced” (re-run) by the Sweaveterm commands if any of the expected output files are missing. This supporting information zip file also includes Rterm code (at the end of the RFM/HTML files) that demonstrates how to use a phyloseqterm data object as an argument to other Rterm functions. In this particular example, the bioenv function from the vegan package [92] is demonstrated.

(ZIP) Click here for additional data file.