Statistical modeling on human microbiome sequencing data

1.   Introduction

The field of genomics has been developed to conduct metagenome of the microbiota over the last two decades ^[1,2,3,4]. Current microbiome studies are mostly motivated by the research topics in which aim to understand the relationship among microbiome, host, and genetic or environmental factors. A variety of studies has identified the association between microbiome and host ^[5,6], and these works have shown how human microbiota affects health and diseases. For example, microbial changes are proved to be linked with Parkinson's disease ^[7], inflammatory bowel disease ^[8,9], diabetes ^[10], and cancers ^[11]. The other field of microbiome studies is to examine the association between microbiome and genetic/environmental variables, in particular, the effects of the interested genome ^[12] and environment covariates ^[13,14] on specific microbiome composition. Research has been conducted on skin conditions ^[15,16], obesity ^[17], and immunity system ^[18].

The microbiome data is quantified by amplified sequences using generic sequence similarity, produced via next-generation sequencing of the 16S ribosomal ribonucleic acid (rRNA) gene ^[19] and bioinformatic pipelines such as QIIME ^[20]. Classification on the sequencing reads based on phylogenetic levels (genus, family, suborder, order, subclass, class, phyla, kingdom, and domain) is referred to as operational taxonomic unit (OTU) counts (Figure 1). These OTU counts provide the foundation and development of research on microbiome while generating challenges of statistical methods due to the features. The OTU counts are usually zero-inflated and over-dispersion since the particular host's microbiome taxonomy reads are context specific (Figure 2) ^[21]. In addition, both the hierarchical structure of microbiome and the sample collection of multiple measurements from related individuals both yields the correlation between different taxa. Another feature in many microbiome studies is that the number of subjects can be smaller than the number of taxa to be explored, which also imposes difficulties in statistical modeling ^[22,23,24].

Figure 1.  An illustration of the taxonomic structure of OTUs.

DownLoad: Full-Size Img PowerPoint

Figure 2.  A histogram of the counts of a single OTU to exemplify the distribution of microbial data.

DownLoad: Full-Size Img PowerPoint

This paper seeks to summarize the recent development regarding statistical methods to analyze human microbiome sequencing data, including assessment and comparison between current popular models and advanced technique for longitudinal microbiome involving serial correlations and hierarchical structure. Diversity among microbiome communities is identified using non-parametric methods. Furthermore, the association between microbial components and genome has been explored, especially including genetic variants on X-chromosome. Finally, potential directions are discussed to consider for further research within this field, such as prediction using microbiome data and classification of OTUs based on the research outcome of interest.

2.   Methods

The relationship between the abundance of single or multiple OTUs and environmental or genetic factors has been investigated extensively yet without standard statistical methods ^[12,13,14]. Classical models such as linear regression and logistic regression models are the most popular approaches, while facing the risk of the violation of normality and constant variance assumptions for linear models, and the loss of information on zero parts hence lower statistical power on (generalized) logistic regression models ^[25]. Non-parametric models such as Wilcoxon rank sum (WRS) test can be used as an alternative approach without the normality assumption, but it cannot adjust for covariate effects ^[26]. Standard t-tests are used for comparison between two relative abundance datasets ^[27,28,29]. When the number of communities is more than two, the one-way analysis of variance (ANOVA) or the Kruskal-Wallis test become appropriate. However, none of these methods above can account for the excess zeros, which is a key feature of microbiome components data, yet currently are still widely used for microbiome studies.

2.1.   Parametric models for Zero-inflated data

In order to deal with excess zeros, mixture models are proposed such as zero inflated (ZI) models and hurdle models (also called two-part models). Typical microbiome data contains OTU counts which can be referred to as the variable $Y_{ij}$ for OTU $j$ in sample $i$. Since the mixture models deal with single OTU, we can ignore the OTU index $j$ in this section. ZI models are usually a mixture of a Poisson (ZIP) or Negative Binomial (ZINB) model with a point mass at zero. These models process data with excessive zeros in two steps. First, structure zeros are distinguished from counts data using a Bernoulli trial with probability $\phi_i$ for each $y_i$. Then the counts data is assumed to follow a Poisson or NB distribution. Specifically, the probability mass function (PMF) of a zero-inflated model for $y_i$ can be written as:

Hurdle models contain two parts: One part is a binomial model to indicate whether a zero or non-zero outcome occurs, and the other part is a truncated-at-zero model for count parts only. In particular, the PMF of the hurdle model for Possion or NB model is written as:

and the $\phi_i$ can be linked to covariates through logit link

where $\boldsymbol{W_i}$ denotes the vector of covariates for $\phi_i$. The two sets of models treat zero part distinctly. ZI models allow structural zeros in the point mass and leave sampling zeros in the parametric model part, while hurdle models handle structural zeros and sampling zeros together in the binomial model part.

Xu. et al. ^[30] performed a comprehensive comparison of different model fitting performance for zero inflated data in human microbiome studies, especially in the aspects of type Ⅰ error and statistical power of the tests. The performance of parameter estimations such as accuracy and efficiency on count parts and zero parts as well as the goodness of fit of the models were also evaluated. Simulations revealed that both hurdle models and ZI models provide a better model fit than standard one part models, resulting in less biased and more efficient parameter estimations while controlling for type Ⅰ errors and maximizing power.

2.2.   Bayesian latent variable models for hierarchical clustered data with repeated measures

The model introduced in the last section follows independent data assumption, which does not always hold for general microbiome data due to the hierarchical clusters and repeated measurements. Xu et al. ^[31] proposed a Bayesian latent variable (BLV) model, which is joint modeling of multiple taxa within a single taxonomic cluster. This model can not only make inference on the genetic and environmental risk factors’ effects within the cluster, but also account for repeated measurements of the microbiome from related family members.

The BLV model can incorporate multiple response variables to account for underlying correlations among the multiple taxa within the taxonomic cluster. Let $\boldsymbol{Y}_{cit} = (y_{cit1}, ..., y_{citJ})^T$ be the $J \times 1$ vector of outcomes measured at the tth time point on the ith individual from the cth family, for $t = 1, ..., T_{ci}$, $i = 1, ..., N_c$, $c = 1, ..., C$. $C$ denotes the total number of families; $N_c$ denotes the number of family members in the cth family; $T_{ci}$ denotes the total number of repeated measurements for the ith individual from the cth family. Let $\boldsymbol{U}_{ci} = (U_{ci1}, ..., U_{ciT_{ci}})^T$ and $\boldsymbol{U}_{z, ci} = (U_{z, ci1}, ..., U_{z, ciT_{ci}})^T$ be the two vectors for the longitudinal latent trait underlying the count components and the structural zero parts.

Note that the model can accommodate a mixture of distributions such as NB and ZINB. This is important in microbiome research since the multivariate OTU outcomes can be mixed counts with and without a zero-inflated feature. Consequently, a portion of $\boldsymbol{Y}_{cit}$ follows $NB(\tau_j, \mu_{citj})$ and the other part of the outcomes follows ZINB with PMF

The BLV model consists of two parts. The first part is the measurement model, which models the latent traits and a portion of the covariates effects directly on the responses using a generalized linear mixed model. In particular, the parameters depend on both latent traits $\boldsymbol{U}$ and $\boldsymbol{U_z}$ and some other covariates $\boldsymbol{W}$ and $\boldsymbol{W_z}$ through canonical links $log(\mu_{citj}) = \gamma_{0j} + \boldsymbol{W}_{cit}^T \gamma_j + \lambda_j U_{cit} + b_{cij}$ and $log(\frac{\phi_{citj}}{1-\phi_{citj}}) = \beta_{0j} + \boldsymbol{W}_{z, cit}^T \beta_j + \lambda_{z, j} U_{z, cit} + b_{z, cij}$. $b_{(z), cij}$ is the family-specific, within-subject random effect and both $b_{cij}$ and $b_{z, cij}$ are assumed follow normal distribution with mean zero and variance $\eta_j^2$ and $\eta_{z, j}^2$, respectively. The second part, known as the structural model, shows the relations between the other portion of covariates and the latent traits with a linear mixed model. Specifically, the second part is written as

and

