R語言怎么處理矯正GEO數(shù)據(jù)批次效應(yīng)

本文小編為大家詳細(xì)介紹“R語言怎么處理矯正GEO數(shù)據(jù)批次效應(yīng)”，內(nèi)容詳細(xì)，步驟清晰，細(xì)節(jié)處理妥當(dāng)，希望這篇“R語言怎么處理矯正GEO數(shù)據(jù)批次效應(yīng)”文章能幫助大家解決疑惑，下面跟著小編的思路慢慢深入，一起來學(xué)習(xí)新知識(shí)吧。

什么是batch effect？

不同平臺(tái)的數(shù)據(jù)，同一平臺(tái)的不同時(shí)期的數(shù)據(jù)，同一個(gè)樣品不同試劑的數(shù)據(jù)，以及同一個(gè)樣品不同時(shí)間的數(shù)據(jù)等等都會(huì)產(chǎn)生一種batch effect 。這種影響如果廣泛存在應(yīng)該被足夠重視，否則會(huì)導(dǎo)致整個(gè)實(shí)驗(yàn)和最終的結(jié)論失敗。

我簡單說下什么叫做batch effect。比對(duì)實(shí)驗(yàn)組和對(duì)照組，不同的處理是患病和不患?。y(cè)序時(shí)，先測(cè)得疾病，然后測(cè)得正常），然后你通過分析，得到很多差異表達(dá)的基因?，F(xiàn)在問題來了，這個(gè)差異表達(dá)的結(jié)果是和你要研究的因素有關(guān)，還是時(shí)間有關(guān)，這個(gè)問題里時(shí)間就會(huì)成為干擾實(shí)驗(yàn)結(jié)果的因素，這個(gè)效應(yīng)就是batch effect。

如何檢測(cè)是否存在這種效應(yīng)呢

最簡單的就是記錄實(shí)驗(yàn)中時(shí)間這個(gè)變量，然后對(duì)差異表達(dá)的基因進(jìn)行聚類，看是否都和時(shí)間相關(guān)，如果相關(guān)就證明存在batch effect。

同樣，如果不同平臺(tái)的數(shù)據(jù)之間存在batch effect ，就不能簡單的整合。

大家可能都會(huì)問標(biāo)準(zhǔn)化，會(huì)不會(huì)處理掉batch effect ？

答案是能夠減弱，不能從根本上消除。如下圖，b是a進(jìn)行過標(biāo)準(zhǔn)化的結(jié)果，從樣本上看都一直，沒有什么問題，但是落實(shí)到基因?qū)用?，c圖中還是有明顯的batch effect，d圖中通過時(shí)間進(jìn)行聚類，很明顯可以看出差異表達(dá)主要是由于時(shí)間引起的。

矯正批次效應(yīng)

假如解決了這個(gè)批次問題，不僅可以讓實(shí)驗(yàn)更可靠，更厲害的是，我們可以做多個(gè)芯片的聯(lián)合分析了。矯正批次效應(yīng)有兩種方法：

1.使用sva的中combat包來校正批次效應(yīng)

下面是舉例子： 安裝必要的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) # 表達(dá)矩陣

看一下pheno里面有54行，4列構(gòu)成，里面記錄了批次信息 

有沒有批次效應(yīng)可以通過使用Hierarchical clustering的聚類方法去看一下聚類的情況：例如下面數(shù)據(jù)中，批次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)

校正批次效應(yīng)：model可以有也可以沒有，如果有，也就是告訴combat，有些分組本來就有差別，不要給我矯枉過正！

#再做一個(gè)分組列，用于批次效應(yīng)中排除項(xiàng)。
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 函數(shù)

還是利用上面的數(shù)據(jù)利用李

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數(shù)據(jù)批次效應(yīng)”文章已經(jīng)介紹完畢，想要掌握這篇文章的知識(shí)點(diǎn)還需要大家自己動(dòng)手實(shí)踐使用過才能領(lǐng)會(huì)，如果想了解更多相關(guān)內(nèi)容的文章，歡迎關(guān)注億速云行業(yè)資訊頻道。