Statistical approaches for differential expression analysis in metatranscriptomics

2.1 Simulating paired metatranscriptomes and metagenomes

We simulated microbial communities based on the 100 most-prevalent species from the healthy human gut microbiome (Lloyd-Price et al., 2017) which we then used as templates for generating synthetic MTX and MGX sequencing count data with spiked phenotypic associations (including simulated DE). The details of this process are summarized here and expanded in Supplementary Note S1. To simulate microbial community compositions, we sampled species in proportion to their empirical prevalence and selected abundances from fitted log-normal distributions. We modeled species pangenomes as bags of 1K genes drawn from a larger pool of 2K gene families. A community strain carriage of a given species was randomly drawn from the species’ 1K genes (gene content thus varies across simulated strains of a species, modeling genomic plasticity). Each gene family was additionally associated with a reference (log-mean) expression level m drawn from N(0, 1) with per-sample expression values drawn from log N(m, 1).

Samples were randomly assigned to a binary phenotype. One or more phenomena (read depth, species abundance, gene presence/absence and gene expression level) were associated with this phenotype by shuffling their values to induce rank-correlations of a pre-specified strength. Gene-phenotype associations involving expression level are positive for DE in subsequent model evaluations. We simulated 11 datasets of N = 100 samples containing different combinations of phenomena spiked to associate with sample phenotype (Table 1 and Supplementary Table S1). For example, the ‘null-bug’ dataset contains species-phenotype associations but no positive DE relationships, while the ‘true-combo-dep-exp’ dataset contains a mix of positive DE relationships and a potentially confounding sequencing depth association.

To simulate MGX count data, reads were sampled up to a sample’s sequencing depth and assigned to sample genes in proportion to the abundance of the gene’s source species. MTX count data were simulated similarly, with reads assigned to transcripts (genes) in proportion to the abundance of their source species and per-sample gene expression value.

2.2 Analytical approaches for DE analysis in metatranscriptomic data

We evaluated a series of linear models for community DE analysis incorporating one of six gene-copy normalization methods and three upstream zero-filtration methods alongside other parameter choices. Models were implemented using the glm function of R’s ‘stat’ package. Unless otherwise stated, species, gene and transcript counts were initially total-sum-scaled (to adjust for sample read depth) and then log-transformed (for variance-stabilization). Zero values surviving filtration were additively smoothed by half the corresponding feature’s smallest non-zero measurement prior to log-transformation.

2.3 Analysis of differentially expressed pilins in IBD

We used MTX and MGX data from the IBD Multi-omics Database (IBDMDB) cohort (Lloyd-Price et al., 2019) to evaluate the real-world performance of a subset of community DE methods. Data were downloaded from https://www.ibdmdb.org/ in December 2020. We focused on 800 samples with paired MTX and MGX profiles from 109 longitudinally-sampled participants: 52 with Crohn’s disease (CD), 30 with ulcerative colitis (UC) and 27 non-IBD controls. Longitudinal CD and UC samples from this cohort had been previously defined as ‘dysbiotic’ if they deviated from the ‘cloud’ of healthy microbiome configurations represented among non-IBD MGX profiles (with ‘deviation’ defined as a median Bray–Curtis distance to non-IBD samples exceeding the 90th percentile of non-IBD comparisons) (Lloyd-Price et al., 2019). IBDMDB had been previously analyzed using MetaWIBELE v0.3.8 (https://huttenhower.sph.harvard.edu/metawibele/): a method that (i) assigns candidate functions and taxonomy to genes in a microbial gene catalog and (ii) prioritizes those functions on the basis of functional, ecological and phenotypic properties. Based on this information, we selected a subset of 113 proteins for DE analysis that were (i) enriched in MGX samples from dysbiotic time points, (ii) taxonomically annotated to Escherichia coli and (iii) functionally annotated as pilin-related proteins on the basis of Pfam domain assignments (release Pfam32) (Finn et al., 2016).

MGX and MTX abundance profiles containing pilin-related proteins were initially community total-sum-scaled and then filtered using the semi-strict method to balance removal of technical zeros and loss of dysbiotic samples. Non-filtered RNA and DNA values were then log-transformed and smoothed following the procedures outlined above for synthetic data. Each pilin feature was analyzed for DE under the ‘dysbiosis’ phenotype using three gene-copy normalization methods (M3, M4 and M6; Fig. 1B). We re-implemented these methods in the lmerTest R package (Kuznetsova et al., 2017) to accommodate subjects’ repeated measures as a random effect, ‘(1|subject).’ ‘Diagnosis’ (CD versus UC versus Non-IBD), ‘antibiotic’ use and consent ‘age’ were added as additional covariates relevant to the IBDMBD cohort. Nominal P-values for UC- and CD-specific pilin DE in dysbiosis were subjected to multiple hypothesis testing correction using the Benjamini–Hochberg method with a FDR threshold of 0.25.