where $\boldsymbol{X_{cit}}$ and $\boldsymbol{X_{z, cit}}$ are the covariates of interest on the latent traits. $g_c$ and $g_{z, c}$ are the environmental random effects for $\boldsymbol{U_{cit}}$ and $\boldsymbol{U_{z, cit}}$ within a family. $\boldsymbol{\alpha_{c}}$ and $\boldsymbol{\alpha_{z, c}}$ are the additive genetic random effects. $\boldsymbol{d_{ci}}$ and $\boldsymbol{d_{z, ci}}$ are the serial random effects. Hence, the risk factors in this model have indirect effects on the response through the latent variable. This Bayesian method is applied for parameter inference on the complicated form of the posterior distribution. Polya-Gamma data augmentation (PGDA) technique is used as a part of the Markov chain Monte Carlo (MCMC) algorithm for the model.

Simulation studies showed that both the direct and indirect effect parameters have small bias and the root mean square errors (RMSEs), suggesting that the proposed model performs well with controlled type Ⅰ error and reasonable power of tests. Random effects estimations mostly are unbiased, except for the variance of serial random effects on the probability of structural zeros and the LVs.

2.3.   Distance based method to access differences between communities

Parametric models introduced in the previous sections can provide parameter estimations and statistical inference on the risk factor association. However, it is difficult to identify and select the optimal distribution of the model, and sometimes the model assumptions may be violated due to complex microbiome data structure. Other than parametric models, distance-based non-parametric models provide another approach to analyze differences between microbiome communities as a whole rather than only incorporate a univariate outcome. In Shestopaloff's paper ^[32], a mixture model has been proposed to model OTUs generated from Poisson models with subject-specific underlying rates. The observed counts are assumed to follow Poisson distribution $X\sim Poisson(r_i)$ and $r_i = q_i N_i$, in which $q_i$ is individual specific relative abundance sampled from underlying population distribution $G_q$, and $N_i$ is the total reads in the sample.

Due to the complexity of the observed data, a set of distribution components is proposed to model the underlying population distribution, including a zero point mass, a set of left-skewed distribution for low rates, a set of Gamma distributions based on the posterior of the Poisson rate $\lambda|n \sim \Gamma(n+1, 1)$, and a truncation point mass $P(X > C) = 1$ for a sufficiently large number C to account for sparsity. The final mixture distribution depends on the specified M gamma distribution $(\alpha_m, \beta_m)$ with weights $\boldsymbol{w} = (w_z, w_1, ..., w_M, w_{C+})$ respectively. The least square optimization function is applied to estimate the optimal weights by taking the difference between the observed and expected agammaegate counts.

The choices of parameters for low rate structure are decided by selecting the model among a set of nested models using a nonparametric bootstrap based on the minimum distance between the expected agammaegate counts and the observed counts. This model selection technique is similar to cross-validation as a way to fit the model with a portion of data to avoid overfitting. The joint mixture distribution estimate $\vec{w} = \sum_l v(l) \vec{w}_l$ can be calculated using the weight $v(l)$ of each model obtained from the bootstrap. Then the probability of observing $n_i$ from the $m^{th}$ Gamma mixture conditioning on the estimated weight and resolution $t_i$ is

And the probabilities of the two point masses are $P(i \in G_z) = I(n_i = 0)$ and $P(i \in G_{C+} = I(n_i > C))$.

Finally, permutation tests can be performed to evaluate the significance of differences between communities after the calculation of pairwise distances. These distances can be treated as the pairwise $L^2$-PDF norms $D_{L^2}(i_1, i_2) = \|P_{i_1} - P_{i_2}\|$ for $i_1$, $i_2 = 1, ..., I$ and $i_1 \neq i_2$, where the PDF for each sample $i$ is $P_i = [P_i(Z), P_i(0), ..., P_i(C), P_i(C+)]$. The permutation method is introduced by Anderson ^[33].

The simulation study showed that the mixture model correctly estimated the true underlying rate distribution and the proportion of structure zeros. When the true model followed certain parametric models, the corresponding parametric methods such as 2P-LOLS and NBH performed well. Non-parametric Kruskal-Wallis tests gave a relatively robust estimation, too. However, these methods can incorporate only a single outcome. When the comparison of multiple OTUs is of interest, distance-based methods become the only applicable approach, and the proposed mixture model performed better than the Manhattan distance method for most scenarios.

