metamicrobiomeR: an R package for analysis of microbiome relative abundance data using zero-inflated beta GAMLSS and meta-analysis across studies using random effects models

Abstract

Background

The rapid growth of high-throughput sequencing-based microbiome profiling has yielded tremendous insights into human health and physiology. Data generated from high-throughput sequencing of 16S rRNA gene amplicons are often preprocessed into composition or relative abundance. However, reproducibility has been lacking due to the myriad of different experimental and computational approaches taken in these studies. Microbiome studies may report varying results on the same topic, therefore, meta-analyses examining different microbiome studies to provide consistent and robust results are important. So far, there is still a lack of implemented methods to properly examine differential relative abundances of microbial taxonomies and to perform meta-analysis examining the heterogeneity and overall effects across microbiome studies.

Results

We developed an R package ‘metamicrobiomeR’ that applies Generalized Additive Models for Location, Scale and Shape (GAMLSS) with a zero-inflated beta (BEZI) family (GAMLSS-BEZI) for analysis of microbiome relative abundance datasets. Both simulation studies and application to real microbiome data demonstrate that GAMLSS-BEZI well performs in testing differential relative abundances of microbial taxonomies. Importantly, the estimates from GAMLSS-BEZI are log (odds ratio) of relative abundances between comparison groups and thus are analogous between microbiome studies. As such, we also apply random effects meta-analysis models to pool estimates and their standard errors across microbiome studies. We demonstrate the meta-analysis examples and highlight the utility of our package on four studies comparing gut microbiomes between male and female infants in the first six months of life.

Conclusions

GAMLSS-BEZI allows proper examination of microbiome relative abundance data. Random effects meta-analysis models can be directly applied to pool comparable estimates and their standard errors to evaluate the overall effects and heterogeneity across microbiome studies. The examples and workflow using our ‘metamicrobiomeR’ package are reproducible and applicable for the analyses and meta-analyses of other microbiome studies.

Background

The rapid growth of high-throughput sequencing-based microbiome profiling has yielded tremendous insights into human health and physiology. However, interpretation of microbiome studies have been hampered by a lack of reproducibility in part due to the variety of different study designs, experimental approaches, and computational methods used [1, 2]. Microbiome studies may report varying results on the same topic. Therefore, meta-analyses examining different microbiome studies are critical to provide consistent robust results. Although many methods for microbiome differential abundance analysis have been proposed, methods for meta-analysis remain underdeveloped. Meta-analysis studies pooling individual sample data across studies for pooled analysis of all samples or processing of all samples together followed by analysis of each study separately have revealed some consistent microbial signatures in certain conditions such as inflammatory bowel disease (IBD) and obesity [3,4,5,6,7,8,9]. Software has been developed for the analysis and meta-analysis of microbiome data [10]. However, these studies do not explicitly model microbiome relative abundance data using an appropriate statistical method and do not examine between-group comparison overall pooled effects in the meta-analysis.

Data generated from high-throughput sequencing of 16S rRNA gene amplicons are often preprocessed into relative abundance. Microbiome relative abundances are compositional data which range from zero to one and are generally zero-inflated. To test for differences in relative abundance of microbial taxonomies between groups, methods such as bootstrapped non-parametric t-tests or Wilcoxon tests (not suitable for longitudinal data and covariate adjustment) [11,12,13] and linear or linear mixed effect models (LM) [14, 15] (suitable for longitudinal data and covariate adjustment) have been widely used. However, these methods do not address the actual distribution of the microbial taxonomy relative abundance data, which resemble a zero-inflated beta distribution. Transformations (e.g. arcsin square root) of relative abundance data to make it resemble continuous data to use in LM has been proposed by Morgan et al. (implemented in MaAsLin software) [16] and has been widely used to test for differential relative abundances [17,18,19,20]. However, this adjustment does not address the inflation of zero values in microbiome relative abundance data.

