This vignette reproduces the real-data analysis of paper §5. The data come from the gnotobiotic mouse experiment of Patnode et al. (2019): a 15-member human gut bacterial consortium under a HiSF diet supplemented with different fibres, with two categorical covariates — diet (HiSF / citrus pectin CPT / pea fibre PEF) and the presence/absence of WH2 (B. cellulosilyticus, OTU 2).
library(BCAIA)
data(mice)
dim(mice$Y)
#> [1] 69 15
mice$species
#> [1] "Bacteroides ovatus ATCC 8483"
#> [2] "Bacteroides cellulosilyticus WH2"
#> [3] "Bacteroides thetaiotaomicron 7330"
#> [4] "Bacteroides thetaiotaomicron VPI-5482"
#> [5] "Bacteroides vulgatus ATCC 8482"
#> [6] "Bacteroides caccae"
#> [7] "Bacteroides finegoldii"
#> [8] "Bacteroides massiliensis"
#> [9] "Collinsella aerofaciens"
#> [10] "Escherichia coli"
#> [11] "Odoribacter splanchnicus"
#> [12] "Parabacteroides distasonis"
#> [13] "Ruminococcaceae sp."
#> [14] "Ruminococcus albus"
#> [15] "Subdoligranulum variabile"mice$Xmean and mice$Xcov are the design
matrices used in the paper (Xcov includes an intercept
column).
Fitting the model
The paper uses K = 8 and niter = 160000
(~22 min on an M1 Mac). Here we run a short chain for illustration.
fit <- bcaia(mice$Y, mice$Xmean, mice$Xcov, K = 8,
niter = 160000, burnin = 80000, thin = 10, seed = 1)#> BCAIA fit (simple model)
#> data: n = 69, J = 15, Pmean = 5, Pcov = 4, K = 8
#> MCMC: 1500 iters, burn-in 750, thin 5 -> 150 saved samples
#> runtime: 0.5 minCondition-specific interactions
Each row of Xcov corresponds to a sample; the six
experimental conditions are the unique rows. We estimate the posterior
correlation matrix for each.
Interaction between a pair of species across conditions
For example, OTUs 9 and 13 (Collinsella aerofaciens and Ruminococcaceae), reported as negatively associated across all conditions in the paper:
xlist <- setNames(lapply(seq_len(nrow(conds)), function(i) conds[i, ]),
paste0("cond", seq_len(nrow(conds))))
posterior_cor_pair(fit, j = 9, k = 13, xlist = xlist)
#> setting median lower upper
#> 2.5% cond1 -0.7815852 -0.8825989 -0.5857827
#> 2.5%1 cond2 -0.6776013 -0.8443246 -0.3762429
#> 2.5%2 cond3 -0.7287531 -0.8593776 -0.5290164
#> 2.5%3 cond4 -0.7188846 -0.8690691 -0.5201368
#> 2.5%4 cond5 -0.5714349 -0.7759705 -0.2950505
#> 2.5%5 cond6 -0.6483336 -0.8474767 -0.1984853Differential abundance
Posterior draws of the mean-regression coefficients beta
quantify how abundance shifts with diet and WH2 status (paper Fig. 5).
For example, the posterior of the difference between two
Xmean columns:
beta <- fit$samples$beta # J x Pmean x nsamp
diff <- beta[, 1, ] - beta[, 2, ] # contrast of columns 1 and 2 of Xmean
data.frame(OTU = mice$species,
median = round(apply(diff, 1, median), 2),
lower = round(apply(diff, 1, quantile, 0.025), 2),
upper = round(apply(diff, 1, quantile, 0.975), 2))
#> OTU median lower upper
#> 1 Bacteroides ovatus ATCC 8483 1.22 0.74 1.68
#> 2 Bacteroides cellulosilyticus WH2 -6.28 -6.76 -5.70
#> 3 Bacteroides thetaiotaomicron 7330 2.04 1.47 2.55
#> 4 Bacteroides thetaiotaomicron VPI-5482 0.22 -0.29 0.69
#> 5 Bacteroides vulgatus ATCC 8482 0.66 0.22 1.21
#> 6 Bacteroides caccae 0.00 -0.50 0.44
#> 7 Bacteroides finegoldii 0.24 -0.22 0.91
#> 8 Bacteroides massiliensis 0.42 -0.19 1.03
#> 9 Collinsella aerofaciens 0.05 -0.66 0.58
#> 10 Escherichia coli 0.36 -0.15 0.81
#> 11 Odoribacter splanchnicus 0.67 0.19 1.14
#> 12 Parabacteroides distasonis 0.34 -0.17 0.81
#> 13 Ruminococcaceae sp. 0.19 -0.34 0.62
#> 14 Ruminococcus albus 0.18 -0.28 0.64
#> 15 Subdoligranulum variabile 0.14 -0.35 0.62With the full chain these recover the significant diet and WH2 effects reported in the paper.