3.1 Statistical models accounting for underlying variation in gene copy number facilitate DE discovery in MTX

We initially evaluated the performance of six approaches to MTX normalization (M1-6, defined in Fig. 1B) using ‘strict’ pre-filtering of MTX zeros (Section 2). For models M3-6 (i.e. those incorporating an estimate of gene copy number), strict per-filtering excludes samples from the analysis of a given feature if the feature’s per-sample RNA abundance or gene-copy estimate were zero. The impact of choosing a less-strict pre-filtering approach is explored in the next section.

All models (M1-6) displayed close to optimal specificity (FPR ∼ 0.05) on the ‘null’ dataset, which contained no spiked phenotype associations (Fig. 2). As expected, model M1 (which incorporates only total community-level MTX normalization) was prone to FP DE detection when RNA counts varied due to confounding changes in gene copy number. This was true for confounding variation in species abundance (the ‘null-bug’ dataset; FPR = 0.110) and gene presence/absence (the ‘null-enc’ dataset; FPR = 0.082). The inability to isolate signals of gene expression from gene-copy variation further limited model M1’s ability to detect even strong positive DE signals (TPR = 0). Surprisingly, the next-most sophisticated model (M2, incorporating within-taxon normalization), suffered from similarly high FP rates in the presence of confounding DNA-level species and gene variation (FPR = 0.110 and 0.085), possibly due to the amplification of noise within species at low absolute RNA abundance. M2 was somewhat more successful than M1 at detecting strong DE signals (TPR = 0.091).

In contrast, the four models incorporating estimates of gene copy number (M3-6) continued to achieve near-optimal specificity even in the presence of confounding signals (Fig. 2). Models M3 and M5, which use taxon-level RNA and DNA (respectively) to estimate copy number, benefit here from strict pre-filtering of RNA zeros (e.g. by avoiding the potential to confuse RNA zeros resulting from gene loss with DE in the ‘null-enc’ dataset). M3 and M5 displayed slightly elevated FPs when sequencing depth was confounded with sample phenotype (‘null-dep’ dataset; FPR = 0.057 for M3 and 0.056 for M5), as did model M6, which uses a per-feature DNA covariate to adjust for gene copy number (FPR = 0.054). Model M4, which is based on the RNA/DNA ratio, was robust to FPs from sequencing depth confounding (FPR = 0.0493) but more likely than M6 to produce a FP given extensive confounding with community structure (‘true-combo-bug-exp’ dataset, FPR = 0.0557 for M4 versus 0.0382 for M6).

Methods M3-6 were broadly similar with respect to sensitivity for positive DE relationships (Fig. 2). Indeed, the strength of DE (from weak to strong across the ‘true-exp-low,’ ‘true-exp-med’ and ‘true-exp’ datasets) was a much stronger determinant of sensitivity than choice of model among M3-6. That said, we did observe a trend for models M3 and M5 to be slightly more sensitive than models M4 and M6 (e.g. TPR of 0.218 and 0.217 versus 0.191 and 0.181 among medium-strength DE trends), possibly due to the greater influence of read-level noise on the per-feature gene-copy estimates used by models M4 and M6. In summary, when using strict pre-filtering of RNA zeros, within-taxon RNA scaling with a taxon-RNA covariate (model M3) provided a good balance of sensitivity and specificity for DE in MTX in the absence of paired MGX data (but requires knowledge of gene-taxon relationships). If paired MGX data are available, models M4 and M6 perform well and do not require knowledge of gene-taxon relationships.

3.2 Prefiltering zero values distinguishes low-expression genes from non-encoding or non-detection events

In MTX, transcript measurements (based on integer read counts) are commonly zero-inflated to an even greater extent than single-organism RNA-seq. This is due to the high ratio of MTX features (∼0.1–1 M transcripts) to sequencing depth (∼10–100 M reads) in most communities and the striking dynamic range among those features (6–8 orders of magnitude) (Lloyd-Price et al., 2017, 2019). Many MTX zeros are due to non-detection events (expression below the limit of detection of a given sequencing depth) or non-encoding events (gene loss within a community’s taxon lineage), instead of only absence/down-regulation of expression.