Various methods for the analysis of differential abundance based have been proposed. For example, the zero-inflated Gaussian distribution mixture model regards zero values as under-sampling and account for it by posterior probability estimates and fit counts after accounting for under-sampling by a log-normal distribution [21]. The Ratio Approach for Identifying Differential Abundance (RAIDA) method uses the ratio between the counts of features in each sample to address possible problems associated with counts on different scales within and between conditions and accounts for ratios with zeros using a modified zero-inflated lognormal (ZIL) model treating the zeros as under-sampling [22]. Other methods adapted from the RNA-seq field that account for zero inflation and utilize Poisson or negative binomial models have shown some promise in differential abundance testing of microbiome datasets [23, 24]. These aforementioned methods treat the dispersion as a nuisance parameter and do not allow the dispersion to depend on covariates. Recently, Chen et al. proposed an omnibus test based on a zero-inflated negative model (ZINB) that allows differential analysis not only for feature abundance but also prevalence and dispersion [25]. However, the downside of these count-based methods is the increased complexity due to modeling the counts.

Here, we developed an R package ‘metamicrobiomeR’ that applies Generalized Additive Models for Location, Scale and Shape (GAMLSS) with a zero-inflated beta (BEZI) family (GAMLSS-BEZI) for the analysis of microbial taxonomy relative abundance data. GAMLSS is a general framework for fitting regression type models in which the response variable can be any distribution [26]. With BEZI family, this model allows direct and proper examination of microbiome relative abundance data, which resemble a zero-inflated beta distribution. In principle, this model is similar to the two-part mixed effect model proposed by Chen et al. [27] in that the presence/absence of the taxon in the samples is modeled with a logistic component and the non-zero abundance of the taxon is modeled with a Beta component. Both logistic and beta components allow covariate adjustment and address longitudinal correlations with subject-specific random effects. The GAMLSS-BEZI is based on the broadly applicable established GAMLSS framework that can be flexibly implemented and applied to different types of data and study designs (e.g. cross-sectional and longitudinal). This is especially useful for later meta-analysis across different studies. The performance of GAMLSS-BEZI was evaluated using simulation studies and real microbiome data. Importantly, the estimates (regression coefficients) from GAMLSS-BEZI are log (odds ratio) of being in the case group (as compared to be in the control group) with changes in relative abundance of a specific bacterial taxon and thus are analogous across microbiome studies and can be directly combined using standard meta-analysis approaches. As such, we apply random effects meta-analysis models to pool the estimates and standard errors as part of the ‘metamicrobiomeR’ package. This approach allows examination of study-specific effects, heterogeneity between studies, and the overall pooled effects across studies. Finally, we provide examples and sample workflows for both components of the ‘metamicrobiomeR’ package. Specifically, we use GAMLSS-BEZI to compare relative abundances of the gut microbial taxonomies of male versus female infants’ ≤6 months of age while adjusting for feeding status and infant age at time of sample collection and demonstrate the application of the random effects meta-analysis component on four studies of the infant gut microbiome.

Implementation

GAMLSS-BEZI for the analysis of bacterial taxa relative abundance and bacterial predicted functional pathway relative abundance data

Relative abundances of bacterial taxa at various taxonomic levels (from phylum to genus or species) are obtained via the “summarize_taxa.py” script in QIIME1 [13]. Bacterial functional pathway abundances (e.g. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway level 1 to 3) are obtained from metagenome prediction analysis using PICRUSt [28]. In the taxa.compare function, all bacterial taxa or pathway data are first filtered to retain features with mean relative abundance ≥ relative abundance threshold (e.g. ≥0.005%) and with prevalence ≥ prevalence threshold (e.g. present in ≥5% of the total number of samples). This pre-filtering step has been shown to improve performance of various differential abundance detection strategies [29]. A filtered data matrix is then modeled by GAMLSS-BEZI and (μ) logit link and other default options using the R package ‘gamlss’ version 5.0–5 [26]. For longitudinal data, subject-specific random effects can be added to the model. We only include subject random intercepts as in practice this is often sufficient to address the longitudinal correlations [30]. However, it is possible to extend the model to include random slopes depending on the specific research content. For performance evaluation, LM and LM with arcsin squareroot transformation (LMAS) were also implemented in the function taxa.compare. In addition, we also implemented different approaches to deal with compositional effects including Centered Log Ratio (CLR) transformation [31] with various zero-replacement options [32] and Geometric Mean of Pairwise Ratios (GMPR) normalization [33]. Multiple testing adjustment can be done using different methods (False Discovery Rate (FDR) control by default). Below is an example call of the taxa.compare function:

taxa.compare (taxtab = taxtab, propmed.rel = “gamlss”, transform = “none”, comvar = “gender”, adjustvar = c(“age.sample”,“feeding”),longitudinal = “yes”, percent.filter = 0.05, relabund.filter = 0.00005, p.adjust.method = “fdr”).

For subsequent meta-analysis, the output from taxa.compare comprises matrices containing coefficients, standard errors, p-values and multiple testing adjusted p-values of all covariates in the models for each bacterial taxon or pathway.

Meta-analysis across studies using random effects models

The adjusted regression coefficient estimates from GAMLSS-BEZI are log (odds ratio) of being in the case group (as compared to be in the control group) with changes in relative abundances of a specific bacterial taxa or a pathway and thus are analogous across microbiome studies. Therefore, standard meta-analysis approaches can be directly applied. In the meta.taxa function, random effects meta-analysis models pooling adjusted estimates and standard errors with inverse variance weighting and the DerSimonian–Laird estimator for between-study variance are implemented to estimate the overall effects, corresponding 95% confidence intervals (CIs) and heterogeneity across studies. A fixed effect meta-analysis model is also implemented for comparison. Meta-analysis is performed only for taxa or pathways observed in ≥ a specified percentage threshold (e.g. 50%) of the total number of included studies. An example call to meta.taxa using the output data matrices combined from multiple calls to the taxa.compare function is shown below:

meta.taxa (taxcomdat = combined.taxa.compare.output, summary.measure = “RR”, pool.var = “id”, studylab = “study”, backtransform = FALSE, percent.meta = 0.5, p.adjust.method = “fdr”).

The output from meta.taxa consists of pooled estimates, standard errors, 95% CI, pooled p-values and multiple testing adjusted pooled p-values of all covariates for each bacterial taxon or pathway. The metatab.show function displays the meta-analysis outputs from meta.taxa as table, heatmap, forest plot or combined dataset to be used by the meta.niceplot function to generate nicer looking integrated heatmap-forest plot.

All implemented functions in the ‘metamicrobiomeR’ package are summarized and illustrated in Additional file 1.

Results and discussion

Performance of GAMLSS-BEZI: simulation studies

Simulation studies were performed to evaluate type I error and power of GAMLSS-BEZI for testing differential relative abundances of microbial taxonomies as compared to linear/linear mixed models with arcsin squareroot transformation (LMAS) (implemented in MaAsLin software [16]). LMAS was chosen for comparison with GAMLSS-BEZI because it is a commonly used approach for microbiome differential relative abundance testing and similarly to GAMLSS-BEZI, it allows covariate adjustment and can be used for longitudinal or non-longitudinal data. Simulations of zero-inflated beta distribution of microbiome relative abundance data were based on the R package “gamlss.dist” version 5.0–3.

In brief, beta distribution (denoted as Beta(μ, ϕ)) has a density function:

$$ f\left(y;\mu, \phi \right)=\frac{\ \boldsymbol{\Gamma} \left(\phi \right)}{\boldsymbol{\Gamma} \left(\mu \phi \right)\boldsymbol{\Gamma} \left(\left(\mathbf{1}-\mu \right)\phi \right)}{y}^{\mu \phi -1}{\left(1-y\right)}^{\left(1-\mu \right)\phi -1},y\in \left(0,1\right) $$

(1)

where 0 ≤ μ ≤ 1, ϕ > 0 and Γ (.) is the gamma function. If y~Beta(μ, ϕ), then E(y) = μ and Var(y) = μ(1 − μ)/(ϕ + 1), in which the variance of the dependent variable is defined as a function of the distribution mean μ and the precision parameter ϕ [34].

