MaAsLin2 validation

To identify model components appropriate for MaAsLin 2’s microbiome per-feature association testing and to objectively benchmark current association methods, we assessed realistic synthetic datasets generated by SparseDOSSA [22,23] (full details of individual association methods, as well as simulation parameters, are described in S1 and S2 Texts and are available online at https://github.com/biobakery/maaslin2_benchmark). Briefly, SparseDOSSA is a synthetic data generation routine that models biologically plausible synthetic data from diverse template microbiome profiles by considering (i) feature-feature, (ii) feature-metadata, and (iii) metadata-metadata correlations, superseding previous efforts by including multiple covariates and longitudinal designs (S1 Text). As compared to previous simulation schemes, SparseDOSSA allows multivariable spike-in both in the presence and absence of repeated measures, as well as arbitrary covariance structure in the metadata design matrix.

For this study, we carried out several spike-in experiments to induce and test controlled associations, as governed by configurable simulation parameters (S1 Fig). When used for this purpose, SparseDOSSA first generates null microbial community features containing no significant association patterns using a Bayesian hierarchical model independently of metadata features (Fig 1B and S1 Text). In addition to varying sample size and feature dimension, a broad range of metadata and experimental designs are then considered, including repeated measures and univariate and multivariate covariates (both continuous and binary) of varying dimension and effect size (S1A Fig). Specifically, in each instance, we varied sample sizes from small (10) to large (200) for a fixed feature size (up to 500), and within each sample size, the effect size parameter was again varied from modest (e.g., <2-fold differences) to strong (10-fold). In each simulation, 10% of features (and 20% of metadata for multivariable scenarios) were modified as an in-silico spike-in (S1 Text). Precision and recall measures (S1B Fig) were averaged over 100 simulation runs. All methods were corrected for multiple hypothesis testing using standard approaches for FDR control, declaring significant associations at a target of FDR 0.05. For a fair comparison, a basic, model-free filtering step to remove low-abundance features was performed before statistical modeling for all methods (S2 Text). Methods unable to process specific simulation configurations due to high computational overhead or slow convergence were omitted for those cases.

To compare the detection power of various methods in identifying true positive feature associations, we first comprehensively evaluated published metagenomic tools and representative methods from bulk RNA-seq literature within each simulation scenario. These methods were combined with several microbiome-appropriate normalization, transformation, and linkage models (S1C Fig and S2 Text). In particular, we considered six distinct categories of methods in our evaluations: (i) published methods specifically designed for microbial communities, such as metagenomeSeq [24], ANCOM [14,25], and ZIB [26,27], (ii) published bulk RNA-seq differential expression methods, such as DESeq2 [28], edgeR [29], and limma VOOM [30,31]; (iii) existing generalized linear model (GLM) count models, such as the negative binomial [17], (iv) methods based on linear models, such as limma [32] and “pure” linear models (LMs); (v) representative zero-adjusted methods from the microbiome and single-cell RNA-seq literature such as the Compound Poisson [18–20] and the ZINB [21,33]; and finally (vi) traditional, simplistic nonparametric methods, such as Spearman correlation and Wilcoxon tests. Of note, many of these methods can only compare two groups (i.e., a single binary metadatum), and not all are compatible with continuous and multivariate metadata, resulting in a distinct set of comparable methods for each experimental design.

Our first consideration in designing MaAsLin 2 for microbiome epidemiology was to ensure that both current statistical theory and practical issues were considered during the analysis of microbiome multi-omics data. To this end, we rigorously characterized finite-sample properties of various association methods focusing on three broadly defined aspects: (i) false discovery, (ii) detection power, and (iii) software implementation, with multiple performance indicators for each category (Figs 1C and S1B). Rather than focusing on a single evaluation metric like the Area Under the Curve (AUC) or the False Positive Rate (FPR), we ranked methods based on a combination of metrics (S1B Fig and S2 Text), many not considered in previous benchmarking. To summarize each evaluation criteria, a normalized continuous score ranging between 0 and 1 was assigned (S2 Text). Methods were then eliminated based on the presence of ‘red flags’ with respect to at least one evaluation criteria, i.e., extreme departure from the best possible value. Finally, metrics that are mainly descriptive rather than quantitative were also evaluated (e.g., the ability to handle complex metadata designs, zero-inflation, or repeated measures) to achieve a final consensus. For simplicity, we thus abbreviate any extreme departure from a metric’s best possible value as a ‘red flag’. This tiered strategy not only allowed us to select a set of "best" methods based on the fewest ‘red flags’ across all scenarios, but also to identify a method that is (i) sufficiently robust to false discovery control and detection power, (ii) scalable to large multi-omics datasets, and (iii) accommodating of most modern epidemiological designs and microbial data types.

Notably, previous benchmarking in this area has only focused on differential abundance testing without the simultaneous consideration of multiple covariates and repeated measures [7–9]. Additionally, with the exception of Hawinkel et al. [7], these benchmarking efforts lacked important considerations to the extent that they either (i) did not consider FDR as a metric of evaluation [9,34,35] or (ii) misreported false positive rate as FDR [8] (S2 Text). While most of these studies made a final recommendation based on the traditional AUC metric or a combination of sensitivity and specificity, we argue that without considering the FDR-controlling behavior of a method, the AUC values alone are too optimistic to draw any meaningful conclusions about its practical utility. In other words, particularly for biological follow-up, high precision among the most confident (lowest recall) predictions is essential. To this end, our large-scale benchmarking enables a progressive unification of traditional and practically important evaluation metrics by providing a comprehensive connected view of microbiome multivariable association methods, especially in the context of modern study designs, multiple covariates, and repeated measures.

Overall, our simulation study revealed that virtually all high-sensitivity methods with an overoptimistic median AUC, especially those targeted to microbial communities, exhibited a highly inflated average FDR (Fig 1C, full results in S1–S8 Data). A similar pattern was observed for other AUC-like measures such as F1 score and Matthew’s correlation coefficient (MCC). On the other end of the spectrum, compositionality-corrected methods such as ANCOM exhibited an extreme departure from ‘good’ performance with respect to several criteria including sensitivity and p-value calibration, as measured by both Conservative and Total Area [7] (S2 Text). Overall, these simulations reveal that while there is no single method that outperforms others in all scenarios, MaAsLin 2 was the only multivariable method tested that controlled FDR with the fewest ‘red flags’ across scenarios (Fig 1C).

This initial phase of our study thus expands the understanding of statistical association methods appropriate for microbial community data under varying study designs, and it especially highlights the inability of many common methods to control false discoveries. This is of critical importance to past and present microbiome association methods, as failure to control the FDR causes uncertainty in both scientific reproducibility and interpretability. Based on these evaluations, a linear model with TSS normalization and log transformation was adopted as the default model in MaAsLin 2, and the software provides several flexible options to apply a combination of other normalization, transformation, and statistical methods to customize specific analysis tasks. The implementation of MaAsLin 2, associated documentation, and example data sets are freely available both as an R/Bioconductor package and a command-line interface tool at https://huttenhower.sph.harvard.edu/maaslin2.

MaAsLin 2 controls false discovery rate while maintaining power in differential abundance analysis

Differential abundance testing for microbial community features (taxa, pathways, chemicals, etc.) is one of the most commonly used strategies to identify features that differ between sample categories such as cases and controls. Despite a large number of developments in the area, a lack of consensus on the most appropriate statistical method has been a major concern [11]. To model experimental designs of this type, we used synthetic count data with spiked-in features differentially abundant between two defined groups of samples. In particular, we multiplied the mean relative abundance of a randomly sampled fraction of 10% of the features with a given effect size (fold change) in one of the groups and renormalized the ensemble of relative abundances to a unit sum to create features differentially abundant between groups. We repeated this procedure for each unique combination of sample size (10, 20, 50, 100, 200), feature dimension (100, 200, 500), and fold change (1, 2, 5, and 10), each time summarizing performance over 100 simulation runs (S1 Text). Before model fitting, features with a low prevalence (<10%) were trimmed from the generated data sets.

As in our overall evaluation (Fig 1C), we observed marked differences between the FDR-controlling behavior of different methods in the simple case of single binary metadatum and non-longitudinal design, in some cases exceeding 75% (Fig 2). Among the methods with good, robust FDR control, only those based on linear models achieved moderate power, whereas, for methods such as DESeq2 and edgeR, the FDR control came at the cost of reduced power. Among other methods, practically all count and zero-inflated models, as well as newer methods based on log-ratios such as ANCOM, struggled to correctly control the FDR at the intended (nominal) level, and the best performance in this class of methods was obtained by metagenomeSeq2, Compound Poisson, and ZINB (as measured by the F1 score). Many of the remaining methods were too liberal, with metagenomeSeq and Negative Binomial standing out with many false positive findings. Overall, linear models (LMs) remained critically the only class of methods tested that has good control of FDR across study designs, and they all exhibited a boost in statistical power with increased sample size and association strength (S2 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. MaAsLin 2 controls false discovery rate while maintaining power in differential abundance analysis of microbial communities.

To assess models’ behaviors during differential abundance analysis, we simulated 100 independent datasets per parameter combination, each containing a single binary metadatum and a fixed number of true positive features (10% of features differentially abundant) for varying association strengths and sample sizes (S1A Fig). We then evaluated the ability of different microbiome association methods to recover these associations using a variety of performance metrics and summarized the results across runs. Both sensitivity and false discovery rates (FDR) are shown for the best-performing method from each class of models (as measured by average F1 score). Compared to zero-inflated and count-based approaches, MaAsLin 2’s linear model formulation consistently controlled false discovery rate at the intended nominal level while maintaining moderate sensitivity (full results in S1–S8 Data). Red line parallel to the x-axis is the target threshold for FDR in multiple testing. Methods are sorted by increasing order of average F1 score across all simulation parameters in this setting.

https://doi.org/10.1371/journal.pcbi.1009442.g002

We also evaluated the average FPR of these methods by recording the fraction of tested unassociated (negative) features that were deemed significant following significance testing. Nearly all methods controlled the FPR well below the imposed level (S3 Fig). Relatedly, we employed a previously proposed metric called “departure from uniformity” (i.e., departure from a uniform distribution of p-values under the null), which, complementary to FPR, quantifies the liberal or conservative area (S2 Text) between observed and theoretical quantiles of a uniform distribution [7]. As expected, methods with high average false discovery rates, including zero-inflated and count models, showed extreme departures from uniformity in the liberal direction, whereas conservative methods such as DESeq2 and edgeR showed the same in the opposite direction, suggesting extreme violation of uniformly distributed p-values under the null hypothesis (S4 Fig). While these results raise potential concerns about the FDR-controlling behaviors of most existing methods, LM-based approaches did not exhibit this trend. In general, tools based on linear models (such as limma) performed very similarly when calibrated with MaAsLin 2’s default model parameters, as expected, but not with their recommended default parameters (S2–S4 Figs). Additionally, their options for handling sparsity and compositionality were generally not appropriate for microbiome data. Amplicon, metagenomic taxonomic, and functional profiles each show distinct count and zero-inflation properties, for example, that are best handled by a multi-model system. In addition to the binary metadata design, we repeated the above simulation experiments for univariate continuous metadata as well, which led to similar conclusions (S5 Fig), further supporting MaAsLin 2’s default model’s performance across metadata types and experimental designs.

As a final evaluation, we assessed the impact of various normalization schemes on the associated statistical modeling, evaluating all combinations of normalizations appropriate for each applicable method (S1C Fig and S2 Text). Focusing on the best-performing linear models, we found that model-based normalization schemes such as relative log expression (RLE [36]) as well as data-driven normalization methods such as the trimmed mean of M-values (TMM [37]) and cumulative sum scaling (CSS [24]) led to good control of FDR, but they also led to a dramatic reduction in statistical power (S2–S5 Figs). In contrast, TSS showed the best balance of performance among all tested normalization procedures, leading to more powerful detection of differentially abundant features. These results have potential implications for other analyses in addition to differential abundance testing, as normalization is usually the first critical step before any analysis of microbiome data, and an inappropriate normalization method may severely impact post-analysis inference. In summary, our synthetic evaluation indicates that TSS normalization, although simplistic in nature, may be superior to other normalization schemes especially in the context of feature-wise differential abundance testing (and more generally for multivariable association testing, as described later), in addition to community-level comparisons as previously described [38].

MaAsLin 2 facilitates multivariable association discovery in population-scale epidemiological studies

Moving beyond univariate comparisons, we next assessed MaAsLin 2’s performance in multivariable association testing in comparison to other multivariable methods. Although widespread in microarray and gene expression literature, multivariable analysis methods have remained underdeveloped in microbial community studies. From an epidemiological point of view, coefficients from a covariate-adjusted regression model are arguably more interpretable than its individual, unadjusted counterparts. As a result, major conclusions from existing benchmarking studies geared towards univariate associations are not generalizable to this broader setting, where challenges such as zero-inflation and multiple testing are likely to be exacerbated, especially in relation to multiple rounds of independently conducted univariate analyses as commonly practiced.

To introduce multivariable associations into synthetically generated microbial feature profiles, we supplemented each “sample” with multiple covariates consisting of both binary and continuous metadata, either independent or correlated (S1A Fig and S1 Text). In each of these datasets, 10% randomly selected features were modified (“spiked”) to be associated with randomly chosen 20% metadata features with a given magnitude (effect size). After spiking in, samples were rescaled to their original (simulated) sequencing depth. As before, we repeated this procedure for each unique combination of sample size (10, 20, 50, 100, 200), feature dimension (100, 200, 500), and effect size (1, 2, 5, 10), each time summarizing performance over 100 simulation runs.

The results from this set of simulations revealed that MaAsLin 2’s default linear model had the highest sensitivity among the methods that controlled the FDR at the target level, which also remained consistent at larger sample sizes and stronger effect sizes (Fig 3). We also observed an improvement in performance when TSS normalization was employed (as compared to no normalization) but did not observe similar improvement for other normalization methods (S6 Fig). As before, zero-inflated and count models failed to control the FDR at the nominal level, in the sense that the actual FDR was always above the nominal threshold used for identifying significant features—a phenomenon that was surprisingly consistent regardless of the metadata covariance structure (S7 Fig). Taken together, these findings further confirm that MaAsLin 2’s default linear model is able to detect relevant associations across a broad range of metadata designs, facilitating population-level analyses of microbial communities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. MaAsLin 2 facilitates multivariable association discovery in large-scale human epidemiological and other microbial community studies.

Synthetic datasets containing five “metadata” with varying types of induced feature associations were analyzed using a variety of multivariable approaches (S1C Fig). As measured by power (recall) and false discovery rate (FDR), MaAsLin 2’s default linear model outperformed other methods in controlling FDR while maintaining power across true-positive fold-change values, regardless of the total number of features. As expected, MaAsLin 2 has better power for stronger effect sizes, eventually attaining the highest power among all FDR-controlling methods (full results in S1–S8 Data). Red line parallel to the x-axis is the nominal FDR. Values are averages over 100 iterations for each parameter combination. The x-axis (effect size) within each panel represents the linear effect size parameter; a higher effect size represents a stronger association. For visualization purposes, the best-performing methods from each class of models (as measured by average F1 score) are shown. Methods are sorted by increasing order of average F1 score across all simulation parameters in this setting.

https://doi.org/10.1371/journal.pcbi.1009442.g003

MaAsLin 2 enables targeted microbiome hypothesis testing in the presence of repeated measures

To further validate MaAsLin 2 for longitudinal (or other repeated measures) microbiome data, we modified our simulation scheme to introduce subject-specific random effects—a key feature of modern microbiome population studies [39]. To this end, we tested MaAsLin 2 and related methods on two types of study designs. The first comprised univariate binary metadata designed to be challenging by the inclusion of non-independence of the data across time points. Second, we also simulated more realistic datasets using multiple independent covariates specific to longitudinal microbiome studies. In both these regimes, realistic data were generated using SparseDOSSA each with five time points, as in previous studies [27], but after introducing within-subject correlations and between-subject random effects drawn from a multivariate normal distribution (S1 Text). It is to be noted that the set of evaluable models is greatly reduced from the previous set of cross-sectional association tests, as methods not capable of assessing repeated measures were discarded.

Using these longitudinal synthetic “microbial communities,” we compared the estimation and inference from MaAsLin 2 with those of the existing methods, which revealed that MaAsLin 2 had much lower false discovery rates than alternatives including ZIB (Figs 4 and S8–S11), a method specifically designed for microbiome longitudinal data. Both ZIB and MaAsLin 2’s linear mixed effects models are capable of identifying covariate-associated features by jointly modeling all time points. However, the computational overhead of ZIB is significantly higher than that of MaAsLin 2, which is prominent even for small datasets (S12 Fig). Notably, although not nearly as severe as count-based and zero-inflated models, MaAsLin 2 had a slightly inflated FDR in the univariate repeated measures scenario (Fig 4A) but not in the multivariable scenario (Fig 4B). Among other methods, methods based on generalized linear mixed models (GLMMs) such as Negative Binomial and Compound Poisson performed similarly to their non-longitudinal counterparts for both normalized and non-normalized counts (S8 and S9 Figs). This remained consistent for both univariate continuous metadata (S10 Fig) as well as multiple, correlated covariates (S11 Fig). Overall, these results suggest that MaAsLin 2’s linear mixed effects model consistently provides lower false discovery rates across metadata designs and can effectively aid in testing differential abundance and multivariable association of longitudinal microbial communities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. MaAsLin 2 enables targeted microbial feature testing in the presence of repeated measures.

Results on simulated data comprising SparseDOSSA-derived compositions with five repeated measures per sample. The FDR is close to the target 0.05 level for MaAsLin 2’s default linear model but not for zero-inflated and count models (full results in S1–S8 Data). As before, MaAsLin 2’s linear model is consistently better powered than both negative binomial and limma VOOM at comparable FDR values, which remains consistent for both univariate continuous metadata (A) and multivariable metadata designs (B). The red line parallel to the x-axis is the given threshold for FDR in multiple testing. Within each panel, methods are sorted by increasing order of average F1 score across all associated simulation parameters in each setting.

https://doi.org/10.1371/journal.pcbi.1009442.g004

Multi-omics associations from the Integrative Human Microbiome Project

We applied MaAsLin 2 to identify relevant microbial features associated with the inflammatory bowel diseases (IBD) using longitudinal multi-omics data from the Integrative Human Microbiome Project (iHMP or HMP2 [39]). The HMP2 Inflammatory Bowel Disease Multi-omics (IBDMDB) dataset included 132 individuals recruited in five US medical centers with Crohn’s disease (CD), ulcerative colitis (UC), and non-IBD controls, followed longitudinally for one year with up to 24 time points each (S3 Text).

Integrated multi-omics profiling of the resulting 1,785 stool samples generated a variety of data types including metagenome-based taxonomic profiles as well as metagenomic and metatranscriptomic functional profiles, producing one of the largest publicly available microbial multi-omics datasets. Metagenomes and metatranscriptomes were functionally profiled using HUMAnN 2 [40] to quantify MetaCyc pathways [41], and taxonomic profiles from metagenomes were determined using MetaPhlAn 2 [42] (S3 Text). For each of these data modalities (i.e., taxonomic profiles and DNA/RNA pathways), independent filtering was performed before downstream testing to reduce the effect of zero-inflation on the subsequent inference. In particular, features for which the variance across all samples was very low (below half the median of all feature-wise variances) or with >90% zeros were removed [39]. To further remove the effect of redundancy in pathway abundances (explainable by at most a single taxon), only features (both DNA and RNA) with low correlation with individual microbial abundances (Spearman correlation coefficient <0.5) were retained.

We first used the IBDMDB to perform an additional semi-synthetic evaluation of association methods’ performance in “real” data, specifically when attempting to associate randomized, null microbial taxonomic profiles to covariates (S3 Text). To this end, we permuted all samples 1,000 times to construct shuffled “negative control” datasets, each time assessing the number of significant associations (unadjusted p <0.05) for each applicable method. These were averaged across iterations to derive the expected number of null associations per method (which should remain near-zero under usual circumstances). In particular, we fit (i) a baseline model as a function of IBD diagnosis (a categorical variable with non-IBD as the reference group) while adjusting for enrollment age (as a continuous covariate) and antibiotic use (as a binary covariate), and (ii) a mixed effects model (with subject as random effects) with IBD dysbiosis state as an additional time-varying covariate in addition to the time-invariant covariates considered in the baseline model. Consistent with prior simulations, we found that several methods produced inflated empirical type I error rates (S13 Fig). This conclusion remained unchanged across varying significance thresholds, and as a result, we did not further apply these methods to the non-permuted data. Relevantly and importantly, linear models did not suffer from this problem, providing additional support for MaAsLin 2’s robustness to false positive findings.

To dissect dysbiotic changes in IBD at greater resolution, we applied MaAsLin 2 to each individual microbial feature type (i.e., species and DNA/RNA pathways) to test association with IBD phenotype while controlling for IBD dysbiosis state, diagnosis, age, and antibiotic use (Fig 5 and S3 Text). Nominal p-values for UC- and CD-specific effects were subjected to multiple hypothesis testing correction using the Benjamini-Hochberg method [43] with an FDR threshold of 0.25. MaAsLin 2 identified a comparable number of significant associations with those initially reported by the IBDMDB [39]. Among microbial species, MaAsLin 2’s default linear model identified 222 significant associations (S9 Data), among which 134 (60.4%) overlapped with the original study. MaAsLin 2 also reported many significant associations that were not discovered in the original study (S14 Fig). For instance, we observed a significant increase in Bacteroides ovatus in both UC and CD dysbiotic patients that was not previously captured, as well as detecting (with MaAsLin 2’s increased power) specific depleted Roseburia species (R. inulinivorans and R. hominis) not captured by the previous analysis. Notably, top hits from both MaAsLin 2 and the original study yielded nearly identical rankings across data types, which broadly manifested as a characteristic increase in facultative anaerobes at the expense of obligate anaerobes, in agreement with the previously observed depletion of butyrate producers such as Faecalibacterium prausnitzii in IBD (Fig 5A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Multi-omics associations from the Integrative Human Microbiome Project.

A) Top 10 significant associations (FDR < 0.25) detected by MaAsLin 2’s default linear model (full results in S9–S14 Data). All detected associations are adjusted for subjects and sites as random effects and for other fixed effects metadata including the subject’s age, diagnosis status (CD, UC, or non-IBD), disease activity (defined as median Bray-Curtis dissimilarity from a reference set of non-IBD samples), and antibiotic usage. B,C,D) Representative significant associations with dysbiosis state from each omics profile are shown: species (B), metagenomic (DNA) pathways (C), and metatranscriptomic (RNA) pathways (D). Values are log-transformed relative abundances with half the minimum relative abundance as pseudo count.

https://doi.org/10.1371/journal.pcbi.1009442.g005

As an additional validation, we next re-analyzed the HMP2 taxonomic and functional profiles using a zero-adjusted model (implemented in MaAsLin 2 as the Compound Poisson [18–20]). While this maintained type I error control in our shuffled data validation (as did the default linear model, S13 Fig), it was generally less desirable due to FDR inflation in simulations (Figs 2–4). In terms of the number of differentially abundant features detected, both the default linear model and the Compound Poisson model performed similarly, with a significant overlap between the top hits identified by each method (S15 Fig). Among other methods, ZIB and limma VOOM also maintained good Type I error control in these experiments (S13 Fig), although again both underperformed along other axes in our simulation studies. These results further strengthen the finding that a combination of controlled parametric simulations and ‘negative control’ experiments based on data shuffling are useful together in identifying methods for real-world applications, as the lack of either can lead to misleading (and irreproducible) conclusions across independent evaluations [7]. This also highlights the flexibility of MaAsLin 2’s multi-analysis framework, wherein researchers are well-served with multiple (i) normalization schemes, (ii) statistical models, (iii) multiplicity adjustments, (iv) fixed and random effects specifications, and (v) in-built visualization and pre-processing options, facilitating seamless application of methods across diverse experimental designs under a single estimation umbrella.

Finally, in addition to taxonomic associations, MaAsLin 2 also detected 399 and 58 significant functional associations for metagenomic (DNA) and metatranscriptomic (RNA) pathways, respectively (S10 and S11 Data), among which 358 (89.7%) and 39 (67.2%) overlapped with the original study. While the original analysis of these data included only community-wide functional profiles, we extended this by considering metagenomic and metatranscriptomic functional profiles at both whole-community and species-stratified levels to quantify overall dysbiotic functions while simultaneously assigning them to specific taxonomic contributors. In particular, this considers a per-feature DNA covariate model [44], in which per-feature normalized transcript abundance is treated as a dependent variable, regressed on per-feature normalized DNA abundances along with other regressors in the model (S3 Text). Surprisingly, bioinformatics and statistics for metatranscriptomics are not yet standardized, and our results indicate that subtle model variations can produce substantially different results, due to the interactions between two compositions: DNA and RNA relative abundances (S12 Data). This novel modeling strategy thus led to the discovery of several novel transcript associations relative to the original study (S13 and S14 Data).

In many of these pathways, functional perturbations were driven by shifts in their characteristic contributing taxa (Fig 5B). For example, the most significant DNA pathways enriched in CD patients were characteristic of facultative anaerobes such as Escherichia coli, which are broadly more abundant during inflammation. These included pathways such as synthesis of the enterobactin siderophore, lipid A, and sulfate reduction. A second set of enriched pathways was depleted due to the loss of microbes such as F. prausnitzii, a particularly prevalent organism that, when abundant, tended to contribute most of all enriched pathways it encodes in this cohort (e.g., synthesis of short-chain fatty acids and various amino acids).

With the increased sensitivity of this analysis for species-stratified pathways, the overwhelming majority of significant metagenomic differences were attributable solely to the most differential individual organisms, as expected (S13 and S14 Data). Essentially every pathway reliably detectable in E. coli was enriched during CD, UC, or both, and most F. prausnitzii pathways depleted, along with many pathways from other gut microbes common in “health” (Bacteroides vulgatus, B. ovatus, B. xylanisolvens, B. caccae, Parabacteroides spp., Eubacterium rectale, several Roseburia spp., and others). Interestingly, since both more potentially causal “driver” pathways, along with all other “passenger” pathways encoded by an affected microbe, are detected by this more sensitive stratified analysis, it can be in many ways more difficult to interpret than the non-stratified, community-wide, cross-taxon metagenomic responses to broad ecological conditions such as immune activity, gastrointestinal bleeding, or oxygen availability.

Conversely, differentially abundant microbe- and pathway-specific transcript levels highlighted a much more specific and dramatic shift toward oxidative metabolism, away from anaerobic fermentation, and towards Gram-negative (often E. coli) growth during inflammation (Fig 5C) [45]. Many of these processes were either more extreme during (e.g., gluconeogenesis) or unique to (e.g., glutathione utilization) active CD, as compared to UC. CD and UC responses were opposed in a small minority of cases (e.g., glutaryl-CoA degradation). When stratified among contributing taxa, these differences were almost universally attributable to a few key species, particularly an increase in E. coli activity during inflammation and decreases of F. prausnitzii transcript representation. Condition-specific transcriptional changes were occasionally contributed (or not) by “passenger” Bacteroides spp. (B. fragilis, B. xylanisolvens, B. dorei) instead. Note that these differences include pathways more likely to be “causal” in some sense, as significant transcriptional changes were generally a subset of those detected due to whole-taxon shifts in DNA content (including housekeeping pathways such as general amino acid or nucleotide biosynthesis). These findings further support the importance of disease-specific transcriptional microbial signatures in the inflamed gut relative to metagenomic profiles of functional potential, suggesting that a potential loss of species exhibiting altered expression profiles in disease may have more far-reaching consequences than suggested by their genomic abundances alone.