To assess the impact of MTX zeros on DE analysis, we tested three zero pre-filtering strategies across the six statistical models introduced above (Supplementary Fig. S1). The first strategy, ‘lenient filtering’, treats a feature as uninformative for DE analysis if it was always zero at the RNA level (all models) or always zero in gene-copy estimates (models M3-6), but otherwise considers all sample values. The second strategy, ‘semi-strict filtering’, excludes a sample from the analysis of a given feature under models M3-6 if the feature’s RNA count and gene-copy estimate are both 0, but includes the sample if one or both values are non-zero. The final strategy, ‘strict filtering’ (as applied in Fig. 2), excludes a sample if either the RNA or gene-copy estimate is zero, thus including only those samples where both values are non-zero (conservatively treating all MTX zeros as probable non-detection events).

Overall, pre-filtering zeros that cannot contribute to significant transcriptional differences improved the performance of statistical models for DE discovery (Fig. 3). This is particularly intuitive for ‘uninformative’ zeros, i.e. samples with no DNA or RNA for a given gene resulting from non-encoding or non-detection events, which can inflate both FPs and false negatives (FNs). Lenient pre-filtering resulted in the greatest expansion in FPs across models, especially as applied to datasets with confounding species abundance and gene presence/absence signals (Fig. 3A and Supplementary Fig. S2); specificity was improved under semi-strict filtering (Fig. 3B and Supplementary Fig. S3) and maximized under strict filtering (Fig. 3C). Sensitivity followed a similar increasing trend, with a notable jump for the models using taxon-level gene-copy estimates (M3 and M5) under strict versus semi-strict filtering. M3 and M5 with semi-strict filtering were additionally prone to FPs in the presence of confounding genomic plasticity (the ‘null-enc’ dataset), as their gene-copy estimates are not informative for RNA zeros resulting from low expression versus gene loss. Notably, the model adjusting RNA for DNA copy number as a covariate (M6) most consistently controlled FPs while maintaining statistical power across the three pre-filtering approaches. Conversely, the taxon-RNA covariate model (M3), i.e. the best-performing approach in the absence of MGX data, was only reliable under strict pre-filtering of zero values.

3.3 Accounting for underlying differences in metagenomic abundance enables DE detection in the absence of taxonomic annotation

An important component of many microbial communities is the pool of DNA or RNA features that cannot be confidently assigned to specific taxa via homology or de novo assembly and binning. While the analyses above presumed that all features could be assigned to specific taxa (a prerequisite for models M2, M3 and M5; Fig. 1B), this assumption will be violated for a fraction of genes in most communities and for potential majorities of genes in poorly characterized communities. To address this phenomenon in the context of community DE analysis, we performed additional evaluations of models M1, M4 and M6 on two synthetic MTX datasets lacking taxonomic provenance (Table 1 and Supplementary Table S1). There, positive DE and confounding signals were spiked consistently over gene orthogroups as they occurred across community species (modeling scenarios where all species contributing a certain function respond in the same way to a given condition). We then performed DE modeling on summed orthogroup abundances.

In the presence of confounding between gene presence/absence and phenotype (the ‘group-null-enc’ dataset), the naïve model (M1) was particularly vulnerable to FPs (FPR = 0.129; Fig. 4). Models M4 (RNA/DNA ratio) and M6 (feature-DNA covariate) achieved the target specificity of FPR ∼ 0.05 for these data under semi-strict or strict filtering of zero values, but M4 experienced elevated FPR with lenient filtering (FPR = 0.091). When applied to positive DE signals (the ‘group-true-exp’ dataset), models M4 and M6 were both highly sensitive and optimally specific independent of pre-filtering. In contrast to within-species DE (Fig. 2), model M1 was more effective at recovering positive community-level DE (TPR = 0.87), but continued to report FPs even in the absence of specific confounding signals (FPR = 0.096). These results stress that modeling adjusted RNA values to account for underlying variation in gene copy number consistently provides lower FPR and effectively aids in DE discovery in microbial communities containing unknown taxonomy.

3.4 Gene-copy normalization improves transcript presence/absence modeling

We separately evaluated the performance of six MTX normalization models (Fig. 1) in a logistic regression framework with an optional sequencing depth covariate. This approach directly models transcript presence/absence—rather than abundance, as above—as a function of sample covariates (including the phenotype of interest). We reasoned that this approach might be effective in detecting DE relationships at the level of transcript detection/non-detection, which proves to be more biologically relevant feature, while more elegantly handling genomic plasticity (gene gain/loss) and zero inflation.

