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

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 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 groups and thus are comparable 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 workflow 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 heterogeneity and overall effects 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 robust results. 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 mirobiome relative abundance data using an appropriate statistical method and do not examined between-group-comparison overall pooled effects in the meta-analysis. For example, a recently published meta-analysis study examined differential relative abundances of bacterial taxonomies of each study separately using Wilcoxon's tests but did not examine overall pooled effects across studies [8]. Moreover, metaanalysis pooling estimates across studies is not practical if studies utilize some methods such as Wilcoxon's tests or Random Forest (RF) as there are no test statistics or estimates that are applicable for pooling.
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 nonparametric 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 microbiome 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. More recently, 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 [21,22]. However, these methods are applicable to absolute abundance, which are count data and largely vary between studies due to the difference in experimental and computational approaches making the estimates from these methods incomparable across studies. Therefore, standard meta-analysis approaches may not be directly applicable to pool estimates and their standard errors from these methods across studies.
Here, we developed an R package "metamicrobiomeR" that applies Generalized Additive Models for Location, Scale and Shape (GAMLSS) [23] with a zero-inflated beta (BEZI) family (GAMLSS-BEZI) for the analysis of microbial taxonomy relative abundance data. With BEZI family, this model allows direct and proper examination of microbiome relative abundance data, which resemble a zero-inflated beta distribution. This model also allows adjustment for covariates and can be used for both non-longitudinal and longitudinal studies. We evaluated the performance of GAMLSS-BEZI using simulation studies and real microbiome data. Importantly, the estimates from GAMLSS-BEZI are log(odds ratio) of relative abundances between groups and thus are comparable across microbiome studies and can be directly combined using standard meta-analysis approaches. As such, we apply random effect 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 effect meta-analysis component on four studies of the infant gut microbiome.

GAMLSS-BEZI for the analysis of bacterial taxa relative abundance and bacterial 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 [24]. 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 [25]. A filtered data matrix is then modeled by GAMLSS-BEZI and (mu) logit link and other default options using the R package gamlss version 5.0-5 [23]. A subject random effect can be added to the model for longitudinal data. For performance evaluation, LM and LM with arcsin squareroot transformation (LMAS) were also implemented in the function taxa.compare.

Meta-analysis across studies using random effects models
The adjusted estimates for bacterial taxa or pathway relative abundances from the GAMLSS models are log(odds ratio) of relative abundances between groups and thus are comparable 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 (95% CI) and heterogeneity across studies. A fixed effect meta-analysis model is also implemented for comparison. Meta-analyses 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",st udylab="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 pvalues and multiple testing adjusted pooled p-values of all covariates for each bacterial taxon or pathway.

Additional microbiome measures
Random Forest (RF) modeling of gut microbiota maturity has been widely used to characterize development of the microbiome over chronological time [12,14,26]. Adapting from the original approach of Subramanian et al [14], in the microbiomeage function, relative abundances of bacterial genera that were detected in the Bangladesh data [14] and in the data of other studies to be included were regressed against infant chronological age using a RF model on a predefined training dataset of the Bangladesh study. This predefined training set includes 249 samples collected monthly from birth to 2 years of age from 11 Bangladeshi healthy singleton infants. The RF training model fit based on relative abundances of these shared bacterial genera was then used to predict infant age on the test data of the Bangladesh study and the data of each other study to be included. The predicted infant age based on relative abundances of these shared bacterial genera in each study is referred to as gut microbiota age (Additional file 1).
All implemented functions in the metamicrobiomeR package are summarized in Table 1. alpha.compare Calculate average alpha diversity indexes for a specific rarefaction depth, standardize and compare alpha diversity indexes between groups microbiomeage Predict microbiota age using Random Forest model based on relative abundances of bacterial genera shared with the Bangladesh study [14] read.multi Read multiple files in a path to R

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 LM 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 relative abundance data were based on the R package "gamlss.dist" version 5.0-3 [27].

Type I error
Three scenarios were simulated to mimic typical case-control microbiome studies of small  (Figure 1).

Performance of GAMLSS-BEZI: application to real microbiome data
Using data from a study of Bangladeshi infants [14], which included 996 stool samples collected monthly from birth to 2 years of life in 50 subjects, we demonstrate that GAMLSS-BEZI outperforms LMAS in detecting differential relative abundances between various grouping variables.

Example 1: Comparison between non-exclusively breastfed (non-EBF) vs. exclusively breastfed
(EBF) infants ≤ 6 months of age Figure 2 showed 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 was observed in non-EBF versus EBF infants. GAMLSS-BEZI was able to detect a significant difference in all four of these phyla whereas LMAS could only detect a significant difference in three phyla (Table 3). Table 3. Results of GAMLSS-BEZI and LMAS: real microbiome data example 1 (non-EBF vs. EBF).  GAMLSS-BEZI detected all three of these differences whereas LMEM could only detect a significant difference in one phylum (Table 4).

GAMLSS-BEZI LMAS
Example 1 and 2 demonstrate the increased sensitivity of GAMLSS-BEZI in detecting bacterial taxa with differential relative abundances as compared to LMAS.  Figure 4 showed 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 was observed in the groups of infants with diarrhea vs. those without diarrhea (Figure 4, upper panel). GAMLSS-BEZI detected 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 was observed between those with diarrhea vs. those without diarrhea ( Figure 4, lower panel) and GAMLSS-BEZI did 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.

