去除芯片和樣本間批次效應(yīng)

老大在群里出的題，說感覺這個熱圖很詭異，然后中間我自己沒有用boxplot查看數(shù)據(jù)的表達(dá)量，對于數(shù)據(jù)不能有正確的認(rèn)識，導(dǎo)致一開始的deg的logFC都沒有達(dá)到正負(fù)1的，最重要的是：

1.是因為我知道什么情況下用log函數(shù)，然后我還在這里面錯誤的用了log函數(shù)；

2.不能用[1:4,1:4]查看數(shù)據(jù)集機構(gòu)；

3.normalizeBetweenArrays和removeBatchEffect函數(shù)的用處。

老大幫助修改了代碼，關(guān)于normalizeBetweenArrays、removeBatchEffect和boxplot的，結(jié)果才清晰明了。

原文圖片

image

下載數(shù)據(jù)+準(zhǔn)備數(shù)據(jù)

rm(list = ls())
options(stringsAsFactors = F)

library(GEOquery)
eSet1 <- getGEO("GSE83521",
                destdir = '.',
                getGPL = F)
#(1)提取表達(dá)矩陣exp
exp1 <- exprs(eSet1[[1]])
exp1[1:4,1:4]
dim(exp1)
# 一定要看boxplot
boxplot(exp1,las=2)
pd1 <- pData(eSet1[[1]])
#(3)提取芯片平臺編號
gpl1 <- eSet1[[1]]@annotation
save(pd1,exp1,gpl1,file = "step1-1output.Rdata")
load("step1-1output.Rdata")
image
eSet2 <- getGEO("GSE89143",
                destdir = '.',
                getGPL = F)
#(1)提取表達(dá)矩陣exp
exp2 <- exprs(eSet2[[1]])
exp2[1:4,1:4]
dim(exp2)
boxplot(exp2,las=2)
image

從上面的箱線圖結(jié)果可以看到，數(shù)值的表達(dá)量并不在同一條水平線上，并且有成敗上千，也有零，很明顯是沒有經(jīng)過log的。這是需要把數(shù)據(jù)log后再用boxplot來看數(shù)據(jù)的分布，用boxplot來看數(shù)據(jù)的分布非常重要。不能僅僅用[1:4,1:4]來查看，因為[1:4,1:4]并不能看到整體的數(shù)據(jù)情況。關(guān)于為什么要log，是因為做差異分析的limma包要求表達(dá)矩陣中的數(shù)據(jù)是經(jīng)過log的。可以參考老大的這篇：關(guān)于limma包差異分析結(jié)果的logFC解釋

exp2 = log2(exp2+1)
boxplot(exp2,las=2)
image

接下來這個函數(shù)厲害了，從上面的圖中可以看到有一個樣本的中位數(shù)和其他樣本明顯不在一條水平顯示，這個normalizeBetweenArrays函數(shù)，可以把他拉回正常水平，normalizeBetweenArrays只能是在同一個數(shù)據(jù)集里面使用。

library(limma)
exp2=normalizeBetweenArrays(exp2)
boxplot(exp2,las=2)
image

從上面的箱線圖可以看到，exp2的數(shù)據(jù)的分布基本在一條水平線上。

接下來將exp2的數(shù)據(jù)保存

#(2)提取臨床信息
pd2 <- pData(eSet2[[1]])
#(3)提取芯片平臺編號
gpl2 <- eSet2[[1]]@annotation

#這些代碼是什么鬼東西，我給你注釋了。

index <- sort.int(pd2$characteristics_ch1,index.return = T)
class(index)
pd2 <- pd2[index$ix,]
exp2 <- exp2[,match(rownames(pd2),colnames(exp2))]

save(pd2,exp2,gpl2,file = "step1-2output.Rdata")
load("step1-2output.Rdata")

探針注釋

## 是同一個平臺，非常棒
gpl2
gpl1
if(T){
  library(GEOquery)
  gpl<- getGEO('GPL19978', destdir=".")
  dim(gpl)
  colnames(Table(gpl)) #查一下列明
  head(Table(gpl)[,c(1,2)])
  ids=Table(gpl)[,c(1,2)]
  ids <- ids[-c(1:2),]
  ids <- ids[-c(1:13),]
  save(ids,file='ids.Rdata')
}

看一下獲得的探針和circ_RNA的對應(yīng)關(guān)系

image

差異分析-去除批次效應(yīng)

x1 <- exp1[rownames(exp1) %in% ids$ID,]
x2 <- exp2[rownames(exp2) %in% ids$ID,]
boxplot(x1,las=2)
boxplot(x2,las=2)
cg=intersect(rownames(x1),rownames(x2))
x_merge=cbind(x1[cg,],x2[cg,])
  • 又有大招了，得到的這個merge后的表達(dá)矩陣x_merge，一定一定要用boxplot看一下，因為我們是將兩個數(shù)據(jù)集通過共同的探針合并了，因為是來自兩個數(shù)據(jù)，所以第一次在實際案例中接觸到了這個高大上的名次-去除批次效應(yīng)。
  • 那么就boxplot來看看！
boxplot(x_merge)
image

上面這張圖，就是非常明顯的看到了，由于后面的三個樣本就是來自另一個數(shù)據(jù)集的。從前面的boxplot(exp2)也可以看到他的表達(dá)量在8以上。

