Phylogenetic factorization of compositional data yields lineage-level associations in microbiome datasets

Phylogenetically-structured compositional data

Microbiome datasets are “phylogenetically-structured compositional data,” compositions of parts linked together by a phylogeny for which only inferences on relative abundances can be drawn. The phylogeny is the scaffolding for the evolution of vertically-transmitted traits, and vertically-transmitted traits may underlie an organism’s functional ecology and response to perturbations or environmental gradients. Performing inference on the edges in a phylogeny driving variation in the data can be useful for identifying clades with putative traits causing related taxa to respond similarly to treatments, but such inferences must account for the compositional nature of the sequence-count data.

A standard analysis of microbiome datasets uses only the distal edges of the tree—the OTUs—and a few edges within the tree separating Linnaean taxonomic groups. However, a phylogeny of D taxa and no polytomies is composed of 2D − 3 edges, each connecting two disjoint sets of taxa in the tree with no guarantee that splits in Linnaean taxonomy corresponds to phylogenetic splits driving variation in our dataset. Thus, instead of analyzing just the tips and a series of Linnaean splits in the tree, a more robust analysis of phylogenetically-structured compositional data should analyze all of the edges in the tree. To do that, we draw on the isometric log-ratio transform from compositional data analysis, which has been used to search for a taxonomic signature of obesity in the human gut flora (Finucane et al., 2014) and incorporated into packages for downstream principal components analysis (Le Cao et al., 2016). However, to the best of our knowledge, the previous literature using the isometric log-ratio transform in microbiome datasets has used random or standard sequential binary partitions, and not explicitly incorporated the phylogeny as their sequential binary partition.

The isometric log-ratio transform of a rooted phylogeny

The isometric log-ratio (ILR) transform was developed as a way to transform compositional data from the simplex into real space where standard statistical tools can be applied (Egozcue et al., 2003; Egozcue & Pawlowsky-Glahn, 2005). A sequential binary partition is used to construct a new set of coordinates, and the phylogeny is a natural choice for the sequential binary partition in microbiome datasets. Instead of analyzing relative abundances, y_i, of D different OTUs, the ILR transform produces D − 1 coordinates, $x_i^{\ast }$ (called “balances”). Each balance corresponds to a single internal node of the tree and represents the averaged difference in relative abundance between the taxa in the two sister clades descending from that node (the difference being appropriately measured as a log-ratio due to the compositional nature of the data; see Supplemental Information 1 for more detailed description of the ILR transform). For an arbitrary node indicating the split of a group, R with r elements from the group, S with s elements, the ILR balance can be written as (1)${\displaystyle x_{\left\{R,S\right\}}^{\ast }=\sqrt{\frac{rs}{r+s}}log\left(\frac{g\left({\mathbf{y}}_R\right)}{g\left({\mathbf{y}}_S\right)}\right)}$ where g(y_R) is the geometric mean of all y_i for i ∈ R.

We refer to the ILR transform corresponding to a rooted phylogeny as the “rooted ILR”. The rooted ILR creates a set of ILR coordinates, $\left\{x_i^{\ast }\right\}$, where each coordinate corresponds to the “balance” between sister clades at each split in the phylogenetic tree. The balances in a rooted ILR transform in equation Eq. (1) can be intuited as the average difference between taxa in two groups, and splits in the tree which meaningfully differentiate taxa will be those splits in which the average difference between taxa in two groups changes predictably with an independent variable. Inferences on ILR coordinates, then, map to inferences on lineages in the phylogenetic tree.

The rooted ILR coordinates provide a natural way to analyze microbiota data as they measure the difference in the relative abundances of sister clades and may be useful in identifying effects contained within clades such as zero-sum competition of close relatives or the substitution of one relative for another across environments. However, if we desire to link the effect of a design variable or an external covariate (e.g., antibiotics vs. no antibiotic treatment) to clades within the phylogeny, the best comparison may not be between sister clades, but instead between all other clades, controlling for any other phylogenetic splits or factors we may know of (e.g., we may compare a lineage within gram-positives with all other gram-positives, once we’ve identified the gram-positive vs. gram-negative split as an important factor for antibiotic susceptibility). We refer to this unrooted approach as ‘phylofactorization.’