Zero-inflated beta distribution is a mixture of beta distribution and a degenerate distribution in a known value c = 0. A parameter α is added to the beta distribution to account for the probability of observations at zero producing a mixture density [34]:

$$ f\left(y;\alpha, \mu, \phi \right)=\left\{\begin{array}{c}\alpha, \kern6.5em if\ y=0,\\ {}\left(1-\alpha \right)f\left(y;\mu, \phi \right),\kern0.75em if\ y\in \left(0,1\right),\end{array}\right. $$

(2)

Type I error

We considered three sample sizes mimicking case-control microbiome studies with small (number of controls [n₁] = number of cases [n₂] = 10), medium (n₁ = n₂ = 100) and large (n₁ = n₂ = 500) scales. For each sample size, relative abundances of a bacterial species were simulated with the same parameters of a zero-inflated beta distribution for case and control groups (μ₁ = μ₂ = 0.5, α₁ = α₂ = 0.5, ϕ₁= ϕ₂ = 5). The simulation was repeated 1000 times. Type I error was calculated for three different alpha levels of 0.01, 0.05 and 0.1. Type I error of GAMLSS-BEZI or LMAS was defined as the proportion of simulations with p-values of GAMLSS-BEZI or LMAS less than the corresponding alpha level over 1000 simulations for each sample size. We noted that Type I errors were well controlled in both GAMLSS-BEZI and LMAS (Table 1).

Table 1 Type I error of GAMLSS-BEZI and LMAS

Receiver operating characteristic (ROC) curve and power

We then evaluated the performance of GAMLSS-BEZI vs. LMAS for identifying bacterial species with differential relative abundance between cases and controls. Two types of simulations were performed. First, relative abundances of 800 bacterial species were simulated in which 400 species had no difference between control and case groups (the same parameters of zero-inflated beta distribution for control and case groups: μ₁ = μ₂ = Uniform [0.0005,0.3], α₁ = α₂ = Uniform [0.1,0.9], ϕ₁= ϕ₂ = 5) and 400 species with a true difference between control and case groups. Specifically, four settings for the 400 species with true differences between control and case groups were considered with 100 species for each setting:

• 1) μ₁ = Uniform [0.0005,0.3] vs. μ₂ = μ₁ + 0.1

• 2) μ₁ = Uniform [0.0005,0.3] vs. μ₂ = μ₁ + 0.2

• 3) μ₁ = Uniform [0.0005,0.3] vs. μ₂ = μ₁ + 0.3

• 4) μ₁ = Uniform [0.0005,0.3] vs. μ₂ = μ₁ + 0.4

Other parameters (α, ϕ) were set the same for control and case groups (α₁ = α₂ = Uniform [0.1,0.9], ϕ₁= ϕ₂ = 5). A sample size of n = 100 for both case and control groups was used.

Performance of GAMLSS-BEZI and LMAS was evaluated based on the receiver operating characteristic (ROC) curve for identifying species with differential abundance between case and control groups. The analysis for the ROC curves and area under the curve (AUC) was done using the R package ‘pROC’ version 1.10.0. Under these settings, GAMLSS-BEZI (AUC = 95.6, 95% CI = [94.2, 97.1%]) significantly outperformed LMAS (AUC = 92.9, 95% CI = [91.1, 94.7%]) (DeLong’s test p-value < 2.2e-16) (Fig. 1a).

Fig. 1

ROC curve and power of GAMLSS-BEZI vs. LMAS. a. ROC curve of GAMLSS-BEZI and LMAS for identifying species with differential abundance between case and control groups. b. Power of GAMLSS-BEZI vs. LMAS for different effect sizes of differential relative abundances between case and control groups. GAMLSS-BEZI: Generalized Additive Models for Location, Scale and Shape (GAMLSS) with a zero inflated beta (BEZI) family; LMAS: linear model with arcsin squareroot transformation (implemented in the software MaAsLin); ROC curve: Receiver operating characteristic curve; AUC: area under the curve

