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