样品不够？结果太low？多芯片联合分析，拯救你的困境

一文学会GEO多芯片批次校正方法

大家好，我是阿琛。在上一期内容中，我们学习了如何用简单的方法，快速完成GEO芯片分析（一文解决GEO芯片数据分析80%的工作！建议收藏！）。解决了主要矛盾以后，我们需要看看另一个限制大家发文章的问题，即样本数量问题。在第一期打赏营中由于小伙伴提出，现场发挥简单讲了讲这个内容，今天来给大家详细的聊聊。

对于GEO数据集而言，尤其是非肿瘤研究的小伙伴们，样本数是一个主要的限制问题。与土豪级别的TCGA数据库不同，很多GEO数据集的样本量往往较少。这样一来，无论是差异分析，还是其他相关性分析，如WGCNA分析要求的至少15个样本数量，由于数量不足导致分析结果可信度不高，或者结果图看起来远远没有那么的高大上。

此时，我们需要把多个数据集进行整合分析，综合其中的样品进行后续分析。然而，即便是相同的实验内容，在不同的时间点操作，或者不同熟练程度的实验员进行，甚至是实验试剂的不同，都会导致分析结果出现批次效应。因此，在GEO多芯片合并时，我们需要对其进行一个十分重要的操作，即去批次效应。

下面，我们一起来看一下如何去除多芯片合并时的批次效应。

1.GEO数据库的下载

关于数据下载的三种方法，就不做具体介绍了。在此，使用GEOquery包进行数据下载，并通过GSE33335和GSE56807两个数据集的合并操作来进行展示。

rm(list=ls()) #清空环境变量options(stringsAsFactors = F)###加载R包library(GEOquery)

1.1 读取数据集一

首先，读取数据集，并提取相应的表达数据和临床信息。

##读取数据集一load("GSE33335_eSet.Rdata")exp_a <-exprs(gset[[1]])pdata_a <-pData(gset[[1]])

使用exprs()函数和pData()函数分别读取表达和临床数据；可以看到，在GSE33335数据集中，共包含了50例样品的相关资料，22011个探针信息。

table(pdata_a$description)

接着，我们来统计一下样本的分组信息，其中正常组织25例，肿瘤组织25例。

1.2 读取数据集二

load("GSE56807_eSet.Rdata")exp_b <-exprs(gset[[1]])pdata_b <-pData(gset[[1]])

相应的，读取数据集二GSE56807的表达数据和临床信息。在该数据集中，共包含了10例样本，17786个探针信息。

table(pdata_b$characteristics_ch1.2)

同样table()函数来统计分组情况；其中肿瘤组织5例，正常组织5例。

1.3 表达矩阵合并

由于两个数据集来自同一平台，因此其基因对于探针编号是相同的，先就不进行探针id的转换。

但是，在前面的信息中可以看到，尽管挑选的两个数据集使用的是同一个芯片平台，但是两者的探针数量却并不相同。因此，在表达矩阵合并之前，我们首先需要获取两者共同检测的探针信息。

probe<- intersect(rownames(exp_a), rownames(exp_b))exp_a<- exp_a[probe,]exp_b<- exp_b[probe,]

这样，此时得到的两个表达矩阵的行名是完全相同的了。但是，这一个操作也间接的使得GSE33335数据集损失了一部分基因的表达信息，这部分基因中也许包含了本来可能的显著差异表达基因。

因此，在这里建议大家在合并的过程中，尽可能选择同一平台的芯片数据集，尽可能的减少基因信息的损失。

exp <-as.data.frame(cbind(exp_a,exp_b))

最后，使用cbind()函数进行按行合并，得到了一个由60例样本，17786个探针组成的表达矩阵。

1.4 设置分组信息

group_list <- c(rep("Normal",25), rep("Tumor",25), rep("Tumor",5), rep("Normal",5))group_list = factor(group_list, levels = c("Normal","Tumor"))

根据前面得到的分组信息，设置新的分组group_list。

2.查看批次信息

在正式的去除批次效应之前，我们先通过三种不同的方法来整体水平看看简单合并后的表达分布情况。

2.1 箱线图

boxplot(exp,outline=FALSE, notch=T,col=group_list, las=2)

结果显示：可以看到，两批样品的表达组间明显存在差异。

2.2 聚类分析

dist_mat <- dist(t(exp))res.hc <- hclust(dist_mat, method = "complete")

使用dist()函数计算样品之间的距离，并进行聚类分析。

plot(res.hc, labels = colnames(exp))

接着，将聚类结果进行可视化展示。

plot(res.hc, labels = group_list)

如果显示样品名不够清晰的话，使用分组信息就可以清晰的看到，在肿瘤组织中混入了3例正常组织。

2.3 PCA分析

library(FactoMineR)library(factoextra)data <-PCA(t(exp), graph = FALSE)fviz_pca_ind(data,geom.ind = "point",col.ind = group_list,addEllipses = TRUE,legend.title = "Groups")

结果显示：两组样品相互混杂，这样会导致差异分析时差异基因数量大大减少。

3.消除批次效应

接下来，我们分别使用两个R包的函数，介绍一下如何消除两个数据集合并时产生的批次效应。

方法一：sva包的ComBat()函数

library(sva)batchType=c(rep(1,50),rep(2,10))model<- model.matrix(~as.factor(group_list))combat_exp<- ComBat(dat = exp,batch = batchType,mod = model,par.prior=TRUE)

根据两个数据集的相应信息，对参数batchType和model进行设置，结合表达数据exp，将其导入到ComBat()函数中，即可完成去批次操作。

由于去批次过程会对基因的表达数据进行相应的调整，因此我们在来对调整后的表达矩阵进行可视化展示。

boxplot(combat_exp, outline=FALSE, notch=T, col=group_list, las=2)

结果显示：

combat_dist <- dist(t(combat_exp))combat_hc <- hclust(combat_dist, method = "complete")plot(combat_hc, labels = colnames(exp))

结果显示：

plot(combat_hc, labels = group_list)

data <-PCA(t(combat_exp), graph = FALSE)fviz_pca_ind(data,geom.ind = "point",col.ind = group_list,addEllipses = TRUE,legend.title = "Groups")

结果显示：经过批次校正后，两个分组的样品之间的差异明显增加，但是仍有部分样品处于异常位置，对此可以进行适当的删除该异常样品。

方法二：limma包的removeBatchEffect()函数

library(limma)batchType=c(rep(1,50),rep(2,10))model<- model.matrix(~as.factor(group_list))rb_exp<- removeBatchEffect(exp,batch = batchType,design = model)

同样的，设置好batchType和model参数后，将其导入到removeBatchEffect()函数中进行分析。

boxplot(rb_exp, outline=FALSE, notch=T, col=group_list, las=2)

结果显示：