We also performed simulations to evaluate power of GAMLSS-BEZI vs. LMAS for different effect sizes of differential relative abundances between case and control groups. Three settings for differential relative abundances (effect sizes) of one bacterial species were considered: 1) μ₁ = 0.5 vs. μ₂ = 0.4; 2) μ₁ = 0.5 vs. μ₂ = 0.3; and 3) μ₁ = 0.5 vs. μ₂ = 0.2. Other parameters were set the same for case and control groups (α₁ = α₂ = 0.5, ϕ₁= ϕ₂ = 5). A sample size of n = 100 for both case and control groups was used and the relative abundance of a bacterial species was simulated in each setting. The simulations were repeated 1000 times. Power of GAMLSS-BEZI or LMAS was calculated as the proportion of simulations with p-values of GAMLSS-BEZI or LMAS < 0.05 over the total number of 1000 simulations. Under these settings, power of GAMLSS-BEZI was better than power of LMAS (Fig. 1b).

Performance of GAMLSS- BEZI: application to real microbiome data

Type I error

We evaluated the type I error of GAMLSS-BEZI and LMAS using published data from a cohort study of 50 healthy Bangladeshi infants, which included longitudinal gut microbiome data from 996 stool samples collected monthly from birth to 2 years of life [14]. We used data from a subset of samples collected around birth as a cross-sectional dataset (50 samples) and data from all samples as a longitudinal dataset (996 samples). For each dataset, we randomly split the samples into two groups (case vs. control) and compared relative abundances of all bacterial taxa at all taxonomic levels (272 taxa from phylum to genus levels in total) between these two random groups using GAMLSS-BEZI and LMAS. The procedure was repeated 1000 times. Type I error was calculated for three different alpha levels of 0.01, 0.05 and 0.1. For each taxon, the type I error of GAMLSS-BEZI or LMAS was defined as the proportion of random splits with p-values of GAMLSS-BEZI or LMAS less than the corresponding alpha level over 1000 random splits. We noted that type I errors were well controlled in both GAMLSS-BEZI and LMAS (Table 2).

Table 2 Type I error of GAMLSS-BEZI and LMAS on real microbiome data

Computation time

The running time of GAMLSS-BEZI for testing all bacterial taxa at all taxonomic levels from phylum to genus (272 taxa in total) on a standard laptop were 6.4 s for the cross-sectional dataset (50 samples) and 12.4 s for the longitudinal dataset (996 samples), respectively. This indicates that the GAMLSS-BEZI algorithm is computationally efficient.

Detecting differential abundance

We evaluated the performance of GAMLSS-BEZI vs. LMAS in detecting differential relative abundances using published data from a cohort study of 50 healthy Bangladeshi infants described above [14]. This study included longitudinal monthly data regarding the infants’ breastfeeding practices (exclusive, non-exclusive), duration of exclusive breastfeeding, infant age (months) at solid food introduction, and occurrence of diarrhea around the time of stool sample collection. We compared the performance of GAMLSS-BEZI vs. LMAS in detecting differential relative abundances between various grouping variables in three examples below.

Example 1: Comparison of longitudinal monthly gut bacterial relative abundances at phylum level between non-exclusively breastfed (non-EBF) vs. exclusively breastfed (EBF) infants from birth to ≤ 6 months of age

Figure 2 (produced using the function taxa.mean.plot of our ‘metamicrobiomeR’ package; more details in Additional file 1) shows the longitudinal monthly average of relative abundance of bacterial phyla in non-EBF and EBF infants from birth to 6 months of age. A higher abundance of Proteobacteria, Firmicutes, and Bacteroidetes as well as a lower abundance of Actinobacteria are observed in non-EBF versus EBF infants. GAMLSS-BEZI is able to detect a significant difference in all four of these phyla whereas LMAS can only detect a significant difference in three phyla (Table 3).

Fig. 2

Relative abundances of bacterial phyla in non-exclusively breastfed vs. exclusively breastfed infants ≤6 months of age. Data from Bangladesh study

Table 3 Results of GAMLSS-BEZI and LMAS: real microbiome data example 1

Example 2: Comparison of longitudinal monthly gut bacterial relative abundances at phylum level between infants from 6 months to 2 years of age introduced to solid food after 5 months vs. before 5 months

