批次效應(yīng)除去-GEO案例分析
最近分析來(lái)自不同平臺(tái)的GEO數(shù)據(jù)，GEO分析推薦學(xué)習(xí)Jimmy大神的B站視頻，也可以學(xué)習(xí)我之前的教程。limma和sva兩種方法進(jìn)行批次效應(yīng)去除。
1、加載軟件和數(shù)據(jù)
rm(list = ls())
options(stringsAsFactors = F)
library(sva)
library(limma)
load("merge.Rdata")#加載上一步中合并不同平臺(tái)的數(shù)據(jù)
處理前
圖中存在明顯的批次效應(yīng)，兩種方法除去批次效應(yīng)，limma和sva兩種方法進(jìn)行批次效應(yīng)去除。
2、limma軟件包中的removeBatchEffect
batch1 <- c(rep('GSE*',30),rep('GSE*',15))
batch1 <- as.factor(batch1)
design <- model.matrix(~0 + batch1)
#校正其實(shí)就一步
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批次除去后
構(gòu)建批次矩陣， removeBatchEffect運(yùn)行后表達(dá)水平基本在同一個(gè)水平上，可以進(jìn)行下游差異分析等。
3、sva軟件包ComBat去除批次效應(yīng)
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)#關(guān)鍵的一步，轉(zhuǎn)換為有向量、數(shù)值或者矩陣
design <- model.matrix(~0 + batch1)
#3.設(shè)置model（可選）
#mod = model.matrix(~as.factor(batch1), data=x_merge1)
#4.校正其實(shí)就一步
combat_edata <- ComBat(dat = x_merge2, batch = batch1)
dim(combat_edata)
boxplot(combat_edata)
save(batch1,x_merge2,design,file = "Combat_edata.Rdata")
Combat去除批次效應(yīng)
sva去除批次效應(yīng)后表達(dá)量基本在一個(gè)水平線上，可以進(jìn)行下游差異分析等。
可以看一下批次效應(yīng)消除前后，樣本聚類的情況，還可以比較一下，兩種批次效應(yīng)處理方法差異基因的異同。
