371_微生物组多中心研究设计¶

# 安装必要包
if (!require("MMUPHin")) BiocManager::install("MMUPHin")
library(MMUPHin)   # 多中心微生物组元分析框架
library(vegan)     # 用于 PERMANOVA 检验
library(ggplot2)   # 可视化

# 模拟加载相对丰度矩阵（样本×物种）
# 行=物种，列=样本（MMUPHin要求格式）
abundance <- read.csv("otu_table.csv", row.names = 1)
metadata  <- read.csv("metadata.csv", row.names = 1)

# 步骤1：检验批次效应是否显著
# PERMANOVA检验：看"中心"能解释多少组间差异
dist_matrix <- vegdist(t(abundance), method = "bray")  # Bray-Curtis距离
permanova_result <- adonis2(
  dist_matrix ~ Center,      # 按"中心"分组检验
  data = metadata,
  permutations = 999         # 置换检验次数
)
print(permanova_result)      # R² 越大，批次效应越强

# 步骤2：用 MMUPHin 校正批次效应
fit_adjust <- adjust_batch(
  feature_abd = abundance,    # 相对丰度矩阵（物种×样本）
  batch = "Center",           # 批次变量列名（对应metadata列）
  covariates = "Disease",     # 需要保留的生物学变量
  data = metadata,            # 元数据
  control = list(verbose = TRUE)
)
abundance_corrected <- fit_adjust$feature_abd_adj  # 取校正后丰度

# 步骤3：荟萃分析差异物种
lm_meta <- lm_meta(
  feature_abd = abundance_corrected,
  batch = "Center",
  covariates = "Disease",
  data = metadata,
  control = list(verbose = TRUE)
)
# 查看显著差异物种
sig_taxa <- lm_meta$meta_fits[lm_meta$meta_fits$qval.fdr < 0.05, ]
print(sig_taxa)

library(lme4)      # 混合效应模型
library(lmerTest)  # 提供 p 值

# 以 Shannon 多样性为例，检验疾病效应
metadata$shannon <- vegan::diversity(t(abundance), index = "shannon")

# 混合模型：疾病为固定效应，中心为随机效应
model <- lmer(
  shannon ~ Disease + Age + BMI + (1 | Center),  # (1|Center)表示中心随机截距
  data = metadata
)
summary(model)   # 查看疾病效应的 p 值