Figure 3 shows the longitudinal monthly average of relative abundance of bacterial phyla in two groups of infants from 6 months to 2 years of age who were introduced to solid food after 5 months vs. those before 5 months of life. Lower relative abundances of Firmicutes, Bacteroidetes and higher relative abundance of Actinobacteria are observed in infants with solid food introduction after 5 months. GAMLSS-BEZI detects all three of these differences whereas LMEM can only detect a significant difference in one phylum (Table 4).

Fig. 3

Relative abundances of bacterial phyla in infants from 6 months to 2 years of age with solid food introduction after 5 months vs. before 5 months. Data from Bangladesh study

Table 4 Results of GAMLSS-BEZI and LMAS: real microbiome data example 2

Example 1 and 2 demonstrate the increased sensitivity of GAMLSS-BEZI in detecting bacterial taxa with observed differential relative abundances as compared to LMAS.

Example 3: Comparison of longitudinal monthly gut bacterial relative abundances at phylum level in infants from 6 months to 2 years of age with vs. without diarrhea stratified by duration of exclusive breastfeeding (EBF)

Figure 4 shows the average of relative abundance of bacterial phyla in groups of infants from 6 months to 2 years of age with vs. without diarrhea around the time of stool sample collection stratified by duration of EBF. In infants who received less than two months of EBF, a higher abundance of Firmicutes and a lower abundance of Actinobacteria is observed in the groups of infants with diarrhea vs. those without diarrhea (Fig. 4, upper panel). GAMLSS-BEZI detects a significant difference in both Firmicutes and Actinobacteria. In contrast, in infants who received more than two months of EBF, no difference in relative abundance of any bacterial phylum is observed between those with diarrhea vs. those without diarrhea (Fig. 4, lower panel) and GAMLSS-BEZI does not report any significant difference (Table 5). This example demonstrates that GAMLSS-BEZI detects differential abundances when there is observed difference and does not report difference when there is no observed difference.

Fig. 4

Relative abundance of bacterial phyla in infants from 6 months to 2 years of age with diarrhea vs. without diarrhea at the time of stool sample collection stratified by duration of exclusive breastfeeding (EBF). Data from Bangladesh study

Table 5 Results of GAMLSS-BEZI and LMAS: real microbiome data example 3

Illustration of meta-analysis examples with real microbiome data from four studies

We used gut microbiome data from four published studies to demonstrate the application of random effects models for meta-analysis across microbiome studies. These four studies include: 1) a cohort of healthy infants in Bangladesh [14] (the data of this study was also used in the three examples demonstrating the performance of GAMLSS-BEZI above); 2) a cross-sectional study of Haiti infants negative for HIV who were exposed or unexposed to maternal HIV [11]; 3) a cohort of healthy infants in the USA (California and Florida [CA_FL]) [12]; and 4) a small cohort of healthy infants in the USA (North Carolina [NC]) [35]. More details about the four studies included in the meta-analysis are described in Table 6. We illustrate the example of meta-analysis comparing relative abundances of gut bacterial taxa and bacterial predicted functional pathways between male vs. female infants ≤6 months of age adjusting for feeding status and infant age at the time of stool sample collection across these four studies (total number of stool samples = 610 [female = 339, male = 271]).

Table 6 Summary of four published microbiome studies included in meta-analysis

Relative abundances of gut bacterial taxa

Meta-analysis results are visually displayed using the functions metatab.show and meta.niceplot of our ‘metamicrobiomeR’ package (Additional file 1). The adjusted estimates (log (odds ratio) of one gender group for changes in relative abundance) from GAMLSS-BEZI for each bacterial taxon of each of the four studies and the pooled adjusted estimates across studies (meta-analysis) are displayed as a heatmap (Fig. 5 left panel). Different significant levels of p-values are denoted for each taxon of each study. The adjacent forest plot displays the pooled adjusted estimates and their 95% CI with different colors and shapes to reflect the magnitude of pooled p-values (Fig. 5 right panel).

Fig. 5