pd1$title
pd2$characteristics_ch1
group_list <- c(rep('tumor',6),rep('normal',6),rep(c('tumor', 'normal'),each=3))
gse <- c(rep('GSE83527',12),rep('GSE89143',6))
table(group_list,gse)
dat <- x_merge
library(sva)
library(limma)
## 使用 limma 的 removeBatchEffect 函數(shù)
dat[1:4,1:4]
batch <- c(rep('GSE83521',12),rep('GSE89143',6))
design=model.matrix(~group_list)
ex_b_limma <- removeBatchEffect(dat,
                                batch = batch,
                                design = design)
dim(ex_b_limma)
  • 這個時候一定要看 boxplot ， 不然就白學(xué)的了?。。。。?/li> 
boxplot(ex_b_limma)
image

  • 從上面的圖可以看到了，去除了批次效應(yīng)后，數(shù)據(jù)的表達(dá)水平。

  • 接下來就做差異分析

{
  library(limma)
  fit=lmFit(ex_b_limma,design)
  fit=eBayes(fit)
  options(digits = 4)
  topTable(fit,coef=2,adjust='BH')
  deg=topTable(fit,coef=2,adjust='BH',number = Inf)
}

得到的deg如下圖

image

??兩張圖主要是為了看logFC的值，我第一次把兩個數(shù)據(jù)全部log并且沒有進行一個normalizeBetweenArrays來去除樣本間的批次差異，而沒有一個logFC值是在>1和<-1的。

火山圖

if(T){
  nrDEG=deg
  head(nrDEG)
  attach(nrDEG)
  plot(logFC,-log10(P.Value))
  library(ggpubr)
  df=nrDEG
  df$v= -log10(P.Value) #df新增加一列'v',值為-log10(P.Value)
  ggscatter(df, x = "logFC", y = "v",size=0.5)

  df$g=ifelse(df$P.Value>0.05,'stable', #if 判斷：如果這一基因的P.Value>0.01，則為stable基因
              ifelse( df$logFC >1,'up', #接上句else 否則：接下來開始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調(diào)）基因
                      ifelse( df$logFC < -1,'down','stable') )#接上句else 否則：接下來開始判斷那些logFC <1.5 的基因，再if 判斷：如果logFC <1.5，則為down（下調(diào)）基因，否則為stable基因
  )
  table(df$g)
  df$name=rownames(df)
  head(df)
  ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
  ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
            label = "name", repel = T,
            #label.select = rownames(df)[df$g != 'stable'] ,
            label.select = head(rownames(deg)), #挑選一些基因在圖中顯示出來
            palette = c("#00AFBB", "#E7B800", "#FC4E07") )
  ggsave('volcano.png')

}
image

熱圖

  • 在畫特圖前，還有一小問題，那么就是這個探針的問題，我們可以看到原文的熱圖上的基因是hsa-circ-0034398,而我們id轉(zhuǎn)換后的表達(dá)矩陣的基因名是，如下圖標(biāo)記所示
image

學(xué)習(xí)群里的小伙伴將一個Alias的對應(yīng)表格分享到群里，把它讀進R里進行進一步的轉(zhuǎn)換，文件是ID.txt。

if(T){
  up <- df[df$g=='up',]
  down <- df[df$g =='down',]
  x <- rbind(up,down)
  x$ID <- rownames(x)
  #上面是為了提取出差異基因的子集
  y <- merge(x,ids,by='ID') #merge函數(shù)可以根據(jù)兩列共有的內(nèi)容進行合并，合并后含有共有的行名
  ex_b_limma2 <- ex_b_limma[match(y$ID,rownames(ex_b_limma)),]
  rownames(ex_b_limma2) <- y$circRNA
  #上面的函數(shù)寫的不咋好，命名變量不好命名，需要得到結(jié)果后看一下
  z <- read.csv('ID.txt',sep = '\t')
  tmp1 <- z[z$circRNA %in% rownames(ex_b_limma2),]
  ex_b_limma3 <- ex_b_limma2[match(tmp1$circRNA,rownames(ex_b_limma2)),]

  rownames(ex_b_limma3) <- tmp1$Alias
  library(pheatmap)
  pheatmap(ex_b_limma3,show_colnames =T,show_rownames = T) 
  n=t(scale(t(ex_b_limma3)))
  n[n>2]=2
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(SampleType=group_list,#這步就是可以實現(xiàn)兩個分組了
                GeoDatabase=gse)
  rownames(ac)=colnames(n) #將ac的行名也就分組信息 給到n的列名，即熱圖中位于上方的分組信息，這步很重要
  pheatmap(n,show_colnames =T,
           show_rownames = T,
           cluster_cols = F,
           annotation_col=ac,
           fontsize = 8,
           filename = 'deg-heatmap.png') #列名注釋信息為ac即分組信息
}
image

重要的如下：

1.理解boxplot的重要性，來看數(shù)據(jù)集是否需要log，以便后面才能用limma包進行差異分析

2.normalizeBetweenArrays只能是在同一個數(shù)據(jù)集里面用來去除樣本的差異，不同數(shù)據(jù)集需要用limma 的 removeBatchEffect函數(shù)去除批次效應(yīng)。
GSE83521/GSE89143去除批次效應(yīng) - 簡書 (jianshu.com)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容