2.4.   X-chromosome association with microbial composition

Besides the research works examined the association between microbiome and diseases, many researchers have interests in the relationship between host factors and intestinal microbial composition. The host factors include age, gender, dietary, and genetic variants. However, among all the studies proved the linkage of the host genome and human microbiome ^{[12,34,35,36]}, none of the research investigated the impact of genetic variants on X-chromosome on the microbiota. Espin-Garcia et al. ^[37] proposed a finite mixture model (FMM) for the analysis of the X-chromosome and the composition of the microbiota. The method incorporates multiple unknown underlying X-chromosome mechanisms (XCMs), including random X-chromosome inactivation (XCI), skewed XCI (XCI-S), and escape of XCI (XCI-E).

This regression-based model consists of two parts: a distribution function and a mixing proportion of XCMs. Let $y_i$ be the OTU counts or the indicator function for the presence/absence of the OTU for subject $i$ depending on the outcome of interest. Let $x_i^k = (1, s_i, g_i^{s, k}, w'_i)$ be the covariates vector for subject $i$ under $k$th mechanism, where $s_i$ be the sex indicator; $g_i^{s, k}$ is the coded SNPs; and $w_i$ is the vector for additional covariates. Then the likelihood for the assumed observed data with a pseudo mechanism is

where $p_k$ specifies the mixing proportion for the $k$th mechanism. $f(\cdot)$ is one of the two distribution functions - either a zero-inflated probability distribution using Poisson or negative binomial distributions for count parts, or a two-part model treating zero and count parts separately.

Comprehensive simulations were constructed to compare the performance of the proposed X-chromosome model with some existing ones such as a Clayton-like model ^[38] and a PLINK-like model ^[39], as well as strategies which assume the same mechanisms for all subjects. Expectation-maximization (EM) algorithm was used to estimate the potential genetic effects, and a score statistic was computed to evaluate hypothesis testing. In conclusion, the FMM provides relatively less biased estimates and competitively higher power while controlling for type Ⅰ error in the comparison of other methods.

2.5.   Association of host genome, microbial composition, and intestinal permeability

As implementation of the statistical models on real human microbiome studies, the heritability of microbial components are explored using a log-normal model with a generalized estimating equation (GEE) algorithm among 270 related individuals, and successfully identified 94 of 249 OTUs with significant additive genetic components with high heritability ^[12]. This result suggested that host genetics is strongly associated with the intestinal microbial composition. In addition, Genome-wide association studies (GWASs) were conducted to identify an association between genetic variants and bacterial taxa. GEE framework was adopted to identify the association between genetic polymorphisms and the relative abundance of heritable taxa, controlling for age, sex, and the top three genetic principal components. The two-part log-normal model was fitted on zero counts and nonzero counts separately by using logistic regression and log-normal model, respectively. External validation confirmed four specific OTUs associated with the host genetic variations.

Analysis of intestinal permeability (IP) is also conducted in healthy first-degree relatives of individuals with Crohn's disease ^[40]. A generalized least squares model was applied to evaluate the association between fecal microbiota and IP. Potential clustering within families was accounted by a compound symmetry correlation matrix. No significant correlation between several levels of bacterial taxa (phylum, family and genus) and IP was found, adjusted by clinical factors such as age, age squared, gender, and the province of origin.

3.   Future research works

In the last few years, the introduced statistical methods have been performed more frequently to model the zero-inflated, over-dispersion, and clustered microbiome data, and investigate the association of microbiota and phenotypes or genotypes. In the following session, some potential directions are introduced for prediction and classification purposes using microbial composition.

3.1.   Prediction of phenotypes based on human microbiome data

Statistical testing methods which are traditionally used for identifying associated OTUs can be performed for prediction purposes, but they are limited by the correctness of regression models. Particularly for microbiome data with multiple measurements for each subject, generalized linear mixed models (GLMMs) are an effective way to estimate the OTUs' effects on clinical outcomes while accounting for both fixed and random effects. By setting certain thresholds for the p-values controlled by false discovery rate (FDR), GLMMs for OTU counts given a phenotype and other covariates can investigate not only the association between microbial composition and the phenotype, but also the potential predictors in which we are interested. Two-part models and zero-inflated models can be adopted to GLMMs to account for excessive zeros.

However, disagreements often occur between statistical significance and predictivity. These limitations lead to the usage of statistical and machine learning methods such as the least absolute shrinkage and selection operator (LASSO) for prediction purposes. LASSO serves feature selection goal by adding the $L_1$ of the regression coefficients as a penalty term to the log likelihood function to achieve the intent of shrinkage. The objective log likelihood function for LASSO given observations ${(y_i, x_i), i = 1, ..., n}$ and a tuning parameter $\lambda$ is

where $\| \beta_k \|_1$ is the $L_1$ penalty of $\beta_k$. Therefore, fitting the model for a phenotype given OTU composition with the LASSO method could be performed to predict the presence/absence of the phenotype. As both LASSO multinomial logistic regression and GLMMs can select OTUs as potential predictors, a combination of the two models becomes a promising method for prediction of phenotypes by first using GLMMs to screen on all the given OTUs then applying LASSO to the subset of OTUs selected previously.

3.2.   Classification on microbiome data

Classification methods as a supervised learning technique, provide another angle to predict the outcome of interest by building a classifier based on the training dataset and predicting for the unknown observations. Some existing classification methods have been developed extensively, such as logistic regression models, discriminant analysis, classification decision trees, and k-nearest neighbor methods. Besides those methods, the distance-based mixture models introduced previously can be extended for a new classification method on microbiome sequencing data.

3.3.   High-dimensional GLM for graphic structure

Generalized linear models (GLMs) are another way of classification under the high-dimensional setting in which the number of predictors exceeds the number of observations. In particular, logistic regression models are popular for classification in microbiome data. Matson et al. ^[41] showed that relevant microbial compositions are naturally correlated, and thus these components within a cluster reflected the association with an interesting clinical outcome simultaneously. Therefore, consideration of such structure within the neighborhood that can incorporate sparsity pattern would boost the correctness of prediction. A combination of the graph structure and GLM through a node-wise penalty which is able to account for neighborhood sparsity can be applied for classification.

4.   Conclusions

In this paper, we review and summarize current statistical methods to analyze zero-inflated and over-dispersion microbiome sequencing data. Other than the traditional linear models, zero-inflated models and hurdle models are recommended which have advantages of handling excess zeros. In addition, repeated measurements for individuals with their relatives are commonly encountered in microbiome studies. A Bayesian latent variable model along with the PGDA technique is proposed to handle the hierarchical clusters and the longitudinal correlation. Meanwhile, the parametric models can only accommodate single OTU as the outcome, but distance-based methods test the differences between communities across the population with multiple microbial compositions. The other difference between the parametric models and non-parametric approach is that covariates can be adjusted in parametric regression models, while distance-based methods are limited to randomized clinical trials due to the lack of ability to handle the potential confounding variables.

Besides treating OTUs as the outcome of some interested clinical variables, some researchers also investigate the association between genome and microbiome. Although not much work had been done to explore the genetic effects of sex-chromosome on the microbiome, one of the papers which we review focuses on the relationship with X-chromosome genetic variants and microbiota, by incorporating the underlying X-chromosome architecture to the finite mixture model. Some applied analysis has also been conducted, specifically on the association of host genome, microbiome, and intestinal permeability.

Prediction for human health using microbial components is a promising direction. In particular, GLMM and LASSO are the models that have good performance when a subset of OTUs is selected to predict the disease status. Furthermore, classification of microbial communities based on disease outcomes also has a potential impact in the medical field. Distance-based classification algorithm and high-dimensional GLM graphic structure models are worth to investigate further to classify the related OTUs given the presence-absence status.

Acknowledgements

W.X. was funded by Canadian Institutes of Health Research (CIHR Grant 145546)，Natural Sciences and Engineering Research Council of Canada (NSERC Grant RGPIN-2017-06672), Crohn’s and Colitis Canada (CCC Grant CCC-GEMIII), and Helmsley Charitable Trust. D.Y. was supported by NSERC Grant RGPIN-2017-06672, CCC Grant CCC-GEMIII, and Edwin S.H. Leong Scholarship.

Conflict of interest

The authors declare that there are no conflicts of interest.