Meta-analysis for the difference in relative abundances of gut bacterial taxa between male vs. female infants ≤6 months of age. a: Phylum level: heatmap of log (odds ratio) (log (OR)) of relative abundances of all gut bacterial phyla between male vs. female infants for each study and forest plot of pooled estimates across all studies with 95% confidence intervals (95% CI). b: Genus level: heatmap of log (OR) of relative abundances of all gut bacterial genera between male vs. female infants for each study and forest plot of pooled estimates across all studies with 95% CI. All log (OR) estimates of each bacterial taxa from each study were from Generalized Additive Models for Location Scale and Shape (GAMLSS) with beta zero inflated family (BEZI) and were adjusted for feeding status and age of infants at sample collection. Pooled log (OR) estimates and 95% CI (forest plot) were from random effect meta-analysis models with inverse variance weighting and DerSimonian–Laird estimator for between-study variance based on the adjusted log (OR) estimates and corresponding standard errors of all included studies. Bacterial taxa with p-values for differential relative abundances < 0.05 are denoted with * and those with p-values < 0.0001 are denoted with **. Pooled log (OR) estimates with pooled p-values< 0.05 are in red and those with false discovery rate (FDR) adjusted pooled p-values < 0.1 are shown as triangles. Missing (unavailable) values are in white. USA: United States of America; CA: California; FL: Florida; NC: North Carolina

The running time for meta-analysis using both random effects and fixed effects models across four studies for all bacterial taxa (328 taxa available in at least 2 studies) from phylum to genus levels was 3.7 s on a standard laptop. This indicates that the meta-analysis algorithm is computationally efficient.