Illustration of meta-analysis workflow with real microbiome data from four studies
We show examples for the comparison of gut microbiomes between male vs. female infants ≤ 6 months of age adjusting for feeding status and infant age at time of sample collection. We illustrate the meta-analyses across microbiome studies with data from four studies (total number of stool samples=610 [female=339, male=271]). These four studies include one in Bangladesh [14], one in Haiti [11], and two in the USA (California and Florida [CA_FL] and North Carolina [NC]) [12,28] (Table 6, Additional file 1). Longitudinal stool sample collection of 6 healthy full term infants with varied age.
Total number of samples: 21 (female=14, male=7) V1-2 /Roche GS FLX Titanium * This study also contain 674 stool samples > 6 months of age which were used in the analyses comparing the performance of GAMLSS-BEZI (Generalized Additive Models for Location, Scale and Shape (GAMLSS) with a zero inflated beta (BEZI) family) vs. LMAS (linear model with arcsin squareroot transformation). Data from this study was downloaded from the authors' website: https://gordonlab.wustl.edu/Subramanian_6_14/Nature_2014_Processed_16S_rRNA_datasets.html. Data from three other studies were obtained directly from the investigators.

Alpha diversity
For each study, the alpha.compare function imports the outputs from "alpha_rarefaction.py" QIIME1 script and calculates mean alpha diversity for different indices for each sample based on a user defined rarefaction depth. Mean alpha diversity indexes are standardized to have a mean of 0 and standard deviation of 1 to make these measures comparable across studies. Standardized alpha diversity indexes are compared between groups adjusting for covariates using LM. Metaanalysis across studies is then done and the results are displayed as a standard meta-analysis forest plot (Additional file 1). Our results showed that alpha diversity (four commonly used indexes Shannon, Phylogenetic diversity whole tree, Observed species, Chao1) was not different between male and female infants ≤6 months of age in the meta-analysis of the four included studies (pooled standardized alpha diversity (Shannon index) difference (DD)= -0.01 standard deviation (sd), 95% CI=[-0.15; 0.12], p-value=0.8288) ( Figure 5, Additional file 1).

Microbiota age
A RF model based on the relative abundance of the shared bacterial genera of the four included studies explained 96% of the variance related to chronological age in the training set and 67% of the variance related to chronological age in the test set of Bangladesh data. This performance is better than the original RF model proposed by Subramanian et al [14]. The predicted infant age in each included study based on relative abundance of the shared gut bacterial genera using this RF model is referred to as gut microbiota age. Standardized gut microbiota age are compared between groups adjusting for covariates using LM and meta-analysis across studies is then done.
Our results showed that standardized microbiota age was significantly different between genders but in opposite directions in two studies with small sample sizes (Haiti and North Carolina).
However, meta-analysis of all four studies revealed no significant difference in gut microbiota age between genders after adjusting for feeding status and infant age at time of sample collection

Bacterial taxa relative abundance
The adjusted estimates (log(odds ratio) of relative abundance between genders) from GAMLSS-BEZI for each bacterial taxon and the pooled adjusted estimates across studies (meta-analysis) are displayed as a heatmap (Figure 7, left panel). The adjacent forest plot (right panel) displays the pooled adjusted estimates and their 95%CI with different colors and shapes to reflect the magnitude of pooled p-values. Our meta-analyses of four included studies showed that, after adjusting for feeding status and infant age at time of sample collection, there was no significant difference in bacterial taxa relative abundance from phylum to genus levels between male vs. female infants ≤6 months of age after adjusting for multiple testing except for genus Coprococcus whose relative abundance was significantly higher in male vs. female infants (FDR adjusted pooled p-value<0.0001) (Figure 7, Additional file 1).

Bacterial functional (KEGG) pathway relative abundance
Meta-analysis of four included studies showed no difference in relative abundances of bacterial KEGG pathways at level 2 and level 3 between genders after adjusting for feeding status and infant age at time of sample collection and adjusting for multiple testing (Figure 8, Additional file 1).
Difference in gut microbial composition between genders has been reported in adults [29][30][31] and in some neonatal studies albeit with small sample sizes [32,33]. However, the reported findings have varied between these studies. Our meta-analyses of four studies showed virtually no difference in gut microbiomes 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 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 [34,35], and the detected difference in our study may provide some insights into the known sex biases in health outcomes.

CONCLUSION
Our metamicrobiomeR package implemented GAMLSS-BEZI for analysis of microbiome relative abundance data and random effect meta-analysis models for meta-analysis across microbiome studies. The advantages of GAMLSS-BEZI are: 1) it directly and properly 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 between groups and thus are directly comparable across studies.
Moreover, standardization of alpha diversity indexes and predicted microbiota age also yields comparable metrics for cross-study comparisons. Random effects meta-analysis models can be directly applied to pool these comparable 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 difference in gut microbial alpha diversity (standardized Shannon index) between male vs. female infants ≤ 6 months old from each study and the pooled effect across four included studies (meta-analysis). Estimates for diversity difference and corresponding standard error from each study were from linear mixed effect models (for longitudinal data) or linear models (for nonlongitudinal data) and were adjusted for feeding status and age of infants at sample collection.  The difference in gut (standardized) microbiota age between male vs. female infants ≤ 6 months old from each study and the pooled effect across four included studies (meta-analysis). Estimates for standardized microbiota age difference and corresponding standard error from each study were from linear mixed effect models (for longitudinal data) or linear models (for nonlongitudinal data) and were adjusted for feeding status and age of infants at sample collection.  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 pooled estimates across all studies with 95% confidence intervals (95% CI) (forest plot). B: Genus level: heatmap of log(OR) of relative abundances of all gut bacterial genera between male vs. female infants for each study and pooled estimates across all studies with 95% CI (forest plot). 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.  Heatmap of log(odds ratio) (log(OR)) of relative abundances of gut microbial KEGG pathways at level 3 between male vs. female infants for each study and pooled estimates of all studies with 95% confidence intervals (forest plot). Only pathways with pooled p-value <0.3 are shown. 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% confidence intervals (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.