Sparse Bayesian Group Factor Model for Multiple Count Tables • SpBGFM

Sp-BGFM (Sparse Bayesian Group Factor Model for Feature Interactions in Multiple Count Tables Data) is a Bayesian semi-parametric model for multiple domains of next-generation sequencing microbiome abundance data. It models cross-domain interaction between microbial features directly using observed count tables, and provides a flexible DP mixture structure for interpretation.

This repository holds both the R package implementing the method and the scripts that reproduce the paper.

Contact: Shuangjie Zhang (shuangjie.zhang AT austin DOT utexas DOT edu)

Citation

Please cite the paper, not the software:

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

@article{Zhang2025SpBGFM,
  title   = {Sparse {B}ayesian Group Factor Model for Feature Interactions in Multiple Count Tables Data},
  author  = {Zhang, Shuangjie and Shen, Yuning and Chen, Irene A. and Lee, Juhee},
  journal = {Journal of the American Statistical Association},
  year    = {2025},
  volume  = {120},
  number  = {550},
  pages   = {723--736},
  doi     = {10.1080/01621459.2025.2449721}
}

The same entry is returned by citation("SpBGFM").

Installation

# install.packages("devtools")
devtools::install_github("shuang-jie/SP-BGFM")

The package compiles C++ through Rcpp and RcppArmadillo, so a working C++ toolchain is required. On macOS install the Xcode command line tools with xcode-select --install.

Quick start

library(SpBGFM)
data(skin)

# 60 samples from 20 subjects, three conditions, two domains (75 bOTUs, 39 vOTUs)
vapply(skin$Y, dim, integer(2))

choose_K(skin$Y)      # data-driven choice of the number of latent factors

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

rho   <- posterior_cor(fit)                     # within- and cross-domain correlations
cross <- posterior_cor(fit, domains = c(1, 2))  # bacteria against viruses
eff   <- posterior_beta(fit, contrast = c(3, 1))# healthy vs pre-treatment
pred  <- posterior_predict(fit, features = 1, x = c(0, 0, 1))

Y may also be a single concatenated count matrix, in which case pass J, the number of features per domain. Omit X to fit without covariates, and omit subject if each sample comes from a different subject.

Main functions

Function	Purpose
`spbgfm()`	Fit the model by MCMC
`spbgfm_control()`	Hyperparameters and MH tuning constants
`choose_K()`	Data-driven choice of the number of latent factors
`posterior_cor()`, `posterior_Sigma()`	Interaction structure, whole matrix or a domain block
`posterior_beta()`	Covariate effects and contrasts between levels
`posterior_predict()`, `normalize_counts()`	Posterior predictive model checking
`simulate_spbgfm()`	Data-generating processes of the simulation studies
`skin`	The multi-domain chronic wound microbiome dataset

Vignettes reproduce Simulation 1, Simulation 2 and the skin microbiome analysis:

vignette("simulation1", package = "SpBGFM")

Reproducing the paper

The original analysis scripts are kept alongside the package. They are self-contained and do not depend on it.

A note on the sampler

Three defects were found and corrected in the MCMC sampler after the manuscript was prepared: the Metropolis-Hastings ratio for the inner mixture weight omega^alpha in update_w_a mixed units across mixture components and used an inverted indicator, and poswrt scaled the size-factor residual by the kernel variance v_m^2 instead of u_r^2. The scripts in this repository include the corrections.

The effect on the reported results is immaterial. On Simulation 1 with identical data and 100,000 iterations, the corrected sampler gives RMSE of the correlations 0.0330 against the 0.031 reported in Table 1, while adhering an order of magnitude more closely to the mean constraints E[alpha] = nu^alpha and E[r] = nu^r that hold by construction under the model.

The exact code used to produce the published figures and Table 1 is preserved and tagged:

git checkout paper-original         # the code behind the published results
git checkout paper-original-2024    # the original February 2024 upload

Environment setup for the scripts

install.packages(c("Rcpp", "RcppArmadillo", "statmod", "GIGrvg", "extraDistr", "abind", "mvnfast",
                   "LaplacesDemon", "Rfast", "Matrix", "ggplot2", "latex2exp", "circlize", "tikzDevice"))

Rfast and Matrix are used by the simulation and real data scripts; ggplot2, latex2exp and circlize by the figure scripts.

For a comparison of the SPIEC-EASI method, please install

library(devtools)
install_github("zdk123/SpiecEasi")
library(SpiecEasi)

For a comparison of the metagenomeSeq method, please install

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

Then, we are able to install metagenomeSeq package and ComplexHeatmap for plotting:

BiocManager::install("ComplexHeatmap")
BiocManager::install("metagenomeSeq")
BiocManager::install("MetaLonDA")

MetaLonDA::normalize performs the cumulative sum scaling (CSS) normalization used in Fig 1.

Please make sure the C++ compiler is correctly installed. For Mac user, please install Xcode from command line tools. Execute the command xcode-select --install on Terminal.

R project reproducing instruction

An user friendly way to reproduce the results of the main text and the algorithm, one can create a R project by following the steps: Open Rstudio - File - New Project - Version Control - Git - enter the github link: https://github.com/shuang-jie/SP-BGFM.

Then, All simulation codes and real data analysis codes are in each sub-folder. All figures in the maintext can be reproduced by the figures-code.

Organization

R, src, man, data, vignettes, tests

The R package.

simulation-code

Running Sim 1-5.R produces the results displayed in Sim 1-5 in the paper and please save the result as RData with this.sim.id.

real-data

Filtered7539OTUs.RData is the multi-domain skin microbiome data from the paper. The same data is available in the package as data(skin).

In the real data, it contains:

Y1 : bacterial microbiome count table. 60 samples $\times$ 75 OTUs. Each row is a sample, and each column is a bacterial OTU.

Y2 : viral microbiome count table. 60 samples $\times$ 39 OTUs. Each row is a sample, and each column is a viral OTU.

Y : combined multi-domain skin microbiome data. 60 samples $\times$ 114(75+39) OTUs.

X : a categorical covariate representing experimental conditions. (1,0,0) pre-treatment & (0,1,0) post-treatment & (0,0,1) healthy condition. That is, $\beta_{mj1}$ , $\beta_{mj2}$ and $\beta_{mj3}$ correspond to the pre-treatment, post-treatment and healthy conditions, respectively.

J : number of OTUs in each domain. (75, 39)

Jsum : number of total OTUs. 114

n : number of samples. 60

S : number of subjects. 20

real-data-code

Real Data.R produces the results displayed in real data in the paper and please save the result as RData with Real Data.RData

figures-code

Produces Figures 1-10 in the main text.

figures

Contains the results from calling the code in the figures-code folder.