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數據批次效應”文章已經介紹完畢，想要掌握這篇文章的知識點還需要大家自己動手實踐使用過才能領會，如果想了解更多相關內容的文章，歡迎關注億速云行業資訊頻道。