A Bayesian framework for identifying consistent patterns of microbial abundance between body sites

2.1 Experimental design

A study suitable to answer the previously described research question can be broken down into the following steps. First, an appropriate subject population exhibiting the disease or target phenotype is identified and biological samples from the two body sites of interest are collected. Multiple samples from the same patient within and across body sites are possible, but not necessarily required. Next, sample preparation and 16S rRNA sequencing are performed. This sequencing technique aims to identify and count hypervariable DNA patterns that are specific to microbial species and OTUs, but that do not exist in human DNA; the rationale being that the DNA content of a group of microbes is approximately proportional to their abundance in the sample. So, by counting how often signatures belonging to a specific OTU are observed, we can obtain an estimate of its abundance relative to how many microbes were observed, in total. After OTUs have been counted, our proposed statistical test is performed individually for each OTU, testing the null hypothesis that there is no PASTA relationship for each specific set of considered microbes. This test is performed by fitting a statistical regression model to the row vector of abundance values corresponding to a target OTU, followed by the calculation of a test statistic Tθ based on the parameter estimates (for example, rate of absence or mean abundance) obtained from said model. This statistic will be small when H0 is true and large when H0 is false. An overview of the experimental design for testing the hypothesis of PASTA can be found in Figure 1.

Figure 1:

Overview of the experimental setup to test for pairwise stratified association (PASTA). Oral and gut samples are obtained from cancer patients and 16S rRNA sequencing is performed on each sample. The resulting microbial abundance data is used to fit a statistical regression model to each observed OTU across all samples. Finally, abundance estimates across strata are used to test whether abundance patterns in disease status are preserved between mouth and gut.

2.2 Data model

In what follows we consider abundance on two taxonomic levels: the genus and the Amplicon Sequence Variant (ASV) level, the latter representing unique biological sequences that were identified from 16S genes (Callahan, McMurdie & Holmes, 2017). In order to make abundance values comparable across samples and bring them to the same scale, raw counts are first transformed into relative abundance values. For a given sample, relative abundance of an OTU refers to the number of times that OTU was observed, scaled by the total number of observed OTUs for that sample. It represents the proportion of times an OTU was observed in a given sample.

Let Yk denote the relative abundance of a specific OTU for sample k. This response can be modeled as a Zero-inflated Beta distribution with probability density fYk(y|pk,ωk,ϕk). This model assumes that the case Yk=0 occurs with probability pk and that given Yk>0, the response Yk follows a Beta distribution with mean ωk and dispersion ϕk. For a given OTU and sample, the probability of absence p defines how likely it is to observe no microbe comprising that OTU within said sample. The mean non-zero relative abundance ω represents the mean relative abundance given that microbes comprising the OTU are actually observed. The mean of Yk, the overall mean relative abundance, is then E[Yk]=μk=ωk(1−pk). The probability density function of this distribution can be expressed as follows:

Before the statistical PASTA test can be performed, a Bayesian ZIBR model is fit to the data Yk utilizing the likelihood fYk and assuming a common dispersion parameter ϕk=ϕ for all samples. The estimated posterior distributions of ω and p resulting from this model are then subsequently used to conduct the test.

Let ω denote the vector of mean relative abundances for all samples, p denote the vector of probabilities of absence for all samples, β,δ denote coefficient vectors, b,d denote random effect vectors and Q,R,W,X represent design matrices. The design matrices code how covariates impact the model parameters via the following link functions:

For our application, the matrices W and X are used to model the strata of body site and disease status, but can additionally be used to adjust for other, fixed covariates (e.g. subject age, gender, smoking status, and/or other potential confounders or sources of variation). On the other hand, the optional inclusion of Q and R permits one to account for correlation structures, such as within-subject correlation when multiple samples are collected from the same patient. If, for example, the probability of absence of a given OTU in sample k from subject jk is assumed to be impacted by disease status gjk, body site sk, age Ajk and within-subject correlation, we would formulate W, Q and our model for the probability of absence as follows:

Here, δ1,gjk,sk captures the effect of disease status gjk and body site sk on the probability of absence, δ2 represents the effect of age on the probability of absence, and d1,jk is a subject specific random intercept. Analogously, if mean relative abundance is believed to be impacted by the same covariates in the same way, except that age was believed to have no effect, we would formulate X, R and our model for non-zero relative abundance such that:

