R语言实现去除批次效应

在整理GEO数据库中的数据时，我们常遇到来源于不同的GSE数据集，为增大样本量，常需要合并数据集，但是由于种种原因：不同技术人员、试剂、仪器、时间等导致，导致相同研究的结果存在数据偏差，此时需要进行校正。

基本方法：

样本量较少

#校正批次效应
library(limma)
y <- matrix(rnorm(20*9),20,9)
y[,1:3] <- y[,1:3] + 5
y[,c(3,6,9)]<y[,c(3,6,9)]+100
batch <- c("A","A","A","B","B","B","C","C","C")
#加入design
grouplist<-rep(c("tum","unkown","nor"),3)
g=factor(grouplist)
g
design1=model.matrix(~g)
#使用removeBatchEffect
y2 <- removeBatchEffect(y, batch)
y3<-removeBatchEffect(y, batch=batch,design=design1)
par(mfrow=c(1,3))
boxplot(as.data.frame(y),main="Original")
boxplot(as.data.frame(y2),main="Batch corrected")
boxplot(as.data.frame(y3),main="Batch-design corrected")
#使用ComBat
library(sva)
C1<-ComBat(dat = y,batch=batch,mod=NULL)
C2<-ComBat(dat = y,batch=batch,mod=design1)
par(mfrow=c(1,2))
boxplot(as.data.frame(C1),main="Batch corrected")
boxplot(as.data.frame(C2),main="Batch-design corrected")
#对比removeBatchEffect和ComBat
library(vioplot)
par(mfrow=c(1,2))
vioplot(as.data.frame(y2),main="removeBatchEffect corrected")
vioplot(as.data.frame(C2),main="ComBat corrected")
#相较而言，效果差不多，ComBat对处理不同批次下的分组效果更加明显

 

 

 样本量较多

#样本数量多
library(limma)
y <- matrix(rnorm(20*90),20,90)
y[,1:30] <- y[,1:30] + 5
y[,31:60] <- y[,31:60] + 50
y[,c(1:15,31:45,61:75)]<y[,c(1:15,31:45,61:75)]+100
batch <- c(rep("A",30),rep("B",30),rep("C",30))
#加入design
grouplist<-rep(c("tum","nor"),45)
g=factor(grouplist)
g
design1=model.matrix(~g)
#使用removeBatchEffect