662 微生物组可重复性分析¶

# R语言：检测批次效应
library(vegan)     # 群落生态学
library(ggplot2)   # 绑图

# 读取数据
otu <- read.csv("otu_table.csv", row.names = 1)  # OTU表
meta <- read.csv("metadata.csv")  # 元数据，含batch列

# 1. PCoA可视化看批次分离
bc_dist <- vegdist(t(otu), method = "bray")  # Bray-Curtis距离
pcoa <- cmdscale(bc_dist, k = 2, eig = TRUE)  # PCoA降维

df_pcoa <- data.frame(
  PC1 = pcoa$points[,1],  # 第一主坐标
  PC2 = pcoa$points[,2],  # 第二主坐标
  Batch = meta$batch,     # 批次
  Group = meta$group      # 实验分组
)

# 画PCoA图
ggplot(df_pcoa, aes(x = PC1, y = PC2, color = Batch, shape = Group)) +
  geom_point(size = 3) +
  theme_minimal() +
  ggtitle("校正前：批次效应检测")

# 2. PERMANOVA量化批次效应
perm_batch <- adonis2(bc_dist ~ batch + group,
                       data = meta,
                       permutations = 999)  # PERMANOVA检验
print(perm_batch)  # 看batch的R²值
# R²越大，批次效应越严重

# 3. PERMDISP检验组内离散度
bd <- betadisper(bc_dist, meta$batch)  # 批次间离散度
permutest(bd)  # 置换检验

# ComBat批次效应校正
library(sva)  # 包含ComBat函数

# 准备数据
otu_matrix <- as.matrix(otu)  # OTU矩阵
batch <- meta$batch            # 批次信息
mod <- model.matrix(~ group, data = meta)  # 保护变量（实验分组）

# CLR转换（ComBat需要近似正态分布的数据）
library(compositions)  # 组成数据分析
otu_clr <- clr(t(otu_matrix) + 0.5)  # CLR转换，加伪计数

# 应用ComBat
combat_data <- ComBat(
  dat = t(otu_clr),    # 数据矩阵（特征×样本）
  batch = batch,        # 批次向量
  mod = mod             # 保护变量设计矩阵
)
# 注意：ComBat假设高斯分布，对微生物组数据可能不理想

# MMUPHin批次效应校正
# BiocManager::install("MMUPHin")
library(MMUPHin)  # 微生物组专用批次校正

# 准备数据
otu_matrix <- as.matrix(otu)  # OTU丰度矩阵

# 批次校正
adjusted <- adjust_batch(
  feature_abd = otu_matrix,  # 特征丰度矩阵（OTU×样本）
  batch = "batch",           # 批次列名
  covariates = "group",      # 需要保护的协变量
  data = meta                # 元数据
)

# 校正后的丰度矩阵
otu_corrected <- adjusted$feature_abd_adj  # 校正后的数据

# 验证校正效果
bc_adj <- vegdist(t(otu_corrected), method = "bray")  # 校正后距离
perm_adj <- adonis2(bc_adj ~ batch + group,
                     data = meta,
                     permutations = 999)
print(perm_adj)  # 校正后batch的R²应该降低

# ConQuR批次效应校正
# install.packages("ConQuR")
library(ConQuR)  # 条件分位数回归

# 应用ConQuR
conqur_result <- ConQuR(
  tax_tab = otu_matrix,      # OTU丰度表（样本×OTU）
  batchid = meta$batch,      # 批次ID
  covariates = meta[, "group", drop = FALSE],  # 协变量
  batch_ref = "Batch1"       # 参考批次
)

# 校正后的计数矩阵
otu_conqur <- conqur_result  # ConQuR输出的是校正后的计数

# ConQuR优势：输出的是零膨胀计数数据，保持了微生物组数据的特性

# 综合评估批次校正效果
library(patchwork)  # 组合多个ggplot

# 1. 计算校正前后的R²
methods <- list(
  "原始" = otu,
  "ComBat" = combat_data,
  "MMUPHin" = otu_corrected,
  "ConQuR" = otu_conqur
)

results <- data.frame()
for (name in names(methods)) {
  d <- vegdist(t(methods[[name]]), method = "bray")  # 计算距离
  perm <- adonis2(d ~ batch + group, data = meta, permutations = 999)
  results <- rbind(results, data.frame(
    Method = name,
    Batch_R2 = perm["batch", "R2"],      # 批次解释的变异
    Group_R2 = perm["group", "R2"],      # 分组解释的变异
    Batch_p = perm["batch", "Pr(>F)"]    # 批次p值
  ))
}
print(results)
# 好的校正：Batch_R2降低，Group_R2保持或增加

# 2. 比较不同距离度量下的效果
distances <- c("bray", "jaccard", "canberra", "manhattan")
for (dist_method in distances) {
  d <- vegdist(t(otu_corrected), method = dist_method)
  perm <- adonis2(d ~ batch, data = meta, permutations = 999)
  cat(sprintf("%s距离 - 批次R²: %.4f, p: %.4f\n",
              dist_method, perm["batch", "R2"], perm["batch", "Pr(>F)"]))
}

# 跨研究整合分析
library(MMUPHin)

# 1. 合并多个研究的数据
study1_otu <- read.csv("study1_otu.csv", row.names = 1)
study2_otu <- read.csv("study2_otu.csv", row.names = 1)

# 确保OTU名称一致
common_otus <- intersect(rownames(study1_otu), rownames(study2_otu))
merged_otu <- cbind(study1_otu[common_otus,], study2_otu[common_otus,])

# 2. 批次校正
merged_meta <- rbind(meta1, meta2)  # 合并元数据
merged_meta$study <- c(rep("Study1", ncol(study1_otu)),
                       rep("Study2", ncol(study2_otu)))

adjusted_merged <- adjust_batch(
  feature_abd = as.matrix(merged_otu),
  batch = "study",
  covariates = "disease_status",
  data = merged_meta
)

# 3. Meta分析差异丰度
meta_result <- lm_meta(
  feature_abd = adjusted_merged$feature_abd_adj,
  batch = "study",
  exposure = "disease_status",
  data = merged_meta
)
print(head(meta_result))  # 差异丰度物种

# 批次效应处理流程
1. 检测：PCoA可视化 + PERMANOVA量化
2. 选择方法：
   - 通用数据 → ComBat
   - 微生物组（零膨胀） → MMUPHin
   - 非参数偏好 → ConQuR
   - 2025最新 → 复合分位数回归
3. 评估：比较校正前后batch R²
4. 验证：确保生物学信号未被消除

# 方法对比
ComBat: 快速，假设高斯，微生物组不理想
MMUPHin: 处理零膨胀，假设零膨胀高斯
ConQuR: 非参数，保持零膨胀计数
2025新方法: 基于负二项，Bray-Curtis/Manhattan最优

# 报告标准
STORMS检查清单：微生物组研究报告标准
Human Microbiome Action：欧洲微生物组研究协调