For the task of linking an external covariate to individual clades in the phylogeny, we examine three features of the rooted ILR that can be improved on by phylofactorization by considering a treatment that decreases the abundance of one and only one clade, B, whose closest relative is clade A. Regression on the rooted ILR coordinates may identify the balance $x_{\left\{A,B\right\}}^{\ast }$ corresponding to the most recent common ancestor of clades A and B as having that strongest response to the treatment, but regression on this coordinate will suggest that clade B decreases relative to A, leading to structured residuals in the original dataset due to an inability to account for the increase in clade B relative to the rest of the OTUs in the data (Fig. 1A). Second, all partitions between the affected clade and the root will be affected. If each balance is tested independently, the rooted ILR may identify many clades that are affected by antibiotics; the correlations between coordinates can yield a high false-positive rate if just one clade is affected (Fig. 1B). Finally, the ILR transformation does not work with polytomies common in real, unresolved phylogenies. Any polytomy will produce a split in the phylogeny between three or more taxa, and there is no general way to describe the balance of relative abundances of three or more parts using only one coordinate.

Nonetheless, the simplicity and theoretical foundations underlying the ILR, and the instant appeal of applying it to the sequential-binary partition of the phylogeny, motivate the rooted ILR as a simple tool for analysis of the phylogenetic structure in compositional data. For that reason, we use the rooted ILR as a baseline for comparison of our more complicated method of phylogenetic factorization.

Phylofactorization

The shortcomings of the rooted ILR can be remedied by modifying the ILR transform to apply not to the nodes or splits in a phylogeny, but to the edges in an unrooted phylogeny. While ILR coordinates of nodes allow a comparison of sister clades, ILR coordinates along edges allow comparison of taxa with putative traits that arose along the edge against all taxa without those putative traits. Traits arise along edges of the phylogeny and so, for annotation of phylogenies, effects in a clade are best mapped to a chain of edges in the phylogeny.

However, the ILR transform requires a sequential binary partition, and the edges don’t immediately provide a clear candidate for a sequential binary partition. In what we refer to as “phylofactorization,” one can iteratively construct a sequential binary partition from the unrooted phylogeny by using a greedy algorithm by sequentially choosing edges which maximize a researcher’s objective function. Phylofactorization consists of 3 steps (Fig. 2): (1) Consider the set of possible primary ILR basis elements corresponding to a partition along any edge in the tree (including the tips). (2) Choose the edge whose corresponding ILR basis element maximizes some objective function—such as the test-statistic from regression or the percent of variation explained in the original dataset—and the groups of taxa split by that edge form the first partition. (3) Repeat steps 1 and 2, constructing subsequent ILR basis elements corresponding to remaining edges in the phylogeny and made orthogonal to all previous partitions by limiting the comparisons to taxa within the groups of taxa un-split by previous partitions.

Explicitly, the first iteration of phylofactorization considers a set of candidate ILR coordinates, $\left\{x_e^{\ast }\right\}$ corresponding to the two groups of taxa split by each edge, e. Then, regression is performed on each of the ILR coordinates, $x_e^{\ast }\sim f\left(X\right)$ for an appropriate function, f and a set of independent variables, X. The edge, $e_1^{+}$, which maximizes the objective function is chosen as the first phylogenetic factor. In this paper, our objective function is the difference between the null deviance of the ILR coordinate and the deviance of the generalized linear model explaining that ILR coordinate as a function of the independent variables. We use this objective function as a measure of the amount of variance explained by regression on each edge because the total variance in a compositional dataset is constant and equal to the sum of the variances of all ILR coordinates corresponding to any sequential binary partition. Consequently, at each iteration there is a fixed amount of the total variance remaining in the dataset, and so at the candidate ILR coordinate which captures the greatest fraction of the total variance in the dataset is the one with the greatest amount of variance explained by the regression. After identifying $e_1^{+}$, we cut the tree in two sub-trees along the edge, $e_1^{+}$.

For the second iteration, another set of candidate ILR coordinates is constructed such that their underlying balancing elements are orthogonal to the first ILR coordinate. Orthogonality is ensured by constructing ILR coordinates contrasting the abundances of taxa along each edge, restricting the contrast to all taxa within the sub-tree in which the edge is found. A new edge, $e_2^{+}$, which maximizes the objective function is chosen as the second factor, the sub-tree containing this edge is cut along this edge to produce two sub-trees, and the process is repeated until a desired number of factors is reached or until a stopping criterion is met. More details on the algorithm, along with a discussion on objective functions, is contained in Supplemental Information 1.

