A hierarchical Bayesian approach for detecting global microbiome associations

2.1 Log ratio model

Let us assume that we are collecting microbiome data from N subjects, and in each sample we are counting the occurrence of S species or taxons. We treat the microbiome data as compositional, i.e., even though we are observing count data, we only consider the proportion of each species. If we observe an S-dimensional count vector z _i for sample i, then the proportion of species s can be estimated as ρ i , s = z i , s ∑ s ′ S z i , s ′ . We can now transform the species proportions for species s ∈ {1, …, S − 1} by using the log-ratio transform (Aitchison 1982):

We assume that the y _i are approximately multivariate Gaussian distributed, and model them as:

In our application, Σ is informed by the phylogenetic tree inferred from the sequencing data. Briefly, under the assumption of a Brownian motion evolutionary model, the covariance of two species is equivalent to the branch length from the root to the most recent common ancestor in the phylogenetic tree. If no phylogenetic tree is available, then we set Σ = I, the (S − 1) × (S − 1) identity matrix. The hyperparameter c, c > 0 is a scaling parameter controlling the desired fit to the data. For the data in this study we found that c = 1 led to satisfactory models that showed no evidence of over- or under-fitting.

2.2 Hierarchical regression model

Having established the log ratio model, we now want to formulate μ _i as a linear expression with information sharing across the species (or response variables). Let:

where x _i is a vector of p observed covariates for subject i, B is an (S − 1) × p matrix of positive linear coefficients and Γ is an (S − 1) × p matrix with entries in {−1, 1}. Estimation of global effects is achieved via a prior on the elements of β _j , the S − 1-dimensional vectors corresponding to the columns in B:

where a _j denotes the global effect of covariate j on all species. A graphical representation of the resulting hierarchical model can be seen in Figure 1. In Eq. (4), N + denotes the positive part of the Gaussian distribution truncated at zero. Note that the global effect only defines the prior for the strictly positive magnitude β _s,j of the species-specific effects. Placing a prior on the signed effect size would not have the desired effect; since the size of each biological sample is fixed, an increase in some species must lead to a decrease in other species, and so the global effect size under such an alternative model would not be reflective of the total change across species.

Figure 1:

Graphical representation of the hierarchical regression model for detecting microbiome associations. Shaded nodes are observed values, unshaded nodes are parameters. Here, y _i denotes the log-ratio transformed species abundances in sample i, Σ is the covariance matrix informed by the phylogenetic relationships, x _i denotes the observed covariates, and a _j denotes the global effect of covariate j on y _i . Please see the text for the remaining parameters, as well as the detailed description of the model.

We also need to define a distribution for the elements γ _s,j of Γ. The natural approach would be to treat them as binary variables following a Bernoulli distribution. However, this requires us to perform inference over discrete variables, which can be onerous, and does not lend itself well to gradient-based inference techniques. Instead, we opt for a Gumbel-Softmax approximation to the Bernoulli distribution (Jang et al. 2016; Maddison et al. 2016). The Gumbel-Softmax is a continuous relaxation of the categorical distribution. In general, for a categorical random variable h with class probabilities π ₁, π ₂, …, π _K , we can represent h as:

where one_hot denotes the one-hot encoding of a categorical variable as a K-dimensional vector (hence the change from h to h). The g _k are sampled from a Gumbel(0,1) distribution. Note that this representation is exactly equivalent to a categorical distribution; to obtain the continuous relaxation we need to replace the argmax by a softmax function, which gives:

for k ∈ {1, …, K}, where δ ∈ Δ ^K−1, the (K − 1)-dimensional simplex. Here τ is a tuning parameter that controls the smoothness, with τ = 0 reverting to the discrete case. In our case, K = 2, which means that Eq. (6) simplifies to:

where π ₋₁ = 1 − π ₁. We drop the subscript k as one value of δ fully defines the two-dimensional vector on the simplex, but introduce the subscript (s, j) to denote the value of δ for species s and covariate j. We use the subscripts 1 and −1 to denote Gumbel draws for γ _s,j = 1 and γ _s,j = −1 respectively. Since δ _s,j is either 0 or 1, we set:

This has the desired property that if δ _s,j ≈ 1, γ _s,j ≈ 1 and if δ _s,j ≈ 0, γ _s,j ≈ −1.

We specify truncated Gaussian N + ( 0,1 ) priors for the global effects a _j , which also serves to regularize the model when the sample size is smaller than the number of parameters, which is the case for almost all the models studied here. We further place a Beta prior on the parameter π ₁, where the hyperparameters can be used to encode any prior beliefs about the prevalence of positive or negative associations; for our study we set both equal to 1.

