来自不同平台GEO数据批次效应除去

批次效应除去-GEO案例分析

最近分析来自不同平台的GEO数据，GEO分析推荐学习Jimmy大神的B站视频，也可以学习我之前的教程。limma和sva两种方法进行批次效应去除。

1、加载软件和数据

rm(list = ls())

options(stringsAsFactors = F)

library(sva)

library(limma)

load("merge.Rdata")#加载上一步中合并不同平台的数据

处理前

图中存在明显的批次效应，两种方法除去批次效应，limma和sva两种方法进行批次效应去除。

2、limma软件包中的removeBatchEffect

batch1 <- c(rep('GSE*',30),rep('GSE*',15))

batch1 <- as.factor(batch1)

design <- model.matrix(~0 + batch1)
#校正其实就一步

ex_b_limma <- removeBatchEffect(x_merge1,
                                batch = batch1)
boxplot(ex_b_limma,las = 2)
save(ex_b_limma,x_merge,x_merge1,batch1,file = "ex_b_limma.Rdata")

removeBatchEffect批次除去后

构建批次矩阵， removeBatchEffect运行后表达水平基本在同一个水平上，可以进行下游差异分析等。

3、sva软件包ComBat去除批次效应

library(sva)

batch1 <- c(rep('GSE*',30),rep('GSE6*',15))

x_merge1 <- as.data.frame(x_merge)
class(x_merge1)
x_merge2 <- as.matrix(x_merge1)#关键的一步，转换为有向量、数值或者矩阵

design <- model.matrix(~0 + batch1)

#3.设置model（可选）
#mod = model.matrix(~as.factor(batch1), data=x_merge1)


#4.校正其实就一步
combat_edata <- ComBat(dat = x_merge2, batch = batch1)

dim(combat_edata)
boxplot(combat_edata)

save(batch1,x_merge2,design,file = "Combat_edata.Rdata")

Combat去除批次效应

sva去除批次效应后表达量基本在一个水平线上，可以进行下游差异分析等。
可以看一下批次效应消除前后，样本聚类的情况，还可以比较一下，两种批次效应处理方法差异基因的异同。

参考

不同矫正批次效应方法的比较 - 生信技能树

GEO 批次效应就靠一个函数搞定