Multi-domain skin microbiome analysis • SpBGFM

The skin dataset is the multi-domain chronic wound microbiome study analysed in the paper: 60 swabs from 20 patients, each sampled from the wound before and after a debridement treatment and from a control site of healthy skin. Two domains were profiled on the same samples, 75 bacterial OTUs (genus level) and 39 viral OTUs (host level).

library(SpBGFM)
data(skin)

vapply(skin$Y, dim, integer(2))
#>      bacteria virus
#> [1,]       60    60
#> [2,]       75    39
colSums(skin$X)                      # 20 samples per condition
#>     pre    post healthy 
#>      20      20      20
head(skin$subject)                   # three consecutive samples share a subject
#> [1] 1 1 1 2 2 2
round(vapply(skin$Y, function(y) mean(y == 0), numeric(1)) * 100, 2)
#> bacteria    virus 
#>    42.98    44.10

The columns of skin$X are indicators for the pre-treatment, post-treatment and healthy conditions in that order, so $\beta_{mj1}$ , $\beta_{mj2}$ and $\beta_{mj3}$ are the effects of those three conditions.

Choosing K

choose_K(skin$Y)
#> [1] 36
#> attr(,"eigenvalues")
#>   [1] 5.890715e+02 1.371683e+02 1.042114e+02 7.613328e+01 6.080027e+01
#>   [6] 5.948550e+01 5.229088e+01 4.529022e+01 4.148690e+01 4.040958e+01
#>  [11] 3.441128e+01 3.195269e+01 3.068491e+01 2.836695e+01 2.729709e+01
#>  [16] 2.454983e+01 2.299181e+01 2.221597e+01 2.086358e+01 1.856415e+01
#>  [21] 1.787167e+01 1.556778e+01 1.519613e+01 1.431937e+01 1.407062e+01
#>  [26] 1.302582e+01 1.234215e+01 1.142411e+01 1.098468e+01 1.024979e+01
#>  [31] 9.820081e+00 9.357950e+00 8.636284e+00 8.187287e+00 7.146183e+00
#>  [36] 7.094049e+00 6.778154e+00 6.488957e+00 6.149683e+00 5.962876e+00
#>  [41] 5.539481e+00 4.936403e+00 4.797765e+00 4.189266e+00 3.973561e+00
#>  [46] 3.756996e+00 3.583981e+00 3.150855e+00 2.811937e+00 2.550467e+00
#>  [51] 2.476530e+00 2.252101e+00 1.949996e+00 1.909034e+00 1.645261e+00
#>  [56] 1.513501e+00 1.283116e+00 1.106884e+00 8.984982e-01 5.948100e-14
#>  [61] 3.975522e-14 3.097517e-14 2.902195e-14 2.197534e-14 2.002696e-14
#>  [66] 1.830765e-14 1.625022e-14 1.463503e-14 1.354989e-14 1.347464e-14
#>  [71] 1.234265e-14 1.090685e-14 9.290496e-15 8.760036e-15 7.314356e-15
#>  [76] 6.865099e-15 6.496032e-15 5.921526e-15 5.464713e-15 4.358076e-15
#>  [81] 3.774840e-15 2.941357e-15 1.931147e-15 1.681818e-15 1.431052e-15
#>  [86] 2.470655e-16 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#>  [91] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#>  [96] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [101] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [106] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> [111] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

The paper uses K = 15.

Fitting the model

fit <- spbgfm(skin$Y, subject = skin$subject, X = skin$X,
              K = 15, niter = 100000, thin = 5, seed = 1)
fit

On an Apple M1 laptop this takes roughly five minutes per 10,000 iterations. The remaining chunks are not evaluated when the vignette is built; run them after fitting.

Within- and cross-domain interactions

rho <- posterior_cor(fit)

## the cross-domain block, bacteria against viruses
cross <- posterior_cor(fit, domains = c(1, 2))
dim(cross)

## Pseudomonas (bOTU 59) with Pseudomonas phage (vOTU 18)
cross[59, 18]

## the bOTUs with a strong association with any other OTU
strong <- which(apply(abs(rho) - diag(nrow(rho)), 1, max) > 0.5)
strong[strong <= skin$J[1]]

Covariate effects

Fig 10 of the paper plots the three pairwise contrasts. Interval estimates that exclude zero are flagged in the significant column.

healthy_vs_pre <- posterior_beta(fit, contrast = c(3, 1))
subset(healthy_vs_pre, significant & domain == 1)

post_vs_pre <- posterior_beta(fit, contrast = c(2, 1))
sum(post_vs_pre$significant)     # few changes between pre and post treatment

Overall the bacterial OTUs tend to be enriched in the healthy condition while the viral OTUs are enriched in the wound conditions, and changes between the pre- and post-treatment samples are small — plausibly because the post-treatment swab was taken shortly after the treatment.

Predictive checking

feat <- c(1, 69, skin$J[1] + 17)     # bOTU 1, bOTU 69, vOTU 17
conds <- list(pre = c(1, 0, 0), post = c(0, 1, 0), healthy = c(0, 0, 1))

obs <- normalize_counts(fit, features = feat, r_pred = 0, log = TRUE)
for (q in seq_along(feat)) {
  dens <- lapply(conds, function(x)
    density(log(posterior_predict(fit, feat[q], r_pred = 0, x = x) + 1)))
  plot(dens[[1]], main = paste("feature", feat[q]), xlab = "log(count + 1)")
  lines(dens[[2]], col = "red", lty = 2)
  lines(dens[[3]], col = "blue", lty = 2)
  points(obs[, q], rep(0, nrow(obs)), pch = 4)
}

Reference

Zhang, S., Shen, Y., Chen, I. A. and Lee, J. (2025). Sparse Bayesian Group Factor Model for Feature Interactions in Multiple Count Tables Data. Journal of the American Statistical Association, 120(550), 723–736. doi:10.1080/01621459.2025.2449721