2.3 Inference

We are interested in performing fully Bayesian posterior inference over the parameters of this model. A convenient framework is Hamiltonian Monte Carlo (HMC), which allows us to take advantage of the differentiability of the Gumbel-Softmax to estimate gradients and speed up the Bayesian inference. We employ the Stan sampler to perform HMC sampling; for more details on the Stan sampler and modelling language see Carpenter et al. (2017).

In microbiome research, it is common to encounter datasets consisting of large numbers of operational taxonomic units (OTUs). Despite the efficiency of HMC sampling, sampling hundreds or thousands of parameters β _s,j and γ _s,j is not currently feasible on a personal computer. We therefore propose an approximate estimation, as an alternative to the full log-ratio model, where we only estimate β _s,j as part of the sampling, and the γ _s,j parameters are estimated in a separate inference step. The approximate inference procedure is as follows:

1. Fit independent linear models y i , s = x i T β s * + ϵ i , s , ϵ i , s ∼ N ( 0 , σ s 2 ) to the log-transformed microbiome counts for all species s. Note that β s * is distinct from the rows of B in the full model.

2. Set γ s , j = sign ( β s , j * ) where the sign() operation returns −1 if β s , j * is negative and 1 otherwise.

We refer to the model with Γ inferred from the data as the full log-ratio model, and to the model with Γ fixed according to the two-step procedure above as the conditional model, since inference is conditional on the estimate of Γ. We can think of the conditional model as an approximation to the full joint model P(Γ, B, θ |Y, X) = P(B, θ |Γ, Y, X)P(Γ|Y, X), where θ collects all parameters not in B or Γ, and X and Y represent the log-ratio response and the design matrix with the observed covariates, respectively. Instead of estimating P(Γ|Y, X), which would require marginalising over all other parameters, we use Γ ̂ = sign ( β ̂ * ) where β ̂ * = ( β ̂ 1 * , … , β ̂ S * ) is the maximum likelihood estimate of the effect sizes in independent linear models of P ( y s | x s , β s * ) for each s. Inference is performed on the conditional model P ( B , θ | Γ ̂ , Y , X ) , which is less computationally expensive.

3.1 Setup

We perform a simulation study under different scenarios to test the performance of the full and conditional log-ratio models. In particular, we study how well the full model can predict the global effect (a _j ) and species-specific effect (β _s,j ) parameters under scenarios with different sample size n, variance parameter σ ², and number of species (or taxonomic units) S and compare the performance to that of the conditional model. In each experiment the variables that are not changing are set to n = 100, σ ² = 0.1 and S = 10. The global effect sizes in the simulation model are set fixed with a _j ∈ {0, 0.5, 1}.

For each unique setting of the simulation parameters, we generate 20 random datasets according to the following procedure. First, a random phylogenetic tree is generated using the ape R package (Paradis and Schliep 2018). Based on this tree, the covariance matrix Σ under the Brownian motion evolutionary model is calculated. Next, we generate an N × p matrix of covariates X with x i , j ∼ N ( 0,1 ) . For each covariate j and species s, γ _s,j are generated from a binomial distribution with p = 0.5. Species-specific coefficients β _s,j are generated according to Eq. (4). Then we simply calculate μ _i according to Eq. (3), and sample y _i according to Eq. (2). The ρ _i,s must sum to 1, and, for s ∈ {1, …, S − 1}, must be equal to ρ _i,S  exp(y _i,s ) by Eq. (1). To achieve this we set ρ i , S = 1 / ∑ s S − 1 exp ( y i , s ) . We can then calculate the ρ _i,s and simulate count data z _i from a multinomial distribution with parameters ρ _i , and number of trials equal to 100*S.

The bench-marking platform used for this study ran R-3.6.1 to generate datasets and invoke rstan. Models were fitted using the Stan software (Carpenter et al. 2017) version 2.19.1. A comparison between popular methods such as the PERMANOVA (generalized uniFrac method using the GUniFrac R package, Chen and Chen 2018), ANOSIM (non-parametric multivariate analysis using the vegan R package, Clarke 1993) and MiRKAT (microbiome regression-based kernel association test using the MiRKAT R package, Zhao et al. 2015) has been done. The sampling was executed in parallel using 4 different chains each with 8000 iteration with a stopping criteria at convergence on 4 Intel Ivy Bridge cores each running at 1.15 GHz speed and 12 GB of RAM memory in total.

