參考教程：
微信公眾號：生信星球
批次效應(yīng)處理實例：combat和removebatcheffect的對比
特別感謝：
人美心善愛護(hù)小白的花花老師
小潔忘了怎么分身
作為一個非常想搞事情的生信小白，我一直想做不同數(shù)據(jù)集合并分析盼玄，下面我就把我的兩個數(shù)據(jù)集合并的心路歷程給大家還原一下
希望大家提出意見或者問題，一起進(jìn)步~
1. 批次效應(yīng)(Batch Effects)
可以理解為：樣本受到檢測的實驗室條件强岸、試劑批次和人員差異的影響锻弓，對結(jié)果的準(zhǔn)確性造成了影響
2.多組數(shù)據(jù)集合并分析的流程
（1）條件
- 取材對象應(yīng)為同一組織
- 本方法適用于同芯片平臺
（2）流程
- 了解數(shù)據(jù)蝌箍，比如判斷原數(shù)據(jù)集是否為標(biāo)準(zhǔn)化數(shù)據(jù)
（在進(jìn)行操作之前，總得弄清數(shù)據(jù)是個什么樣子的吧） - 提取老三套：表達(dá)矩陣妓盲、臨床表型（也就是你的分組信息）、平臺數(shù)據(jù)
- 各數(shù)據(jù)集組間校正
（有時候雖然作者對數(shù)據(jù)進(jìn)行了標(biāo)準(zhǔn)化悯衬，但是不是我們想要的弹沽，還得需要再次normalization一下） - 合并矩陣，批次矯正
兩種方法：
limma removeBatchEffect
sva combat - 再次查看合并且去除批次差異的數(shù)據(jù)
（沒有前后對比怎么能確定批次差異矯正有效果呢）
3.我的流程
首先炸渡，了解數(shù)據(jù)
我想合并的兩個數(shù)據(jù)集，取材相同丽已，同一芯片平臺
（GEO網(wǎng)站上都能找得到）
然后，讀入數(shù)據(jù)沛婴，提取數(shù)據(jù)
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
library(stringr)
gse = "GSE79704"
eSet1 <- getGEO("GSE79704",
destdir = '.',
getGPL = F)
eSet2 <- getGEO("GSE83582",
destdir = '.',
getGPL = F)
#(1)提取表達(dá)矩陣exp
exp1 <- exprs(eSet1[[1]])
exp1[1:4,1:4]
exp2 <- exprs(eSet2[[1]])
exp2[1:4,1:4]
#將兩個數(shù)據(jù)集行數(shù)調(diào)整一致,不然沒辦法合并
table(rownames(exp1) %in% rownames(exp2))
#TRUE
#25560
length(intersect(rownames(exp1),rownames(exp2)))
#[1] 25560
exp1 <- exp1[intersect(rownames(exp1),rownames(exp2)),]
exp2 <- exp2[intersect(rownames(exp1),rownames(exp2)),]
#(2)提取臨床信息
pd1 <- pData(eSet1[[1]])
pd2 <- pData(eSet2[[1]])
#(3)提取芯片平臺編號
gpl1 <- eSet1[[1]]@annotation
gpl2 <- eSet2[[1]]@annotation
#兩個gpl均為GPL19983
兩者的boxplot如下：
exp1.png
exp2.png
然后我仔細(xì)思考了要不要先進(jìn)行組間校正再合并
所以如下我進(jìn)行了兩條路線
（1）直接批次矯正
【1】合并數(shù)據(jù)
#開始合并
exp = cbind(exp1,exp2)
#合并group_list
group_list1 = c(rep('NN',20),rep('PP',12),rep('GPP',32))
group_list2 = c(rep('NN',12),rep('PP',12),rep('NN',8),rep('EP',30),rep('IP',40))
group_list = c(group_list1,group_list2)
table(group_list)
#group_list
#EP GPP IP NN PP
#30 32 40 40 24
boxplot如下：
#boxplot
boxplot(exp,las=2,main='exp-before')
exp-before.png
cluster如下：
#有沒有批次效應(yīng)，看一下聚類的情況
dist_mat <- dist(t(exp))
clustering <- hclust(dist_mat, method = "complete")
plot(clustering,labels = group_list)
cluster-exp-before.png
PCA如下：
#PCA
dat=as.data.frame(t(exp))
library(FactoMineR)#畫主成分分析圖需要加載這兩個包
library(factoextra)
# pca的統(tǒng)一操作走起
dat.pca <- PCA(dat, graph = FALSE)
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
#palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
PCA-exp-before.png
【2-1】limma方法
#處理批次效應(yīng):limma
library(limma)
#?removeBatchEffect()
batch <- c(rep("A",64),rep("B",102))
y2 <- removeBatchEffect(exp, batch)
boxplot如下：
exp-limma.png
cluster如下：
cluster-exp-limma.png
PCA如下：
pca-exp-limma.png
【2-2】combat方法
#處理批次效應(yīng)(combat)
library(sva)
?ComBat
batch <- c(rep("A",64),rep("B",102))
mod = model.matrix(~as.factor(group_list))
y = ComBat(dat=exp, batch=batch, mod=mod)
#Found2batches
#Adjusting for4covariate(s) or covariate level(s)
#Standardizing Data across genes
#Fitting L/S model and finding priors
#Finding parametric adjustments
#Adjusting the Data
boxplot如下：
exp-combat.png
cluster如下：
cluster-exp-combat.png
PCA如下：
pca-exp-combat.png
（2）先組間校正再批次校正
【1】組間校正
#歸一化
exp1 = limma::normalizeBetweenArrays(exp1)
exp2 = limma::normalizeBetweenArrays(exp2)
boxplot(exp1,las=2,main='exp1-normalization')
boxplot(exp2,las=2,main='exp2-normalization')
boxplot如下：
exp1-normalization.png
exp2-normalization.png
【2】合并數(shù)據(jù)
boxplot如下：
exp-n-before.png
cluster如下：
cluster-exp-n-before.png
PCA如下：
PCA-exp-n-before.png
【3-1】limma方法
boxplot如下：
exp-n-limma.png
cluster如下：
cluster-exp-n-limma.png
PCA如下：
pca-exp-n-limma.png
【3-2】combat方法
boxplot如下：
exp-n-combat.png
cluster如下：
cluster-exp-n-combat.png
PCA如下：
pca-exp-n-combat.png
