R语言学习笔记-ComBat函数的使用

ComBat 是 R 中用于校正批次效应（batch effect）的强大工具，尤其在组学数据（如代谢组学、转录组学）分析中广泛应用。以下是关于其作用和用法的详细说明：

​​一、ComBat 的核心作用​​

解决什么问题​​
消除因实验批次（不同时间、不同操作员、不同仪器等）引入的技术变异，保留真实的生物学差异。

算法原理​​

基于经验贝叶斯方法（Empirical Bayes），对每个特征（如代谢物）进行：

1. ​​均值校正​​：调整不同批次间的基线偏移
2. 方差校正​​：缩小不同批次的离散度差异

​​适用场景

1. 多批次实验数据整合
2. 纵向研究中不同时间点的数据对齐
3. 跨平台数据标准化

​​二、基本用法（以代谢组学数据为例）​​

​​1. 输入数据要求​​

• ​​dat​​：数值矩阵，​​行是特征​​（如代谢物），​​列是样本​​
• ​​batch​​：向量，指定每个样本所属的批次（必须为因子或字符型）
• ​​mod​​（可选）：设计矩阵，指定需要保留的生物学协变量（如疾病分组）

​​2. 基础代码示例​​

library(sva)



# 假设数据已准备：

# metab_data: 代谢物矩阵（行=代谢物，列=样本）

# batch: 批次信息（如 c("Batch1","Batch1","Batch2","Batch2")）

# group: 生物学分组（如 c("Case","Control","Case","Control")）



# 无协变量校正（仅去批次效应）

corrected_data <- ComBat(

  dat = metab_data,

  batch = batch,

  par.prior = TRUE  # 使用参数先验（默认）

)



# 包含协变量校正（保留组间差异）

design <- model.matrix(~ group)  # 设计矩阵

corrected_data <- ComBat(

  dat = metab_data,

  batch = batch,

  mod = design,     # 指定需保留的生物学变量

  par.prior = TRUE

)

​​三、关键参数详解​​

​​四、注意事项​​

​​避免协变量与批次共线性​​

错误示例：若所有Case样本都在Batch1，则group与batch完全共线，会导致报错：

Error: At least one covariate is confounded with batch!

解决方案：检查交叉表 table(batch, group)，确保各分组在所有批次中均有分布。

数据预处理​​

• 建议先进行对数转换（如 log2(x+1)）使数据更接近正态分布
• 去除缺失值过多的特征（如 >20% 缺失）

结果验证​​

• PCA图​​：校正后批次间样本应混合分布
• ​​箱线图​​：查看批次间分布是否对齐

五、完整实战案例​​

# PCA验证示例

pca <- prcomp(t(corrected_data), scale. = TRUE)

plot(pca$x[,1:2], col = batch, pch = 16)

# 加载数据

data <- read.csv("metabolites.csv", row.names = 1)  # 行=代谢物，列=样本

batch <- c(rep("Batch1", 10), rep("Batch2", 10))    # 20个样本的批次

group <- rep(c("Case","Control"), 10)               # 病例/对照分组



# 1. 缺失值处理

data <- data[rowSums(is.na(data)) < 5, ]  # 保留缺失<5个样本的代谢物



# 2. 对数转换

data_log <- log2(data + 1)



# 3. 检查批次-分组共线性

table(batch, group)  # 确保各分组在批次中分布均衡



# 4. ComBat校正

design <- model.matrix(~ group)  # 保留组间差异

corrected <- ComBat(

  dat = as.matrix(data_log),

  batch = batch,

  mod = design,

  par.prior = TRUE

)



# 5. 验证效果

library(ggplot2)

pca <- prcomp(t(corrected), scale. = TRUE)

pca_df <- data.frame(PC1 = pca$x[,1], PC2 = pca$x[,2], Batch = batch, Group = group)

ggplot(pca_df, aes(PC1, PC2, color = Batch, shape = Group)) +

  geom_point(size = 3) +

  stat_ellipse(aes(group = Batch)) +

  theme_bw()

​​六、常见问题解答​​

​​Q1：如何选择是否使用mod参数？​​

• 若已知某些协变量（如年龄、性别）影响表达，需在mod中指定以保留这些效应
• 若仅需去除批次效应，可省略mod

​​Q2：校正后数据出现负值怎么办？​​

• ComBat 校正可能产生负值（尤其对原始计数数据）
• 解决方法：校正前进行对数转换，或对结果取exp()转换（视分析需求而定）

​​Q3：与limma::removeBatchEffect的区别？​​

• ComBat：基于经验贝叶斯，更适合小样本数据
• removeBatchEffect：基于线性模型，更适合大样本且批次效应较弱的情况