3.2 Results

In this section we investigate and compare the performance of:

1. the full log-ratio model,

2. the conditional log-ratio model which uses a linear marginal model for each species to determine the species-specific effect signs γ _s,j and then infers the size of the species-specific effects β _s,j and the global effect a _j conditional on the γ _s,j ,

3. the ”oracle conditional model”, where the γ _s,j parameters are set to the true values used for simulating the data.

Figure 2 shows the performance of the models for estimating the global effect size in three different experiments: changing the sample size, noise variance on the β _s,j parameters, and number of species. Dashed lines represent the true values of the global effect (a _j ) in the simulation model. It can be seen that increasing the number of samples n improves the accuracy of prediction and reduces the variance. On the other hand, increasing the variance parameter for the effect values and increasing the number of species (or taxonomic units) both make estimation of the global effect parameters harder. We see that there is not much difference between the performance of the full model and that of the oracle conditional model, but as could be expected, the full log-ratio model produces less biased estimates compared to the conditional log-ratio model, with the oracle log-ratio model outperforming both. The bias in the conditional log-ratio model effect size estimates can be explained by the inherent model misspecification due to the estimation of γ _s,j via independent linear regressions.

Figure 2:

Performance comparison between the models. (a) Comparison of the performance of the full, conditional and oracle conditional log-ratio models on simulated data when varying the sample size (top row), noise variance of the species-specific effects (middle row) and number of species or taxonomic units (bottom row). The linear simulation model consists of the intercept term and two covariates with global effect sizes 1, 0.5 and 0 respectively. The dashed horizontal lines show the size of the true effect, and the boxplots represent the posterior mean estimates obtained from 20 independent simulated datasets. (b) Comparison of the performance of the full, conditional and oracle conditional log-ratio models on simulated data when estimating the species-specific effect size (β _s,j ). Each point corresponds to the estimated effect size for one species in one model. The linear simulation model consists of the intercept term and two covariates with global effect sizes 1, 0.5 and 0 respectively. The solid lines represent y = x (note that the x-axes are on different scales).

Figure 2 also reveals that even for large sample sizes, a small bias remains when estimating the global effect size for covariate 1, whose simulated value was set to zero. We observe that the simulated species-specific effect sizes tend to be small but non-zero. In the absence of a strong regularizing prior on the global effect, this leads to a small non-zero global effect. We chose not to impose such a prior to avoid over-regularising.

Figure 2b shows a comparison between the performance of the three models when predicting the species-specific effect size β _s,j . We note that species-specific effects for covariate 1 and the intercept are predicted well across all three models. For covariate 2, we again see a significant bias in the estimation for the conditional log-ratio model. Looking at both Figure 2a and b we can see that in most cases the full log-ratio model outperforms the conditional model.

We have also checked the sensitivity of our model to differences in the variance parameters for the species-specific effect sizes β _s,j (see Appendix A), and note that despite the assumption of a common variance parameter, our model is robust to varying noise levels across species. Furthermore, in Figure 6c (see Appendix A), we can see the distribution of the estimated species-specific effects when changing the number of species. It can be seen that both models converge to a small non-zero value, however, in the full model sampler converged in a few iterations but the conditional model took longer to explore.

Figure 3b shows the p-values obtained using the PERMANOVA (generalized UniFrac), ANOSIM and MiRKAT methods applied to simulated data with n = 100 samples, σ ² = 0.1 and S = 10. The number of covariates increases from p = 10 to p = 60. We note that empirically, for a large number of covariates (in this case, when the number of covariates exceeds 30), the statistical power becomes insufficient to reject the null hypothesis of no global effect.

Figure 3:

Efficiency and estimation performance of the proposed models. (a) Comparison of the computational time to convergence of the full and conditional log-ratio models, on the simulated data when increasing sample size n and number of species (or taxonomic units) S. The error bars show the 95% confidence interval calculated over 20 independent simulated dataset. (b) Comparison of the performance of the ANOSIM, MiRKAT and PERMANOVA models on simulated data when increasing the number of covariates. (c) Comparison of the performance of the full and conditional models on simulated data when increasing the number of covariates. The y axis represents the estimated global effects with the dashed line representing the mean value.

The PERMANOVA, ANOSIM and MiRKAT models quantify the effect of a covariate via the variability in microbial compositions, and do not try to estimate a global effect size. Hence we cannot directly compare with our log-ratio models. Figure 3c shows the performance of the full and conditional log-ratio models when increasing the number of covariates. We see that increasing the number of covariates results in a slight negative bias in the estimated effect size. However, even for the maximum of p = 60 covariates, the estimated global effect size remains far from zero under our models, indicating a lower false negative rate compared to the methods in Figure 3b.