Across the four studies, there is a large heterogeneity in the difference (log (odds ratio)) of gut bacterial taxa relative abundances between male vs. female infants ≤6 months of age after adjusting for feeding status and age of infants at sample collection (Fig. 5, Additional file 1). For example, at the phylum level, relative abundance of Actinobacteria is significantly higher in male vs. female infants in two studies with small sample sizes (Haiti and North Carolina) while two other studies with larger sample size (Bangladesh and US (CA_FL) shows non-significant results in opposite directions. In addition, differential relative abundance of Proteobacteria is significant in two studies but in opposite directions (higher in male infants in the USA (CA_FL) study while lower in male infants in the Haiti study as compared to female infants). Moreover, at the genus level, each study shows significant differential relative abundances of different bacterial genera between male vs. female infants and the effects of many genera are in opposite directions between studies. Since the results are heterogeneous or opposite between studies and thus difficult to interpret, meta-analysis across studies is necessary to evaluate the overall consistent effects.

On the other hand, there are also some consistent effects across studies. For example, phylum Bacteroidetes is consistently decreased in male vs. female infants across four studies. However, the decrease is not significant in any study (Fig. 5a). Therefore, meta-analysis across studies is also important to evaluate if there is an overall significant effect.

Meta-analysis of the four studies shows no significant differential relative abundance of any bacterial phylum between male vs. female infants (Fig. 5a). At the genus level, meta-analyses show four genera with significant consistent differential relative abundances (pooled p-value < 0.05) between male vs. female infants. After adjusting for multiple testing, only genus Coprococcus remains significantly higher in male vs. female infants (FDR adjusted pooled p-value< 0.0001) (Fig. 5b).

Relative abundances of bacterial predicted functional (KEGG) pathways

Across the four studies, there is also a large heterogeneity in the difference (log (odds ratio)) of relative abundances of gut bacterial predicted functional KEGG pathways between male vs. female infants ≤6 months of age after adjusting for feeding status and age of infants at sample collection (Fig. 6). For example, at level 2 of KEGG pathway, the USA (CA_FL) study (with relatively large sample size) shows many pathways with significant differential relative abundances between male vs. female infants. The other three studies varyingly show significantly differential relative abundances in some of these pathways. However, the effects of almost all of these pathways in the USA (CA_FL) study are in opposite directions with the effects of these pathways in any of the other three studies. Therefore, it is difficult to interpret the results regarding male vs. female pathway differential relative abundances. As such, meta-analysis across studies is important to examine the overall consistent effects. Meta-analysis of four included studies shows only one KEGG pathway at level 2 with significant consistent differential relative abundance between male vs. female infants (pooled p-value < 0.05). However, after adjusting for multiple testing, no KEGG pathway (at both level 2 and level 3) remains significantly different between genders (Fig. 6, Additional file 1).

Fig. 6

Meta-analysis for the difference in relative abundances of gut microbial KEGG pathways between male vs. female infants ≤6 months of age. Heatmap of log (odds ratio) (log (OR)) of relative abundances of gut microbial KEGG pathways at level 2 between male vs. female infants for each study and forest plot of pooled estimates of all studies with 95% confidence intervals (95% CI). All log (OR) estimates of each pathway from each study were from Generalized Additive Models for Location Scale and Shape (GAMLSS) with beta zero inflated family (BEZI) and were adjusted for feeding status and age of infants at sample collection. Pooled log (OR) estimates and 95%CI (forest plot) were from random effect meta-analysis models with inverse variance weighting and DerSimonian–Laird estimator for between-study variance based on the adjusted log (OR) estimates and corresponding standard errors of all included studies. Pathways with p-values for differential relative abundances < 0.05 are denoted with * and those with p-values < 0.0001 are denoted with **. Pooled log (OR) estimates with pooled p-values< 0.05 are in red and those with false discovery rate (FDR) adjusted pooled p-values < 0.1 are shown as triangles. KEGG: Kyoto Encyclopedia of Genes and Genomes; USA: United States of America; CA: California; FL: Florida; NC: North Carolina

Difference in gut microbial composition between genders has been reported in adults [36, 37] and in some neonatal studies albeit with small sample sizes [38, 39]. However, the reported findings have largely varied between these studies. Our analyses also showed heterogeneous results among the four studies included. This highlights the importance of meta-analyses to evaluate overall consistent results across studies. Our meta-analyses of four studies showed virtually no difference in gut bacterial community and predicted functional pathways between male vs. female infants’ ≤6 months of age after adjusting for feeding status and infant age at time of sample collection as well as after adjusting for multiple testing. There was one exception: relative abundance of Coprococcus was significantly higher in male vs. female infants. Coprococcus has been implicated in many conditions including hypertension and autism [40, 41], and the detected difference in our study may provide some insights into the known sex differences in health outcomes.

In addition, random effects meta-analysis models can also be generally applied to other microbiome measures such as microbial alpha diversity and microbiome age. To make the estimates for these positive continuous microbiome measures comparable across studies, these measures should be standardized to have a mean of 0 and standard deviation of 1 before between-group-comparison within each study. Random effects meta-analysis models can then be applied to pool the “comparable” estimates and their standard errors across studies. Meta-analysis results of these measures can be displayed as standard meta-analysis forest plots (Additional file 1).

Conclusion

Our metamicrobiomeR package implemented GAMLSS-BEZI for analysis of microbiome relative abundance data and random effects meta-analysis models for meta-analysis across microbiome studies. The advantages of GAMLSS-BEZI are: 1) it directly address the distribution of microbiome relative abundance data which resemble a zero-inflated beta distribution; 2) it has better power to detect differential relative abundances between groups than the commonly used approach LMAS; 3) the estimates from GAMLSS-BEZI are log (odds ratio) of relative abundances of bacterial taxa between comparison groups and thus are directly analogous across studies. Random effects meta-analysis models can be directly applied to pool the adjusted estimates and their standard errors across studies. This approach allows examination of study-specific effects, heterogeneity between studies, and the overall pooled effects across microbiome studies. The examples and workflow using our “metamicrobiomeR” package are reproducible and applicable for the analysis and meta-analysis of other microbiome studies. The R package ‘metamicrobiomeR’ we developed will help researchers to readily conduct microbiome meta-analysis appropriately.

Availability and requirements

Project name: metamicrobiomeR.

Project home page: https://github.com/nhanhocu/metamicrobiomeR

Operating system(s): Platform independent.

Programming language: R.

Other requirements: R 3.4.2 or higher.

License: GNU GPL v. 2.

Any restrictions to use by non-academics: none.