参考教程：
微信公众号：生信星球
批次效应处理实例：combat和removebatcheffect的对比
特别感谢：
人美心善爱护小白的花花老师
小洁忘了怎么分身
作为一个非常想搞事情的生信小白，我一直想做不同数据集合并分析，下面我就把我的两个数据集合并的心路历程给大家还原一下
希望大家提出意见或者问题，一起进步~
1. 批次效应(Batch Effects)
可以理解为：样本受到检测的实验室条件、试剂批次和人员差异的影响，对结果的准确性造成了影响
2.多组数据集合并分析的流程
（1）条件
- 取材对象应为同一组织
- 本方法适用于同芯片平台
（2）流程
- 了解数据，比如判断原数据集是否为标准化数据
（在进行操作之前，总得弄清数据是个什么样子的吧） - 提取老三套：表达矩阵、临床表型（也就是你的分组信息）、平台数据
- 各数据集组间校正
（有时候虽然作者对数据进行了标准化，但是不是我们想要的，还得需要再次normalization一下） - 合并矩阵，批次矫正
两种方法：
limma removeBatchEffect
sva combat - 再次查看合并且去除批次差异的数据
（没有前后对比怎么能确定批次差异矫正有效果呢）
3.我的流程
首先，了解数据
我想合并的两个数据集，取材相同，同一芯片平台
（GEO网站上都能找得到）
然后，读入数据，提取数据
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)提取表达矩阵exp
exp1 <- exprs(eSet1[[1]])
exp1[1:4,1:4]
exp2 <- exprs(eSet2[[1]])
exp2[1:4,1:4]
#将两个数据集行数调整一致,不然没办法合并
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
然后我仔细思考了要不要先进行组间校正再合并
所以如下我进行了两条路线
（1）直接批次矫正
【1】合并数据
#开始合并
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如下：
#有没有批次效应，看一下聚类的情况
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的统一操作走起
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方法
#处理批次效应: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方法
#处理批次效应(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】合并数据
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