We also investigate the computational efficacy of the full and conditional log-ratio models. Figure 3a shows the run time to convergence under each of the models. We see that the sampler for the conditional model takes longer for large S; this is a consequence of setting the γ _s,j parameters fixed which results in less efficient exploration of the sample space and longer convergence times. The sudden drop in time for S ≥ 100 is explained by the fast convergence of the chains in the HMC sampler to a small non-zero value (when the true global effect value is set to zero). A comparison is shown in Figure 7, where we can see the estimated species-specific effects when changing the number of species, across all different HMC chains. It can be seen that when S ≥ 100 the sampler in the full model converged in a few iterations but the conditional model took longer to explore.

Taken together, Figures 2a, b and 3a indicate that for small sample sizes (n < 100) one could use either the full or conditional log-ratio models, as long as the number of covariates is not too large. For large sample sizes the conditional model could be deployed to improve computational efficiency, but there is a trade-off in terms of estimation accuracy. For large number of species S (S > 50), the full log-ratio model seems more efficient than the conditional model, perhaps because it can take advantage of Hamiltonian paths in the posterior landscape that are closed off when conditioning on Γ. In general, we would recommend to use the full log-ratio model where possible, and only resort to the conditional log-ratio model when resource constraints and large sample sizes make use of the full model impractical.

4.1 Detecting associations

Using metabolomics as a measure of health, or to predict response to pharmacologics is a very attractive prospect, as measuring metabolites in urine, for example, may provide a non-invasive and potentially inexpensive screening method that has potential for application to a plethora of pathologies.

We apply the full log-ratio model to a dataset collected as part of a pilot study into inflammatory bowel disease (IBD) (Beamish 2017). IBD is an umbrella term used to describe two conditions (Crohn’s disease and ulcerative colitis) that are characterized by chronic inflammation of the gastrointestinal tract. 42 gut bacteria samples were obtained from patients undergoing colonoscopy at Lancashire Teaching Hospitals Trust (LTHTR) and University Hospitals of Morecambe Bay Trust (UHMBT), and metagenomic profiling was performed to obtain the bacterial species counts. The patients were additionally asked to provide a urine sample, and metabolite concentrations found in those samples were recorded. One objective of the study was to investigate associations between the gut microbiome and the urinary metabolite profiles.

After quality control, our cohort consists of n = 40 people (IBD cases and controls), for which we have measured a set of 3 clinical covariates: age, sex, and the hospital where the sample was collected. Note that we do not investigate any associations with IBD status, as number of cases was too small to make any reasonable inferences. Instead, we focus on detecting associations between the microbiome and individual metabolites. We consider the set of the S = 21 most abundant operational taxonomic units (OTUs or species) at the phylum level. As a result, our design matrix X is a 40 × 5 matrix consisting of the intercept, the metabolite level, age, sex and sample collection location. As the data is limited, we do not try to infer the variance parameter σ ², but keep it constant and equal to 0.1. Finally, we have a total of 261 metabolites.

Figure 4a shows the posterior distribution of estimated global effect size a _j for the metabolites with the largest posterior mean effect estimate. We see that out of these 21 metabolites, saccharopine, allantoin, vanillate and ornithine have the strongest association with the microbiome. Some metabolites show interesting bimodal patterns in their posterior effect distribution; this is particularly obvious for glycine. Such patterns could potentially be indicative of more than one mechanism that links the microbiome to metabolite levels.

Figure 4:

Application of the log-ratio model to a microbiome-metabolome association study (Beamish 2017). (a) Density plots showing the posterior distribution of the global effect size associating each metabolite to changes in the microbiome. (b) Posterior mean estimates of the phylum-specific effect size for each metabolite across the different bacterial phyla.

Figure 4b shows the posterior mean estimate of the species-specific effect sizes β _s,j , specifically the association of each phylum with each metabolite. We note that for the metabolites with the largest posterior effects, the majority of species-specific effects are strongly negative, indicating that these metabolites are associated with a decrease in most phylum counts, compensated by an increase in just a couple of phyla.