While one could use other methods of amalgamating abundances along edges, the conceptual importance of using the ILR transform is twofold: the ILR transform has proven asymptotic normality properties for compositional data to allow the application of standard multivariate methods (Egozcue et al., 2003), and the ILR transform serves as a measure of contrast between two groups. The log-ratio used in phylofactor is an averaged ratio of abundances of taxa on two sides of an edge (see Supplemental Information 1 for more detail), thus phylofactorization searches the tree for the edge which has the most predictable difference between taxa on each side of the edge or, put differently, the edge which best differentiates taxa on each side. Thus, each edge that differentiates taxa and their responses to independent variables is considered a phylogenetic “factor” driving variation in the data.

The output of phylofactorization is a set of orthogonal, sequentially “less important” ILR basis elements, their predicted balances, and all other information obtained from regression. After the first iteration of phylofactorization, we are left with an ILR basis element corresponding to the edge which maximized our objective function and split the dataset into two disjoint sub-trees, or sets of OTUs that we henceforth refer to as “bins,” and we have an estimated ILR balancing element, ${\hat{x}}_1^{\ast }\left(X\right)$, where X is our set of independent variables. Subsequent factors will split the bins from previous steps, and after n iterations one has n factors that can be mapped to the phylogeny, n + 1 bins for binning taxa based on their phylogenetic factors, n estimates of ILR balancing elements, and an orthonormal ILR basis that can be used to project the data onto a lower dimensional space. The sequential splitting of bins in phylofactorization ensures sequentially independent inferences: having already identified group B as hyper-abundant relative to group A in the example illustrated in Fig. 2, downstream factors must analyze sub-compositions entirely within B and within A.

Computational tools

Phylofactorization was done using the R package “phylofactor” available at https://github.com/reptalex/phylofactor. The R package contains detailed help files that demo the use of the package, and the exact code used in analyses and visualization in this paper are available in Supplemental Information 1. The rooted ILR transform was performed as described in Egozcue & Pawlowsky-Glahn (2005) where the sequential binary partition was the rooted phylogeny.

Power analysis of rooted ILR and phylofactorization

To compare the ability of phylofactorization and the rooted ILR to identify clades of OTUs with shared associations with independent variables, we simulated random communities of D = 50 OTUs and p = 40 samples by simulating random absolute abundances, N_i,j, such that log(N_i,j) were i.i.d Gaussian random variables with mean μ = 8 and standard deviation σ = 0.5. The OTUs were connected by a random tree (the tree remained constant across all simulations), and then either 1 or 3 clades were randomly chosen to have associations with a binary “environment” independent variable with p = 20 samples for each of its two values to represent an equal sampling of microbial communities across two environments.

For simulations with one significant clade, the abundances of all the OTUs within that clade increased by a factor a in the second environment where a ∈ {1.5, 3, 6}. For simulations with three significant clades, the three clades were drawn at random and randomly assigned a fold-change from the set {π^b, 0.5^b, exp(−b)} in a randomly chosen environment where b ∈ {1, 2, 5}. For each fold-change, 500 replicates were run to compare the power of the rooted ILR and phylofactorization in correctly identifying the affected clades.

Regression of rooted ILR coordinates was performed and the coordinates were ranked by the difference between their null deviance and the model deviance. The ability of a rooted ILR coordinate to identify the correct 1 clade or 3 clades was measured by the percent of its top 1 or 3 ILR coordinates, respectively, which corresponded to the node on the tree from which the affected clade(s) originated. The ability of phylofactor to identify the correct 1 clade or 3 clades was measured by the percent of the factors that correctly split an affected clade from the rest (e.g., the percent of factors corresponding to edges along which a trait arose).

For the 3 clade simulations, we also compared the amount of variance explained by three factors in phylofactorization with the amount of variance explained by the top 3 ILR coordinates in the rooted ILR. The amount of variance explained was measured as the difference in the null deviance and the model deviance, summed across all three factors or the top 3 ILR coordinates.

KS-based stopping function for phylofactor