Here, β1,gjk,sk captures the effect of disease status gjk and body site sk on the mean non-zero relative abundance and b1,jk is a subject specific random intercept. The specific models used in our analysis of the pancreatic cancer dataset are presented in the Results Section are further discussed in Section 2.7.

Let t,i∈N denote placeholder indices for any potential coefficient specified in the above ZIBR model. Then all posterior distributions were estimated, utilizing the following independent prior distributions:

Priors were chosen to be weakly informative, with the exception of ϕ being restricted to values larger than or equal to 1 in an attempt to stabilize estimation of means. Under these priors and for some integer vectors T,I1,I2 the posterior distribution of parameters will then satisfy the following:

There are generally no analytical solutions for the posterior distributions of the coefficients when random effects are present. Regardless, whether the model structure is a special case that allows for analytical calculation of posterior distributions or whether we employ a more complex model where this is not possible, posterior distributions can be estimated via Markov chain Monte Carlo (MCMC) methods. Briefly, MCMC procedures allow one to draw arbitrarily large samples from a posterior distribution that will numerically approximate the said distribution, as the number of draws increases. Models were fit via this method in the software OpenBUGS (version 3.2.3 rev 1012) via the R (version 3.4.0) package “R2OpenBUGS” (version 3.2.3.2).

2.3 Formal definition of pairwise stratified association (PASTA)

In order to understand how the Bayesian regression model can be used to conduct the desired hypothesis test, we will first provide a formal definition of PASTA. Let s denote a grouping variable for which two groups are to be compared. For our purposes, this grouping variable represents body sites: s=1 denotes gut and s=2 denotes mouth. Let g denote another grouping variable with three or more distinct categories. This grouping variable will represent different types of disease status, more specifically cancer-subtype. Let θsg be a population parameter of the response for a given body site s and disease status g. The population parameter represents fundamental properties of the distribution of the response. For the here considered ZIBR model, p,ω and μ are relevant candidates for θ. If PASTA holds for a given OTU, then either p,ω or μ will associate between the two body sites, because they all relate to the magnitude of abundance.

We thus define: The parameter θ exhibits PASTA with respect to s and g if there exists an increasing function h(x) such that θ1g=h(θ2g) holds for all g∈{1,2,…,G}, where G≥3. Conceptually, this definition says that as we move from one disease status group to another, if θ increases in oral samples, it will also increase in gut samples. Analogously, if θ decreases from one disease status group to another in the mouth, it will also decrease in the gut. A visualization is provided in Figure 2.

Figure 2:

Visualization of pairwise stratified association (PASTA). Let θ represent a population parameter of interest, for example the mean relative abundance of a particular OTU. Each column of sub-figures below the table are examples of a PASTA relationship, i.e. of h being an increasing function. The first row plots parameter values of mouth and gut side-by-side and demonstrates that a variety of different scenarios are covered by this definition. In the second row, plotting parameter values of gut against parameter values of mouth reveals their association through a trend. T denotes Pearson correlation values between gut and mouth.

2.4 Testing for PASTA

Let T(x,y)∈[−1,1] denote a correlation statistic between two numerical vectors x,y; for example, the Pearson or Spearman correlation statistic. Under this definition, Tθ=T(θ1,θ2) denotes the correlation statistic calculated for the two parameter vectors corresponding to s=1 (e.g. parameters for disease status groups in the mouth) and s=2 (e.g. parameters for disease status groups in the gut). Generally, if a PASTA relationship holds between θ1 and θ2, this statistic should assume a larger value compared to cases where such a relationship does not hold. This means that we are able to formulate our desired test by rejecting H0 if Tθ is larger than a specific threshold and fail to reject H0 if it is less than said threshold.

In summary, assume that θ1g=h(θ2g) implies Tθ>tc for some −1<tc<1. The constant tc represents a meaningful degree of association. For example, a value of tc=0 would mean that any tangible degree of association is meaningful, where a value of tc=0.5 would mean that a moderate degree of association is meaningful. This definition is useful because Tθ can score the degree of association without explicitly having to specify the shape of h. Considering the complexity of the biology underlying the samples, specifying h in advance may not only be hard to justify, but strong deviations of a chosen h from the true h could also result in missing promising associations. Instead, our regression model will allow each stratum (s,g) to have an independent effect on the response, leading to a unique, agnostic posterior distribution of each θsg. These unique posteriors are then in turn used to calculate the posterior distribution of Tθ and conduct the hypothesis test.

