R语言怎么处理矫正GEO数据批次效应

本文小编为大家详细介绍“R语言怎么处理矫正GEO数据批次效应”，内容详细，步骤清晰，细节处理妥当，希望这篇“R语言怎么处理矫正GEO数据批次效应”文章能帮助大家解决疑惑，下面跟着小编的思路慢慢深入，一起来学习新知识吧。

什么是batch effect？

不同平台的数据，同一平台的不同时期的数据，同一个样品不同试剂的数据，以及同一个样品不同时间的数据等等都会产生一种batch effect 。这种影响如果广泛存在应该被足够重视，否则会导致整个实验和最终的结论失败。

我简单说下什么叫做batch effect。比对实验组和对照组，不同的处理是患病和不患病（测序时，先测得疾病，然后测得正常），然后你通过分析，得到很多差异表达的基因。现在问题来了，这个差异表达的结果是和你要研究的因素有关，还是时间有关，这个问题里时间就会成为干扰实验结果的因素，这个效应就是batch effect。

如何检测是否存在这种效应呢

最简单的就是记录实验中时间这个变量，然后对差异表达的基因进行聚类，看是否都和时间相关，如果相关就证明存在batch effect。

同样，如果不同平台的数据之间存在batch effect ，就不能简单的整合。

大家可能都会问标准化，会不会处理掉batch effect ？

答案是能够减弱，不能从根本上消除。如下图，b是a进行过标准化的结果，从样本上看都一直，没有什么问题，但是落实到基因层面，c图中还是有明显的batch effect，d图中通过时间进行聚类，很明显可以看出差异表达主要是由于时间引起的。

矫正批次效应

假如解决了这个批次问题，不仅可以让实验更可靠，更厉害的是，我们可以做多个芯片的联合分析了。矫正批次效应有两种方法：

1.使用sva的中combat包来校正批次效应

下面是举例子： 安装必要的R包并加载，comat就在sva包中。

# 安装包,提前添加镜像，加快安装速度
if (!requireNamespace("BiocManager", quietly=TRUE)){
  install.packages("BiocManager")
}
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
local({r <- getOption("repos")  
r["CRAN"] <- "http://mirrors.tuna.tsinghua.edu.cn/CRAN/"   
options(repos=r)}) 
BiocManager::install("sva")
BiocManager::install("bladderbatch")
library(sva)
library(bladderbatch)
library('Biobase')
library('GEOquery')
data(bladderdata)
#bladder 的属性是EsetExpressionSet，所以可以用pData和exprs方法
pheno <- pData(bladderEset) # 注释信息
edata <- exprs(bladderEset) # 表达矩阵

看一下pheno里面有54行，4列构成，里面记录了批次信息 

有没有批次效应可以通过使用Hierarchical clustering的聚类方法去看一下聚类的情况：例如下面数据中，批次1中cancer样品与normal有混合，需要矫正一下：

dist_mat <- dist(t(edata))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow=c(2,1))
plot(clustering, labels = pheno$batch)
plot(clustering, labels = pheno$cancer)

校正批次效应：model可以有也可以没有，如果有，也就是告诉combat，有些分组本来就有差别，不要给我矫枉过正！

#再做一个分组列，用于批次效应中排除项。
pheno$hasCancer <- as.numeric(pheno$cancer == "Cancer")  
#分组模型
model <- model.matrix(~hasCancer, data = pheno)
combat_edata <- ComBat(dat = edata, batch = pheno$batch, mod = model)
dist_mat_combat <- dist(t(combat_edata))
clustering_combat <- hclust(dist_mat_combat, method = "complete")
par(mfrow=c(2,1))
plot(clustering_combat, labels = pheno$batch)
plot(clustering_combat, labels = pheno$cancer)

2.使用 limma 的 removeBatchEffect 函数

还是利用上面的数据利用李

dat = edata, batch = pheno$batch

library("limma")
design <- model.matrix(~pheno$batch+pheno$hasCancer)
limma.edata <- removeBatchEffect(edata,
                                batch = pheno$batch,
                                design = design)

dist_mat_limma <- dist(t(limma.edata))
clustering_limma <- hclust(dist_mat_limma, method = "complete")
par(mfrow=c(2,1))
plot(clustering_limma, labels = pheno$batch)
plot(clustering_limma, labels = pheno$cancer)

读到这里，这篇“R语言怎么处理矫正GEO数据批次效应”文章已经介绍完毕，想要掌握这篇文章的知识点还需要大家自己动手实践使用过才能领会，如果想了解更多相关内容的文章，欢迎关注亿速云行业资讯频道。