While a researcher can iterate through phylofactorization until a full basis of D − 1 ILR coordinates is constructed, the researcher may be interested in stopping the iteration before the full basis is constructed and focus their analysis and interpretations on a conservative subset of the true number of phylogenetic factors. We implemented a stopping function based on a Kolmogorov–Smirnov (KS) test of the distribution of P-values from analyses of variance of the regressions on candidate ILR coordinates. If there is no phylogenetic signal, we anticipate the true distribution of P-values to be uniform (albeit with some dependence among the P-values due to overlap in the OTUs used in the ILR coordinates). Thus, we tested the ability of phylofactor to correctly identify the number of clades if phylofactorization is stopped when a KS test of the P-values produces its own P-value P_KS > 0.05.

We simulated 300 replicate communities with M clades for each M ∈ {1, …, 10}. For simulations with M clades, D = 50 and p = 40 communities were simulated as above and fold changes, c, were drawn as log-normal random variables where log(c_k) were i.i.d Gaussian random variables with μ = 0 and σ = 3 for k = 1, …, M. The number of clades identified by phylofactor for a given true number of clades, K_M,r, was tallied for r = 1, …, 300. We calculate the mean $\overline{K_M}$ across all replicates and, for visualization purposes, interpolate the α = 0.025 and α = 0.975 quantiles by finding the best fit of a logistic function to the cumulative distribution of ${\left\{K_{M,r}\right\}}_{r=1}^{r=300}$ for each M.

Analysis of fecal/oral microbiome data