Based on this scoring definition of PASTA, we formulate our hypotheses in the following way:

• H0:Tθ≤tc, i.e. θ1 and ​θ2 do NOT exhibit PASTA

• H1:Tθ>tc, i.e. θ1 and θ2 DO exhibit PASTA

While deriving analytical solutions of the distribution of Tθ|Y=T(θ1|Y,θ2|Y) will depend on the data model and may be difficult or even impossible to obtain depending on the modeling scenario, a general testing procedure can still be derived. As decribed earlier, MCMC methods allow one to conveniently obtain a large sample of posterior draws of each θsg, even when obtaining analytical solutions of posterior distributions is not possible. Furthermore, plugging the posterior draws of each MCMC iteration into T allows one to obtain posterior draws from Tθ itself. Let α denote the target credibility threshold, H0 is then rejected if the lower bound tQα of the one-sided credible interval of Tθ|Y exceeds tc. This is equivalent to rejecting H0 if the estimated probability of no association exceeds α, i.e. Pr(Tθ|Y≤tc)<α. In detail, the step by step process for testing H0 is as follows:

1. Specify a likelihood for the response data Y and prior distributions for the parameters 𝜽.

2. Utilize a MCMC sampling scheme to draw a large number of samples from the posterior distributions of the parameters θ|Y. One draw from the Markov Chain contains a unique draw for each θsg.

3. Calculate Tv∗=T(θ1∗v,θ2∗v) where θ∗v denotes the vth MCMC draw. Then T∗ is a large sample of the posterior distribution Tθ|Y.

4. Calculate the α⋅100% sample quantile tQα of T∗. If the Markov Chain is sufficiently long, the sample quantiles of T∗ will closely approximate the quantiles of the true posterior distribution. The value tQα is thus the lower bound of the (1−α)⋅100% one-sided credible interval of Tθ|Y.

5. Reject H0​ if the lower bound tQα is larger than tc.

This process is generally applicable regardless of the data model or the parameter being tested, as long as each θsg can be estimated without constraining them to a parameter space that implies PASTA.

2.5 Pancreatic cancer patient dataset

In order to evaluate validity of the approach in the context of microbiome data, analyses were performed based on a biological 16S rRNA sequencing dataset first published in del Castillo et al. (2019). This dataset contained samples of various gut and oral sites from 77 patients with pancreatic cancer with age range 31–86 years. Sequencing was performed utilizing the Illumina MiSeq System and alignments were performed using BLASTN against a reference library combining sequences from HOMD (version 14.5), Greengenes Gold and the NCBI 16S rRNA reference sequence set. OTU counts were obtained utilizing the QIIME (Quantitative Insights Into Microbial Ecology 22) software package version 1.9.1, while the unique Amplicon Sequence Variant (ASV) counts were calculated using the QIIME2 software package release 2018.4. The former was used to obtain taxonomic genus level counts, whereas the latter was used to obtain rarefied ASV level information, calculated based on sequencing data rarefied at a sampling depth of 1200. Both genus level and ASV level counts were considered for analysis. Before fitting statistical models to the data, relative abundance values of less than 0.01 were treated as noise and set to 0. To ensure inference was based on sufficient signal, OTUs and ASVs were only tested if more than 5of all samples exhibited non-zero values.

The dataset was used to both guide simulation studies (described in the next section) and to deploy models to identify potential microbes that may exhibit a PASTA pattern.

2.6 Simulation studies

