昨天收到一个读者的消息，聊到在处理tpm数据时为什么需要去除批次效应，早上5点多就睡不着，就这个读者问题来跟大家聊一下为什么做GEO数据挖掘前需要去除批次效应，以及如何选择不同的去除批次效应方法。

为什么需要去除批次效应？

在做基因表达分析时，我们会遇到一个常见的问题——批次效应。简单来说，批次效应是指由于实验条件、时间、操作人员、平台等因素的不同，导致不同批次或平台间的样本表现出系统性的差异。这些差异有时候并不是生物学上真实存在的，而是实验过程中产生的伪差异。批次效应如果不处理好，可能会干扰我们的分析结果，甚至误导我们的结论。

举个例子，假设我们有来自不同实验或平台的数据，如果这些数据没有经过批次效应的处理，就可能看到样本之间存在非生物学上的偏差。比如，某些样本的基因表达水平看起来特别高或特别低，但其实这些差异是由实验条件引起的，而不是由于样本的生物学特征（如疾病状态、分组等）造成的。

所以，去除批次效应是非常重要的，尤其是在我们把多个数据集合并或者需要比较不同组别时。常见的去批次效应的方法包括数据的标准化、归一化，或者使用一些专门的统计方法，比如我们常见limma包中的normalizeBetweenArrays和 removeBatchEffect 函数。

接下来我们来看下什么是用normalizeBetweenArrays去批次效应，什么时候用removeBatchEffect 函数。

先说结论：单个数据集使用normalizeBetweenArrays，多个数据集使用removeBatchEffect

下面我们用GSE89143和GSE83521这两个数据集，配合箱线图，来给大家展示去批次效应的步骤和去批次效应前后的箱线图对比。

单个数据集使用normalizeBetweenArrays

我们通过 GEOquery 包下载数据集GSE89143 。然后从中提取了表达矩阵，再绘制未去批次效应前的箱线图

library
(GEOquery)

library
(limma)


# 下载并加载数据集 GSE89143

GSE89143 <- getGEO(
"GSE89143"
, destdir = 
"."
, getGPL = 
FALSE
)

exp1 <- exprs(GSE89143[[
1
]])


range(exp1) 
# 85 65535


exp1 <- log2(exp1 + 
1
) 
# 对数据进行log2转换


range(exp1) 
# 6.426265 16.000000

# 未去批次效应前绘制的boxplot

boxplot(exp1, las = 
2
)

我们可以看到，上面是还没去批次效应前的箱线图，数值的表达量并不在同一条水平线上，接着我们使用normalizeBetweenArrays 对单个数据集进行去批次效应

exp1 <- normalizeBetweenArrays(exp1)

range(exp1) 
# 6.626236 15.445942

# 已经去批次效应前绘制的boxplot

boxplot(exp1, las = 
2
) 

可以看到去批次效应后，箱线图的数值的表达量已经在同一条水平线上了

多个数据集使用removeBatchEffect

接下来，我们重复同样的步骤下载并处理GSE83521数据集

GSE83521 <- getGEO(
"GSE83521"
, destdir = 
"."
, getGPL = 
FALSE
)

exp2 <- exprs(GSE83521[[
1
]])


range(exp2)

# 绘制 boxplot 以观察数据分布

boxplot(exp2, las = 
2
)

exp2 <- normalizeBetweenArrays(exp2)

boxplot(exp2, las = 
2
)

再合并两个数据集的数据，并且绘制合并后的箱线图

# 合并GSE89143 和 GSE83521的表达数据

common_genes <- intersect(rownames(exp1), rownames(exp2))

exp_merge <- cbind(exp1[common_genes, ], exp2[common_genes, ])


# 绘制合并后的boxplot

boxplot(exp_merge)

可以看到，上面是合并后的数据集 exp_merge 显示出明显的批次效应，两个数据集的样本仍然呈现出不同的分布。因此，我们需要进一步去除批次效应。

为了消除合并后的数据集中的批次效应，我们使用了 removeBatchEffect 函数，来根据样本的批次信息进行批次效应的调整。

# 定义批次信息

batch <- c(rep(
'GSE89143'
, 
6
), rep(
'GSE83521'
, 
12
))


# 定义分组信息

group_list <- c(rep(c(
'tumor'
, 
'normal'
), each = 
3
), rep(
'tumor'
, 
6
), rep(
'normal'
, 
6
))


# 创建设计矩阵

design <- model.matrix(~ group_list)


# 多个数据集合并使用limma的removeBatchEffect去除批次效应

exp_merge <- removeBatchEffect(exp_merge, batch = batch, design = design)


range(exp_merge)

# 绘制去除批次效应后的boxplot

boxplot(exp_merge)

从上图可以看到，去除了批次效应后，两个数据集的表达量已经在同一条水平线上了。

到这一步，数据清洗已经完成了，接下来就可以根据自个的需求，做后续的差异表达分析、富集分析等。

写在最后

总结就一句话，去批次效应时，单个数据集使用normalizeBetweenArrays，多个数据集使用removeBatchEffect

如果大家不想写代码，可以用这个零代码在线生信分析平台，支持单细胞转录组、TCGA转录组、通用Bulk转录组分析以及多种绘图功能。立即体验，无需R环境配置。