16S amplicon sequencing data from Caporaso et al. (2011) were downloaded from the MG-RAST database (http://metagenomics.anl.gov/) along with associated metadata. QIIME (Caporaso et al., 2012) was used to trim primers from these data, and to cluster OTUs with the Greengenes reference database (May 2013 version; http://greengenes.lbl.gov). Longer sequence lengths in the greengenes database (∼1,400 BP) compared to the original Illumina sequences (∼123 BP) allows more informative base pairs for phylogenetic tree construction. We used the phylogenetic tree that is included with the greengenes database for all analyses. The resulting OTU table was rarefied to 6,000 sequences per sample.

A total of 10 time points were randomly drawn from each of the male tongue, female tongue, male feces and female feces datasets, giving a total of n = 20 samples at each site. Taxa present in fewer than 30 of the 40 samples were discarded, and phylofactorization was done by adding pseudo-counts of 0.65 to all 0 entries in the dataset (Aitchison, 1986), converting counts in each sample to relative abundances, and then regressing the ILR coordinates against body site. The complete R script is available in the file “Data Analysis pipeline of the FT microbiome.”

Complete phylofactorization of this dataset was performed by stopping the algorithm when a KS-test on the uniformity of P-values from analyses of variance of regression on candidate ILR-coordinates yielded P_KS > 0.05. These results were compared with a standard, multiple hypothesis-testing analysis of CLR-transformed data. The summary of the taxonomic detail at the first three factors is provided in the results section, and a full list of the taxa factored at each step is available in the supplement and can be further explored using the R pipeline provided.

Analysis of soil microbiome data

The soil microbiome dataset from Ramirez et al. (2014) was included to illustrate the ability of phylofactor to work on bigger microbiome datasets with continuous independent variables and multiple regression. Details on sample collection, sequencing, meta-data measurements and OTU clustering are available in Ramirez et al. (2014). The phylogeny was constructed by aligning representative sequences using SINA (Pruesse, Peplies & Glöckner, 2012), trimming bases that represented gaps in ≥20% of sequences, and using fasttree (Price, Dehal & Arkin, 2010).

The complete dataset contained 123,851 OTUs and 580 samples. Data were filtered to include all OTUs with on average 2 or more sequences counted across all samples, shrinking the dataset to D = 3,379 OTUs. The data were further trimmed to include only those samples with available pH, C and N meta-data, reducing the sample size to n = 551.

Phylofactorization was done by adding pseudo-counts of 0.65 to all 0 entries in the dataset (Aitchison, 1986), converting counts in each sample to relative abundances, and performing multiple regression of pH, C and N on ILR coordinates. The first three factors are used for ordination-visualization. To determine the relative importance of each abiotic variable in driving phylogenetic patterns of microbial community composition, we used the lmg method from the R package ‘relaimpo’ (Grömping, 2006) which averages the sequential sums of squares over all orderings of regressors to obtain a measure of relative importance of each regressor in the multivariate model.

Power analysis and conservative stopping of phylofactorization

Phylofactorization remedies the structured residuals from the rooted ILR regression on data with fold-changes in abundances within clades. Phylofactorization also remedies the problem of high false-positive rates arising from the nested-dependence and correlated coordinates of the rooted ILR transform, as sequential inferences in phylofactorization are independent. Phylofactorization out-performs the rooted ILR in identifying the correct clades with a given fold-change in abundance (Figs. 3A and 3B), and can be paired with other algorithms assessing residual structure to stop factorization when there is no residual structure and thus accurately identify the number of affected clades (Fig. 3C). Finally, by focusing the inferences on edges instead of nodes in the phylogeny, this algorithm works on trees with polytomies and doesn’t require a forced resolution of polytomies to construct a sequential binary partition of the OTUs.

Soil microbiome

The soil microbiome dataset was much larger—3,379 OTUs and 580 samples—and a much smaller fraction of the variation could be explained by the dominant factors resulting from phylofactorization. Phylofactorization confirmed that the pH of the environment plays a dominant role in the microbial community composition, consistent with previous analyses based on Mantel tests (Ramirez et al., 2014). Dominance analysis of the generalized linear models associated with each factor determined pH to account for approximately 92.87%, 89.78%, and 92.94% of the explained variance in the first, second, and third factor, respectively. C and N were relatively unimportant, and the dominance of pH in the first three factors can be visualized by ordination-visualization plots of the ILR coordinates of the first three factors (Fig. 5A).

The first factor splits a group of 206 OTUs in two classes of Acidobacteria from all other bacteria: class Acidobacteriia and class DA052 are shown to decrease in relative abundance with increasing pH. The second factor split 31 OTUs in the order Actinomycetales (some from the family Thermomonosporaceae and the rest unclassified at the family level) from the remainder of all other bacteria, and these monophyletic Actinomycetales also decrease in relative abundance with increasing pH. The third factor identified another clade within the phylum Acidobacteria to decrease with pH: 115 bacteria from the classes Solibacteres and TM1.

The fourth factor identifies a monophyletic clade of 193 OTUs in the remainder of phylum Acidobacteria (i.e., those Acidobacteria not mentioned above in factors 1 and 3) as having relative abundances that increase with pH (dominance analysis: 94.79% of explained variance attributable to pH). Unlike the previous three factors above which were acidophiles, this monophyletic group of Acidobacteria consists of alkaliphiles, which includes the classes Acidobacteria-6, Chloracidobacteria, S053 and three OTUs unclassified at the class level.

The first four factors define 5 bins of OTUs that we refer to as “binned phylogenetic units” or BPUs: a monophyletic group of Acidobacteria (classes Chloracidobacteria, Acidobacteria-6, and S035), another monophyletic group of Acidobacteria (classes Solibacteres and TM1), a monophyletic group of several families of the order Actinomycetales, a monophyletic group of Acidobacteria (classes Acidobacteriia and DA0522), and a paraphyletic amalgamation of the remaining taxa. Binning the OTUs based on these BPUs tells a simplified story of how pH drives microbial community composition (Fig. 5B).

Overview

We have introduced a simple and generalizable exploratory data analysis algorithm, phylofactorization, to identify clades driving variation in microbiome datasets. Phylofactorization integrates both the compositional and phylogenetic structure of microbiome datasets and produces outputs that contain biological information: effects of independent variables on edges in the phylogeny, including the tips of the tree traditionally analyzed. The output of phylofactorization contains a sequence of “factors,” or splits in the tree identifying sub-groups of taxa which respond differently to treatment relative to one-another. The splits identified in phylofactorization need not be splits in the Linnaean taxonomy but can identify strong responses in clades of unclassified taxa. The researcher does not need to choose a taxonomic level at which to perform analysis—those taxonomic levels are output based on whichever clades maximize the objective function, and so researchers will be able to identify multiple taxonomic scales of importance.

Phylofactorization outputs an isometric log-ratio transform of the data with known asymptotic normality properties, coordinates that can be analyzed with standard multivariate methods (Pawlowsky-Glahn & Buccianti, 2011). The resulting coordinates correspond to particular edges between clearly identifiable nodes in the tree of life, allowing researchers to annotate a given phylogenetic tree with correlations between clades and various environmental meta-data, sample categories, or experimental treatments.

The phylogenetic inferences obtained by phylofactor would be difficult to obtain with other pre-existing methods. Three iterations of phylofactorizaiton on the oral/fecal dataset yielded the three major splits in the phylogeny, all of which are consistent with known distributions of taxa, none of which would be revealed directly from a taxonomic-based analysis. Algorithms such as phylosignal (Keck et al., 2016), which track P-values up the tree, yield inferences with nested dependence and identify clades with common significance, yet not necessarily clades with common signal or direction of trend—it is a common signal such as a shared habitat association, not a common significance such as the existence of a habitat association, which better indicates a putative trait driving shared responses in microbes. In the 142 factors above, phylofactor identified numerous clades with common significance yet different signals. Phylogenetic kernel-based methods (Lozupone & Knight, 2005; Purdom, 2011; Ning & Beiko, 2015) use the phylogeny to modify distances/dissimilarities between samples and can easily differentiate between oral and fecal body sites, but such nominally similar methods classify and differentiate sites, which is a conceptually different from classifying and predicting OTU abundances given their place in the phylogeny. Similarly, phylogenetic comparative methods (PCMs) (Felsenstein, 1985; Harvey & Pagel, 1991) such as phylogenetic generalized least squares (Grafen, 1989; Martins & Hansen, 1997) are nominally similar but conceptually distinct from phylofactor; while PCMs aim correct for the dependence of observations of evolutionarily related taxa, our goal is not to correct for such dependence but to infer precisely where on the tree such dependence appears to arise. While many tools exist for the nominally similar task of “using the phylogeny to analyze microbiome data”, the utility of phylofactorization lies in its ability to construct variables corresponding to edges in the phylogeny along which putative functional ecological traits may have arisen, constructed so as to avoid nested dependence and overlapping comparisons that frequent the analysis of hierarchically structured, and constructed in a manner appropriate for compositional data.

Future work

The generality of phylofactorization opens the door to future work employing phylofactorization with other objective functions. As we showed with the human oral/fecal microbiomes, phylofactorization is not restricted to basal clades, but includes the tips as possible clades of interest, but the objective function we used minimized residual variance in the whole community and thereby may prioritize deeply rooted edges or abundant taxa with weaker effects over individual OTUs with stronger effects. Other objective functions could be constructed to meet the needs of the researcher. If researchers are interested in identifying basal lineages, their objective function can weight edges based on distance from the tips. A researcher interested in fine-tuning the evolutionary assumptions in phylofactorization can define an objective function that increases explicitly with the length of the edge being considered to reflect an assumption that the probability of a trait arising increases with the amount of time elapsed.

Each edge identified in phylofactorization corresponds to two bins of taxa on each side of the edge, and consequently phylofactorization brings in two complementary perspectives for analyzing the data: factor-based analysis and bin-based analysis. Factor-based analysis looks at the each factor as an inference on an edge in the phylogeny, conditioned on the previous inferences already made, and indicating that taxa on one side of an edge respond differently to the independent variable compared to taxa on the other side of the edge. Bin-based analysis, on the other hand, looks at the set of clades resulting from a certain number of factors—what we call a “binned phylogenetic unit” (BPU). These bins will create a lower-dimensional, compositional dataset and can be freed from the underlying ILR coordinates for different analyses on these amalgamated clades. BPU-based analysis can inform sequence binning in future research aimed at controlling for previously-identified phylogenetic causes of variation, and combine the effects of multiple up-stream factors for predictions of OTU abundance. See Supplemental Information 1 for a more detailed discussion of factor-based and bin-based analyses.

One conceptual challenge with phylofactorization cross-validation and BPU-based analysis is the treatment of paraphyletic “remainder” clades which occur after factorization of edges within groups that split a monophyletic clade from a “remainder,” paraphyletic clade or BPU. These paraphyletic BPUs are informative as they are the necessary contrasts with the monophyletic group, conditioned on previously identified phylogenetic factors, and thus serve as a means of identifying the presence or absence of a trait in the monophyletic group. For instance, phylofactorization of vertebrates based on their relative abundance in the air should identify two monophyletic groups: birds and bats. The monophyletic groups are defined by the presence of a specialized adaptation—wings—and the paraphyletic remainder is defined by the absence of wings; any future treatment that differentially affects winged vs. non-winged organisms will differentially affect the abundances of these monophyletic groups relative to the paraphyletic remainder. Conversely, “legs” in the order Squamata disappeared along edges corresponding to legless lizards and snakes; in this case, the monophyletic clades are defined by the absence of a specialized adaptation, legs, and the paraphyletic remainder is defined by the presence of legs. Whether a trait can be flexibly defined as either the presence of absence of a specialized adaptation is an important theoretical consideration, but considering monophyletic BPUs and their complementary paraphyletic BPUs as equally likely to have the presence or absence of a trait and as the necessary comparison for cross-validation is an important empirical consideration. Phylofactorization pinpoints the contrasts to be made for microbial genomic and physiological studies aiming to identify the causes of differential abundances or responses to treatments across OTUs. Genomic/physiological investigations and cross-validation of phylofactorization must understand and grapple with the choice of appropriate comparisons of taxa, and that may require a fair contrast of a monophyletic clade to a paraphyletic one.

Phylofactorization will benefit from community discussion and further research overcoming general statistical challenges common to greedy algorithms and analysis of phylogenetically-structured compositional data. For instance, the log-ratio transform at the heart of phylofactorization requires researchers deal with zeros in compositional datasets. While there are many methods for dealing with zeros (Aitchison, 1986; Martín-Fernández, Barceló-Vidal & Pawlowsky-Glahn, 2003; Pawlowsky-Glahn & Buccianti, 2011), it’s unclear which method is most robust for downstream phylofactorization of sparse OTU tables. Second, phylofactorization as presented here only sequentially infers ILR elements and does not allow for simultaneous inference of ILR basis elements—the set of factors identified after n iterations may explain less variation combined than an alternative set of factors that did not maximize the explained variance at each iteration. This limitation may be overcome by running many replicates of a stochastic greedy algorithm and choosing that which maximizes the explained variance after n factors. Third, the researcher must choose an objective function which matches her question, and future research can map out which objective functions are appropriate for which questions in microbial ecology. Fourth, like any method performing inference based on phylogenetic structure, phylofactorization assumes an accurate phylogeny. Accurate statistical statements about a researcher’s confidence in phylofactors must incorporate the uncertainty in our constructed phylogeny. Fifth, phylogenetic-based methods may be sensitive to the binning and filtering methods for sequence-count data. Our filtering methods here were chosen to allow somewhat standard and simplified analysis of common OTUs present in many of the samples, but these filtering protocols might not be optimal for researchers with other research questions and objective functions. Future work can investigate the sensitivity of phylofactorization to myriad binning and filtering methods given the objectives and objective function of the researcher. Finally, future research can investigate the unique kinds of errors in phylofactorization: in addition to the multiple-hypothesis testing of edges, phylofactorization may propagate errors in the greedy algorithm, and, even when taxa are correctly factored into the appropriate functional bins, the presence of multiple factors in the same region of the tree can lead to uncertainty about the exact edge along which a putative trait arose (see Supplemental Information 1 for more discussion on the uncertainty of which edge to annotate).

Incorporating that phylogenetic structure into the analysis of microbiome datasets has been a major challenge (Martiny et al., 2015), and now phylofactorization provides a general framework for rigorous exploration of phylogenetically-structured compositional datasets. The soil dataset analyzed above, for instance, contains 3,379 OTUs and 580 samples, and phylofactorization of the clades affected by pH in the soil dataset yielded not just the three dominant factors used for ordination-visualization, but 2,430 factors in all, each with an intricate phylogenetic story. Many Acidobacteria are acidophiles, but some—Chloracidobacteria, Acidobacteria-6, S035, and some undescribed classes of bacteria factored here—appear to be alkaliphiles. By incorporating the phylogenetic structure of microbiome datasets, the big data of the modern sequence-count boom just got bigger, and future research will need to consider how to organize, analyze and visualize the large amounts of phylogenetic detail that can now be obtained from the analysis of microbiome datasets.