老大在群里出的題，說感覺這個(gè)熱圖很詭異，然后中間我自己沒有用boxplot查看數(shù)據(jù)的表達(dá)量，對于數(shù)據(jù)不能有正確的認(rèn)識，導(dǎo)致一開始的deg的logFC都沒有達(dá)到正負(fù)1的，最重要的是：
1.是因?yàn)槲抑朗裁辞闆r下用log函數(shù)，然后我還在這里面錯(cuò)誤的用了log函數(shù)；
2.不能用[1:4,1:4]查看數(shù)據(jù)集機(jī)構(gòu)；
3.normalizeBetweenArrays和removeBatchEffect函數(shù)的用處。
老大幫助修改了代碼，結(jié)果才清晰明了。
原文圖片
GSE83521/GSE89143中的tup
下載數(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')
exp1的boxplot圖eSet2 <- getGEO('GSE89143',
destdir = '.',
getGPL = F)
#(1)提取表達(dá)矩陣exp
exp2 <- exprs(eSet2[[1]])
exp2[1:4,1:4]
dim(exp2)
boxplot(exp2,las=2)
exp2的boxplot圖從上面的箱線圖結(jié)果可以看到，數(shù)值的表達(dá)量并不在同一條水平線上，并且有成敗上千，也有零，很明顯是沒有經(jīng)過log的。這是需要把數(shù)據(jù)log后再用boxplot來看數(shù)據(jù)的分布，用boxplot來看數(shù)據(jù)的分布非常重要。不能僅僅用[1:4,1:4]來查看，因?yàn)閇1:4,1:4]并不能看到整體的數(shù)據(jù)情況。關(guān)于為什么要log，是因?yàn)樽霾町惙治龅膌imma包要求表達(dá)矩陣中的數(shù)據(jù)是經(jīng)過log的?？梢詤⒖祭洗蟮倪@篇：關(guān)于limma包差異分析結(jié)果的logFC解釋
exp2 = log2(exp2 1)
boxplot(exp2,las=2)
exp2進(jìn)行l(wèi)og后仍需normalizeBetweenArrays接下來這個(gè)函數(shù)厲害了，從上面的圖中可以看到有一個(gè)樣本的中位數(shù)和其他樣本明顯不在一條水平顯示，這個(gè)normalizeBetweenArrays函數(shù)，可以把他拉回正常水平，normalizeBetweenArrays只能是在同一個(gè)數(shù)據(jù)集里面使用。
library(limma)
exp2=normalizeBetweenArrays(exp2)
boxplot(exp2,las=2)
exp經(jīng)過normalizeBetweenArrays后的boxplot圖從上面的箱線圖可以看到，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')
探針注釋
## 是同一個(gè)平臺，非常棒
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)系
ids
差異分析-去除批次效應(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,])
boxplot(x_merge)
需要用removeBatchEffect去除批次效應(yīng)上面這張圖，就是非常明顯的看到了，由于后面的三個(gè)樣本就是來自另一個(gè)數(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(ex_b_limma)
ex_b_limma{
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如下圖
deg??這張圖主要是為了看logFC的值，我第一次把兩個(gè)數(shù)據(jù)全部log并且沒有進(jìn)行一個(gè)normalizeBetweenArrays來去除樣本間的批次差異，而沒有一個(gè)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')
}
火山圖熱圖
在畫特圖前，還有一小問題，那么就是這個(gè)探針的問題，我們可以看到原文的熱圖上的基因是hsa-circ-0034398,而我們id轉(zhuǎn)換后的表達(dá)矩陣的基因名是，如下圖標(biāo)記所示
ID.txt
學(xué)習(xí)群里的小伙伴將一個(gè)Alias的對應(yīng)表格分享到群里，把它讀進(jìn)R里進(jìn)行進(jìn)一步的轉(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)容進(jìn)行合并，合并后含有共有的行名
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,#這步就是可以實(shí)現(xiàn)兩個(gè)分組了
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即分組信息
}
熱圖重要的如下：
1.理解boxplot的重要性，來看數(shù)據(jù)集是否需要log，以便后面才能用limma包進(jìn)行差異分析
2.normalizeBetweenArrays只能是在同一個(gè)數(shù)據(jù)集里面用來去除樣本的差異，不同數(shù)據(jù)集需要用limma 的 removeBatchEffect函數(shù)去除批次效應(yīng)。