Modeling trends incorporating logistic regression were broadly similar to those achieved using linear regression (Supplementary Fig. S4). Specifically, models incorporating gene copy-number normalization (M3-6) outperformed unadjusted within-community (M1) and within-taxon (M2) MTX normalization. Sensitivity was globally weaker after reducing continuous MTX abundance values to presence/absence calls; thus, we consider the logistic regression approach to be a supplement (but not alternative) to the previously-benchmarked linear regression models. The feature-DNA covariate model (M6) proved to be particularly specific in this domain, while models M3 and M5 (which use taxon-level estimates of gene copy number) were less robust to confounding gene presence/absence signals. Surprisingly, the inclusion of a sequencing depth covariate did little to improve specificity and tended to reduce power (sensitivity) across models, possibly due to the much lower degree to which it influences presence/absence than abundance for most features.

3.5 Adjusting RNA with DNA copy number as a covariate effectively identifies DE’ed inflammation-associated Enterobacteriaceae pilins in the human gut

Finally, we compared the behavior of selected well-performing models from the synthetic evaluations (M3, M4 and M6) on 800 stool samples with paired MTX and MGX data from the IBDMDB cohort (Lloyd-Price et al., 2019). These samples longitudinally span 109 subjects with three primary phenotypes (27 non-IBD controls, 52 Crohn’s disease or CD and 30 ulcerative colitis or UC). A subset of samples from CD/UC subjects were previously defined (Lloyd-Price et al., 2019) as ‘dysbiotic’ if they deviated from the cloud of healthy microbiome configurations covered by controls (Section 2). When evaluating genes from this dataset for potential community DE, we filtered samples using the semi-strict approach to maximize the number of dysbiotic time points available for analysis. Based on our findings above, M6 was expected to be more robust to this choice than M3 or M4 (Fig. 3).

We focused our analysis on 113 pilin-family proteins contributed by E.coli with potential DE in the dysbiosis state. E.coli is frequently implicated in IBD pathogenesis (Darfeuille-Michaud, 2002; Rolhion and Darfeuille-Michaud, 2007) where its pilin-family proteins play roles in bacterial adhesion, co-aggregation and biofilm formation (Hasegawa et al., 2009; Park et al., 2005). The pilin-family proteins selected here had been previously prioritized on the basis of their increased MGX abundance in dysbiotic samples (Section 2), and we hypothesized that they might exhibit additional increases in functional activity during inflammation (e.g. if they were playing roles as ‘drivers’ of dysbiosis versus ‘passengers’ in more-abundant E.coli genomes).

DE analysis of IBDMDB data using the feature-DNA covariate model (M6) highlighted a specific and dramatic transcriptional increase in E.coli genes encoding pilin-related proteins in dysbiotic samples, thus supporting their potential active role in IBD-associated inflammation (mixed-effects linear model, FDR q < 0.25; Fig. 5A). Since ground-truth DE status is not known for these features, we focused our evaluation on comparisons between results returned by models that had performed well on synthetic data (M3, M4 and M6). Among the 16 significantly up-regulated pilin-related proteins predicted by M6, the taxon-RNA covariate model (M3) tended to agree that these proteins were up-regulated, but only two achieved statistical significance (Fig. 5B). M3 additionally called 11 pilins as significantly up-regulated in dysbiosis and 3 as significantly down-regulated (Supplementary Fig. S5). While these discrepancies may indicate unique biological signals detected by M3, it is also possible that they arise from violations of model M3’s assumption of constant gene copy number, which is particularly inappropriate for highly strain-variant clades like E.coli.

Surprisingly, while the RNA/DNA ratio model (M4) tended to agree with M6 in its assessments of the statistical significance of pilin DE (Fig. 5C), it did so projecting trends with opposite (negative) sign for most pilins (Supplementary Fig. S5). This behavior is counter to biological intuition, in which inflammation-associated E.coli would be more likely to upregulate pilus structures. Where E.coli’s MGX abundance is also expanded in dysbiotic samples, RNA/DNA ratios for E.coli genes may be depressed to an excessive degree, leading to spurious evidence of down-regulation. Indeed, such interactions between species abundance and small RNA values (allowed under semi-strict filtering) were observed for M4 in the synthetic evaluations (Fig. 3B), where they provided an uncommon point of differentiation between models M4 and M6 (with M6 proving more robust to small values generally, consistent with its behavior in this application).