Comprehensive evaluation of methods for differential expression analysis of metatranscriptomics data

Abstract

Understanding the function of the human microbiome is important but the development of statistical methods specifically for the microbial gene expression (i.e. metatranscriptomics) is in its infancy. Many currently employed differential expression analysis methods have been designed for different data types and have not been evaluated in metatranscriptomics settings. To address this gap, we undertook a comprehensive evaluation and benchmarking of 10 differential analysis methods for metatranscriptomics data. We used a combination of real and simulated data to evaluate performance (i.e. type I error, false discovery rate and sensitivity) of the following methods: log-normal (LN), logistic-beta (LB), MAST, DESeq2, metagenomeSeq, ANCOM-BC, LEfSe, ALDEx2, Kruskal-Wallis and two-part Kruskal-Wallis. The simulation was informed by supragingival biofilm microbiome data from 300 preschool-age children enrolled in a study of childhood dental disease (early childhood caries, ECC), whereas validations were sought in two additional datasets from the ECC study and an inflammatory bowel disease study. The LB test showed the highest sensitivity in both small and large samples and reasonably controlled type I error. Contrarily, MAST was hampered by inflated type I error. Upon application of the LN and LB tests in the ECC study, we found that genes C8PHV7 and C8PEV7, harbored by the lactate-producing Campylobacter gracilis, had the strongest association with childhood dental disease. This comprehensive model evaluation offers practical guidance for selection of appropriate methods for rigorous analyses of differential expression in metatranscriptomics. Selection of an optimal method increases the possibility of detecting true signals while minimizing the chance of claiming false ones.

Keywords: benchmark; differential expression; early childhood caries; logistic-beta; metagenomics; metatranscriptomics.

Figure 1

Column A: Parameter estimates of baseline ZILN distributions obtained from 300 randomly selected genes in ZOE2.0 with the three-dimensional scatter plot on the top row and each of the subsequent rows representing estimates being within 0.03 from 0.9, 0.6 and 0.3. Column B: Disease effect estimates based on ZILN models obtained from the ZOE2.0 data in absolute values Column C: Batch effect estimates based on ZILN models obtained from the ZOE2.0 data in absolute values .

Figure 2

Goodness of fit (Kolmogorov–Smirnov, KS) test results for Beta, Log-normal and Gamma distributions (rows) with different scaling/transformation methods (columns). Histograms of the number P-values of the KS test, based on randomly select 300 genes in each of the three datasets (A) ZOE2.0, (B) ZOE-pilot, (C)IBD. The IBD data are available only in a compositional form, and thus we do not consider RPK in IBD. Lower rejection rate suggests better model fitting.

Figure 3

Type I error rates (under the Null D1) and FDR (Under the Alternative of mean, D2) for ZILN models. Columns correspond to sample sizes and evaluation criteria, rows are different tests, the axis represents baseline distributions and colors indicate batch effects. The dotted horizontal lines denote the significance level (5%). A failure in evaluation is marked as to be discerned from zero. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2, LFE = LEfSe, ALDEX = ALDEx2.Because p-values of LEfSe are not available, we do not obtain its FDR, also seen in Methods for Simulation I.

Figure 4

Sensitivity under alternative ZILN distributions for a small sample (). Columns and rows correspond to tests and alternative distributions, respectively, the axis represents baseline distributions and colors represent batch effects. The dotted horizontal lines denote the significance level (5%). A failure in evaluation is marked as to be discerned from zero. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2, LFE = LEfSe, ALDEX = ALDEx2.

Figure 5

Sensitivity under alternative ZILN distributions for a large sample (). Columns and rows correspond to tests and alternative distributions, respectively, the axis represents baseline distributions and colors represent batch effects. The dotted horizontal lines denote the nominal significance level (5%). A failure in evaluation is marked as to be discerned from zero. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2, LFE = LEfSe, ALDEX = ALDEx2.

Figure 6

Sensitivity curves of DE tests according to different cutoff values ranging from 0 to 0.2 and a few baseline and disease-effects scenarios of the ZILN model. No batch effects are simulated in these scenarios. The gray solid diagonal lines denote the nominal significance level.

Figure 7

Sensitivity of the DE tests according to different effect sizes for a subset of the baseline scenarios and . No batch effects are simulated. A failure in evaluation is marked as to be discerned from zero. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2, LFE = LEfSe, ALDEX = ALDEx2.

Figure 8

Performances of the analysis methods under semi-parametric simulations. 1000 (first three columns) or 10 000 genes (last three columns) were randomly selected for simulation. 10% of these genes, as specified in the panel heads, were given artificial disease effects. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2.

Figure 9

Type I error rates under the global Null condition generated by permuting the disease labels of samples in each of the three studies. 10 000 genes were tested. 100 permutations were generated. DS2 = DESeq2, DS2ZI = DESeq2-ZINBWaVE, ANCOM = ANCOM-BC2.

Figure 10

Application to the ZOE2.0 data analysis results. A. Histogram of the P-values of the log-normal models. B. Histogram of the joint P-values of the logistic Beta models (x-axis is for the Beta part and y-axis is for the logistic part). C. Histogram of the global P-values of the logistic Beta models (Wald test statistics). D. Scatter plot of the coefficients (disease effects on nonzero proportions and nonzero means) of the LB models, with the circled dots representing the most significant genes—Wald test statistic for the three types of Wald tests. Blue is for the global test, red is for the logistic part and green is for the Beta part. NA on the y-axis indicates that the logistic part was not estimated.

Figure 11

Venn diagram of DE genes in LN and the two parts of LB, at a P-value cutoff of in ZOE2.0 metatranscriptome data.

Figure 12

Application to the IBD data analysis results. A. Histogram of the P-values of the log-normal models. B. Histogram of the joint P-values of the logistic Beta models (logistic and Beta parts). C. Histogram of the global P-values of the logistic Beta models (Wald test statistics). D. Scatter plot of the coefficients of the LB models, with the circled dots representing the most significant genes—Wald test statistic . The NA results around indicate that the logistic part was not estimable.

Figure 13

Venn diagram of genes with P-values less than for each evaluated model in the IBD data.