Before simulations were performed, an empirical approach was pursued in order to obtain sampling distributions of the parameters p,ω,ϕ that would be representative of biological microbiome data. First, a marginal, unstratified ZIBR model was fit to the pancreatic cancer dataset that assumed all samples of relative abundance for a given OTU originated from the same distribution. These model fits yielded a single estimate of p,ω and ϕ for each OTU. These estimates were then assumed to be representative of or approximate the true distribution of parameters in biological data. In the next step, the estimates were used to obtain smooth probability distributions that parameters could be sampled from during the simulation studies. For both p and ω, individual Beta distribution models were fit to the marginal estimates in order to obtain their smooth sampling distributions. On the other hand, log⁡ϕ was sampled via a Normal distribution through an observed linear relationship between log⁡ϕ and log⁡ω that was present on the ASV level and the genus level. More specifically, since our models assumed fixed dispersion among all groups, dispersion was sampled from (log⁡ϕ|minsg{log⁡ωsg})∼N(aminsg{log⁡ωsg}+b,σ2), where a,b,σ2 differed between genus and ASV level.

After the smooth sampling distributions were obtained, the performance of PASTA tests was evaluated via simulations. Let t denote a target, fixed degree of association, n denote the number of observations in each stratum (s,g) and tc=0 denote the tested degree of association. A single simulation run was carried out by first randomly drawing all θsg parameters from the representative sampling distributions, until |Tθ−t|<0.001 was satisfied. This process yields parameters that are both representative and that also exhibit a target degree of association (within a small error margin). Next, the drawn parameters satisfying this condition were plugged into the likelihood of the ZIBR data model, which was in turn used to draw a random sample of relative abundance values. This simulated pseudo-data was then used to fit the Bayesian ZIBR model and conduct our hypothesis test. Each considered scenario was simulated 1000 times and statistical power for given t,n and tc was then estimated as the proportion of times H0 (i.e. Tθ≤0) was rejected. We specifically considered Pearson correlation as choice for T(x,y) in this simulation.

An additional restriction was put in place for sampling pseudo-data in order to prevent rare cases of sparse datasets with insufficient signal to perform the analysis. If a generated pseudo-dataset contained more than three sub-strata (s,g) in which all observations exhibit a response value of either all 0 or all 1, then it was rejected and a new pseudo-dataset was sampled.

2.7 Model fitting

Let jk denote the unique identifier index for the subject and sk denote the body site that sample k originated from. Also, let and gjk denote the disease status for subject jk, let Xk be the log of total sample abundance for sample k and let bjk denote the random intercept for subject jk. The three different models that were utilized in this study are shown below:

1. logit(ωk)=βsk,gjk&logit(pk)=βsk,gjk

2. logit(ωk)=βsk,gjk+bjk&logit(pk)=βsk,gjk+bjk

3. logit(ωk)=βsk,gjk+bjk&logit(pk)=β1,sk,gjk+Xkβ2+bjk

Model A was utilized in the simulation studies. Model B was utilized for fitting ASV level data, while Model C was utilized for fitting genus level data. This choice was made because scaling OTU counts to relative abundance will only make non-zero relative abundance comparable between samples, but not the rate of absence. This is due to the fact that, even if the true probability of absence p for a specific OTU is very high, if more microbes are overall observed in sample 1 than in sample 2, then the probability of observing none of the microbes belonging to the target OTU in sample 1 is much lower than in sample 2. For example, if a total of 1,000,000 microbes live in a body site and 100 of them belong to the genus Prevotella, then if we randomly extract 1000 microbes from this body site with our sample, we would expect to only rarely find one of these 100 microbes in our sample. However, if our sample randomly extracts 100,000 microbes from the body site, it would be rare to find none of the 100 microbes in it that belong to the genus Prevotella. So since the genus level data was not rarefied, the total sample abundance differed between samples and an adjustment was necessary, whereas the ASV level data was rarefied and did not require adjustment for total sample abundance.

In order to achieve potentially better convergence behaviour and to simplify and speed up the model fitting, the logistic regression component of the model was fit independently of the Beta regression component, in all cases. The resulting posterior chains of p and ω were then used to calculate the posterior chain of ω. This approach is justified under the assumption that p and ω are independent after adjusting for covariates, but may be inadequate when there are confounders affecting both parameters not accounted for in the model.

3.1 Simulation studies

Performance of our proposed approach was first evaluated using series of simulation studies. In an attempt to obtain sampling distributions of parameters that would approximate biological distributions, unstratified ZIBR models were fit to each OTU in the pancreatic cancer dataset (see Methods for details of this dataset). Unstratified parameter estimates were then used to obtain smooth sampling distributions of ω,p,ϕ. Finally, these sampling distributions were used to generate many pseudo-datasets satisfying H1 and performance was evaluated when applying the previously described testing approach to the simulated dataset.