We note that our analysis identified several phyla-associated metabolites that have previously been shown to be linked to markers of intestinal health or nutrition, demonstrating potential impacts on health that are reflective of microbiota changes. First, our study suggests that adipate is strongly associated with changes in the microbiome. Adipate, the salts and esters of adipic acid, are increased in the urine of breast-fed versus formula-fed neonates (Dessì et al. 2016), a now well-known essential factor in the establishment of a diverse, healthy microbiome. We observe effect sizes that are indicative of an association between adipate and an increase in Firmicutes, Proteobacteria and Bacteroidetes, the major components of a healthy adult intestinal flora. Therefore, clinical measurement of adipate may have potential as biomarkers for intestinal health.

Epicatechin is a polyphenolic flavonoid abundantly present in plants, cocoa and that offers protection against diabetes by its ability to mimic the effects of insulin and has strong antioxidant properties (Kanehisa and Goto 2000; Samarghandian et al. 2017). The association we observed between epicatechin and proteobacteria abundance may prove useful in the prediction of pathologies associated with proteobacteria expansion, such as IBD or nonalcoholic fatty liver disease (Rizzatti et al. 2017). Of course, any potential mechanisms that underpin these associations will require further investigation but collectively this demonstrates the broad and significant potential for clinical application of our model. Several metabolites, including isoleucine, thymidine, ornithine and vanillate, were associated with an increase in cyanobacteria (Figure 4b). Cyanobacteria are unusual intestinal microflora (Ley et al. 2005), in that they perform oxygenetic photosynthesis and cyanobacteria-based supplements, such as Spirulina, regarded as functional foods, have been shown to support immunological function (Finamore et al. 2017). These associations may allow us to get a clearer picture of the functionality of the gut microbiome, by assessing the functions of these metabolites.

4.2 Estimating the global effect of diet on the gut microbiome

It is well-known that diet has a profound effect on the composition of the gut microbiome (David et al. 2014; Singh et al. 2017). However it is difficult to quantify this effect using traditional microbiome analysis methods, as we either need to look at individual species, or rely on a pre-specified distance metric. In Turnbaugh et al. (2009), the authors took the latter approach to characterise the differences in microbiota in mice whose gut tracts were colonised by donor bacteria of human or mouse origin, and who received either a Western or low-fat plant-based (LF/PP) diet. The dataset is described in Table 1. Donor bacteria either originated from human microbiome samples (fresh or frozen), or from mice who had previously been implanted with (fresh) human samples (HMouse) and had received either the LF/PP or Western diet. The dataset consists of n = 444 fecal samples, obtained at various time points from 15 mice, and we use the bacterial abundance data for S = 160 taxonomic units at the genus level. Note that we did not have access to the information on sample times or sample source, and so were unable to take this aspect into account in our model.

Using the UniFrac distance measure, Turnbaugh et al. (2009) found that recipient mice consuming a Western diet showed a similar microbiome composition that was distinct from mice on the LF/PP diet, regardless of donor. We applied the full log-ratio model to the dataset to determine whether this was reflected in the global microbiome effect of diet on the microbiome, and to study the individual (genus-specific) effects. We include sex, donor type and diet as covariates. The reference taxonomic unit for the log-ratio transform was the Dorea genus, although investigation of alternative reference taxonomic units showed little change in the results Appendix.

The results of our study are shown in Figure 5. We observe that overall there is a significant (non-zero) effect of diet, and a strong effect of using the frozen samples. This is concordant with the findings in Turnbaugh et al. (2009). We also note that receiving a bacterial implant from a donor mouse on the LF/PP diet seems to have a stronger effect than receiving an implant from a donor mouse on the Western diet, compared to having a human donor. Ruminococcaceae Incerta Sedis is strongly associated with the change in diet, a finding which is backed up by the literature (Clarke et al. 2013). Similarly, Eubacteria colonization of the gut has been shown to be affected by diet (Simmering et al. 2002). Finally, we see that Enterobacteria is associated with the change in diet; previous studies in rats have shown that a high-fat diet results in considerably more propionate and acetate producing species, including Enterobacteria (Singh et al. 2017). Interestingly, this bacterial genus also shows a strong association with implantation by frozen donor bacterial samples.

Figure 5:

Application of the log-ratio model to a gut microbiome study (Turnbaugh et al. 2009). The reference taxonomic unit for the log-ratio transformation is Dorea. (a) Density plots showing the posterior distribution of the global effect size associating each covariate to changes in the microbiome. Note that the reference category for diet is LF/PP, meaning that the global effect is the change with respect to that diet. Similarly the reference category for bacterial donor is fresh (human), and the reference category for sex is female. (b) Posterior mean estimates of the genus-specific effect size for each metabolite across the different bacterial genera. We display here the 10 OTUs with the highest effect size for the association with diet.