为什么要对数据做批次矫正?
对数据进行批次矫正是数据预处理中常见的步骤,主要是为了消除在实验过程中因批次差异而导致的系统性误差。这种差异可能来源于多个因素,如实验条件、样本处理、设备校准、技术人员之间的差异等。批次效应可能在数据分析时掩盖实际的生物学信号,从而影响结果的可靠性和可重复性。
常见的批次矫正
1.sva包中的Combat
2.limma包中的removeBatchEffect
两种的区别
- ComBat 更加强调在多批次和复杂数据集的情况下消除批次效应,且能提供更强的统计模型基础。并且更加适合多源数据集
- removeBatchEffect 则更适用于简单场景,快速有效,但不那么灵活。
Combat
下载SVA包并加载
install.packages("BiocManager") # 如果尚未安装 BiocManager
BiocManager::install("sva")
library(sva)
加载文件并提取信息
# 读取数据
data <- read.csv("C:\\Users\\Deskto\\AAA.csv")
# 提取ID列和sample_group列
sample_id <- data[, 1] # ID列
sample_group <- data[, ncol(data) - 1] # sample_group列(倒数第二列)
batch <- data[, ncol(data)] # batch列(倒数第一列)
准备两个矩阵——基因表达量矩阵,设计矩阵
# 提取基因表达数据
#基因表达量矩阵中只能有基因的表达数据,不能存在字符型所以需要去掉第一列的样本ID,倒数一二列的分组和批次信息
gene_expression <- data[, 2:(ncol(data) - 2)] # 基因表达量数据
#将 sample_group 变量转换为因子类型(factor)
sample_group <- factor(sample_group)
# 使用model.matrix生成设计矩阵
design_matrix <- model.matrix(~sample_group)
#将 batch 变量转换为因子类型(factor)
batch <- factor(batch)
这里设计矩阵中包含分组信息是为了在进行批次矫正时将分组信息作为协变量一些考虑,这样会更加准确,一般来说批次矫正时需要考虑分组信息
进行Combat并且保存数据
# 使用Combat进行批次效应消除
#这里加一个t()是为了将基因表达矩阵转置,因为我的数据是基因ID为列名,样本ID为行名
combat_data <- ComBat(dat=t(as.matrix(gene_expression)), batch=batch, mod=design_matrix, par.prior=TRUE, prior.plots=FALSE)
# 将结果合并回原始数据
combat_datat <- data.frame(ID = sample_id, sample_group = sample_group, batch = batch, t(combat_data))
# 保存结果
write.csv(combat_datat, "C:\\Users\\Desktop\\new\\AAA_combat.csv", row.names = FALSE)
removeBatchEffect
下载并加载limma包
install.packages("BiocManager") # 如果尚未安装 BiocManager
BiocManager::install("limma")
library(limma)
expression_data_t<-t(expression_data )
# 注意:这里不再转置expression_data,因为removeBatchEffect通常直接应用于表达矩阵
# 如果你的数据有特殊需求需要转置,请确保在转置前后理解数据的结构
corrected_data <- removeBatchEffect(expression_data_t, batch = batch)
corrected_data<-t(corrected_data)
其余大多代码和combat一致