Sampling distributions for parameters were similar for both genus and ASV level. However, for ω, the mean non-zero relative abundance, distributions tended to be slightly further concentrated toward 0.0 on the ASV level as compared to the genus level. Further, distributions of p tended to be slightly more concentrated toward 1.0 on the ASV level as compared to the genus level. In both cases a linear relationship was observed between log⁡ω and log⁡ϕ which was ultimately used to sample ϕ conditionally on ω (Figure 3).

Figure 3:

Observed relationships between marginal distributions of ω and ϕ estimated from the pancreatic cancer dataset. For both the genus and the ASV level, parameters were estimated marginally for each OTU across all observations without any stratification. When plotting marginal parameter estimates of ω and ϕ a linear relationship can be observed on the log scale. This relationship was utilized to sample ϕ conditionally on ω in the simulation studies.

In summary, the following sampling distributions were obtained:

As expected, simulations of biological data revealed that analyses on the genus level were overall more powerful than on the ASV level, regardless of which population parameter was investigated (Figure 4). Assuming tc=0, four disease status groups, 95% credible intervals and utilizing Pearson correlation, the highest power was achieved when testing PASTA of ω. Under a moderate degree of association of Tθ=0.537 a target power of 0.8 was reached for 5 samples per stratum on the genus level and 15 samples per stratum on the ASV level. Type 1 error rates appeared adequately calibrated to the 5% significance level ranging from 0.03 to 0.056 on the genus level and from 0.032 to 0.06 on the ASV level. Despite the relatively modest within group sample size needed to detect a moderate degree of association for ω with adequate statistical power, there appeared to be considerably less power for tests of p. Under a high degree of association of Tθ=0.834 a target power of 0.8 was reached for 40 samples per stratum on the genus level. On the ASV level, utilizing as many as 80 samples per stratum resulted in a power of only 0.59 for the same Tθ. Type 1 error rates also appeared mostly calibrated in this scenario, but showed deflatation for smaller sample sizes, assuming a value of 0.027 on the genus level and 0.009 on the ASV level.

Figure 4:

Results of the simulation studies. Power plots are displayed for testing PASTA of various population parameters with tc=0 at both ASV and genus level. The term “n per group” refers to the number of samples available in each of the eight sub-group combinations resulting from two body sites and four different levels of disease status. H0 was rejected if Pr(Tθ|Y≤0)<0.05. Type 1 error rates are displayed in white colored boxes with black fonts. Power values less than 0.8 are colored blue, values larger than 0.9 are colored red and values between 0.8 and 0.9 are colored orange. Genus level pseudo data generally has higher statistical power than the ASV level. High performance is achieved by the non-zero mean ω, while an increased sample size is required for the probability of absence p. Tests of the overall mean μ result in low performance, when only mildly constraining sparsity.

Testing PASTA of the overall mean μ=ω(1−p) was also investigated. While improving with increasing size of effect and sample size, the power for this parameter was lower than when considering ω,p individually. Even when considering the large degree of association Tθ=0.834 and using 100 samples per stratum, the genus level scenario achieved a power of only 0.546. Notably, type 1 error rates were consistently deflated, ranging from 0.005 to 0.018 on the genus level and 0.002–0.008 on the ASV level. Type 1 error rates were deflated across all simulated scenarios, reaching values of less than or equal 0.018 or less.

Discrepancies in performance were found to be directly related to precision of parameter estimates. When plotting the posterior means of Tθ against their true simulated values across various simulation runs, the variation around the identity line consistently increased from ω to p, aswell as from genus to ASV level (Figure 5). Analogously, posterior distributions of Tθ were found to on average become more diffuse and more biased towards 0, when moving from ω to p or from genus to ASV level. When performing simulation runs of a scenario with low relative precision, in which p was sampled from a Uniform(0.85,0.95) distribution, the posterior distribution of Tθ was on average almost perfectly centered at zero and highly diffuse.

Figure 5:

Effects of the relative precision of parameter estimates on the posterior distribution of Tθ. The first row shows the average point estimate of posterior quantiles of Tθ across simulation runs for various simulation scenarios. The second row shows the associated plots of the parameters’ posterior means versus their true values across simulation runs. As the relative precision of parameter estimates decreases, the posterior distribution of Tθ becomes more diffuse and more biased towards 0.

Poor power when testing μ was also found to be related to two additional factors. Detailed results of simulations accounting for these factors are displayed in Additional File 1. The deflated type 1 error rates when utilizing 95% credible intervals, lead to overly conservative tests that negatively affected power. Calibrating type 1 errors to 5% by adjusting lower bounds of the credible intervals of Tθ for each considered sample size, lead to a consistent improvement in power, reaching a value of 0.73 for Tθ=0.922 and 80 samples per stratum on the genus level. The second factor that affected performance was the employed liberal three sub-strata rule, which allowed up to three strata to exhibit exclusively zeroes. Since high rates of absence were simulated, this case often naturally occurred leading to the three respective posterior estimates being imputed with the vague prior distribution, which is very imprecise. A follow-up simulation restricting all strata to have at least one non-zero observation, lead to a consistent increase in power, reaching a value of 0.82 for Tθ=0.922 and 80 samples per stratum on the genus level. In both settings, overall performance for testing μ was consistently lower than for testing p regardless of taxonomic levels. When both calibrating type 1 error and restricting non-zero observations at the same time, power increased further but was not consistently better than for testing p.

Additionally, three sets of simplified supplementary simulations were also performed to showcase how the PASTA testing approach can be analogously utilized in other data models. A simple Poisson regression model and a log ratio Aitchison model both achieved performance metrics slightly less performant but overall comparable to testing PASTA of the mean non-zero relative abundance via Beta regression. In particular, the Beta regression model appeared to achieve higher power for small sample sizes than the other two approaches and the log-ratio Aitchison model appeared to perform slightly worse than the Poisson regression model. On the other hand, testing PASTA of the overall mean in a zero-inflated Poisson model, utilizing the same smooth sampling distribution of zero-inflation rate p as for the ZIBR model on the genus level, achieved performance metrics comparable to testing PASTA of the overall mean in the ZIBR model. Even though minor differences with respect to calibration of type 1 error and statistical power were observed across the different models, the testing approach was overall viable regardless of the scenario. A detailed summary of these simulations is provided in Additional File 2.

3.2 Applying the approach to biological data

The 10th revision of the International Statistical Classification of Diseases and Related Health Problems (ICD-10), is a systematic classification of medical conditions provided by the World Health organization. Disease status information in this dataset was available via ICD-10 codes for each subject. The clinical pancreatic cancer dataset contained four predominant cancer-types: C24.x, C25.x, K86.2 and other, where “.x” denotes a further sub-type that could differ by subject and “other” refers to pancreatic cancer in various other categories or other diseases of the foregut. OTUs exhibiting significant PASTA with respect to disease status were successfully identified for both genus and ASV level in this dataset. For analysis we considered coding cancer sub-type in two ways: four group coding as described above and three group coding, which collapsed “K86.2” and “other” into one group. On the genus level, when coding disease status into four groups three genera exhibiting PASTA between mouth and gut were identified: Fusobacterium, Haemophilus and Veillonella (Table 2). After substratifying oral sites into saliva, tongue, buccal and gum, these association were found to be preserved for some of the site pairs: Fusobacterium also exhibited PASTA between gut an saliva sites; Haemophilus also exhibited PASTA between gut and gum, as well as gut and tongue; Veillonella also exhibited PASTA between gut and gum. Several genera also exhibited PASTA between individual mouth sites (Table 3). Two genera exhibited PASTA between four pairs of mouth sites: Fusobacterium and Actinomyces. Three genera exhibited PASTA between two pairs of mouth sites: Atopobium, Haemophilus and Prevotella. Six genera exhibited PASTA in only one pair of mouth sites.

On the ASV level, two ASVs exhibited PASTA with respect to p between mouth and gut when disease status was coded into four groups. When coding disease status into three groups, the same two ASVs as before and three additional ASVs exhibited PASTA with respect to p between mouth and gut. Notably, among these additional ASVs was a candidate belonging to the Fusobacterium genus. Further details of the ASV level analysis are discussed in Chung et al. (2019).