老大讓復(fù)現(xiàn)的是一張韋恩圖研叫，來(lái)自3個(gè)數(shù)據(jù)集，分別找每個(gè)數(shù)據(jù)集的差異基因，然后取UP和DOWN的交集媒峡。原文鏈接：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6054764/#.
原文中的圖片如下
image-20191121212758565
想先用循環(huán)下載含有表達(dá)矩陣和臨床信息的一個(gè)list，然后試著用循環(huán)
library(GEOquery)
ls<-c('GSE38959','GSE45827','GSE65194')
for (i in 1:length(ls)){
gset[[i]] <- getGEO(ls[i], #系列編號(hào)
destdir = '.', #當(dāng)前目錄
getGPL = F)
}
但是下載下來(lái)的就只是表達(dá)矩陣葵擎，并不是含有臨床信息的list谅阿，下載后的表達(dá)矩陣如下
image-20191121212900565
接著想再嘗試循環(huán)讀取上面的表達(dá)矩陣
tmp<-dir(pattern = ".gz")
tmp
dat<-list()
for (i in 1:length(tmp)) {
dat[[i]]<-read.table(tmp[i],comment.char = '!',sep = '\t',header = T,row.names = 1)
}
上面就能得到含有三個(gè)data.frame的list了，上面兩個(gè)循環(huán)就只會(huì)到這里酬滤，因?yàn)檫€要再去下載數(shù)據(jù)集的臨床信息签餐，所以就分三個(gè)數(shù)據(jù)集分別下載，得到上下調(diào)基因然后做韋恩圖盯串。
image-20191121212816589
- 獲得第一個(gè)數(shù)據(jù)集的差異基因氯檐，下面是代碼
下載數(shù)據(jù)
rm(list = ls())
options(stringsAsFactors = F)
f='GSE38959_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38959
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE38959', destdir=".",
AnnotGPL = F,
getGPL = F)
save(gset,file=f)
}
load('GSE38959_eSet.Rdata')
class(gset)
length(gset)
class(gset[[1]])
a=gset[[1]]
dat=exprs(a) #a現(xiàn)在是一個(gè)對(duì)象，取a這個(gè)對(duì)象通過(guò)看說(shuō)明書(shū)知道要用exprs這個(gè)函數(shù)
dim(dat)
dat[1:4,1:4]
pd=pData(a) #通過(guò)查看說(shuō)明書(shū)知道取對(duì)象a里的臨床信息用pData
tnbc=rownames(pd[grepl('triple',as.character(pd$title)),])
normal=rownames(pd[grepl('normal',as.character(pd$title)),])
dat=dat[,c(tnbc,normal)]
group_list=c(rep('tnbc',length(tnbc)),
rep('normal',length(normal)))
table(group_list)
dat[1:4,1:4]
if(F){
# library(GEOquery)
# #Download GPL file, put it in the current directory, and load it:
# gpl <- getGEO('GPL4133', destdir=".")
# colnames(Table(gpl))
# head(Table(gpl)[,c(1,10)]) ## you need to check this , which column do you need
# probe2gene=Table(gpl)[,c(1,10)]
b=read.table('GPL4133-12599.txt',
sep = '\t',quote = '',
fill=T,header = T)
colnames(b)
probe2gene=b[,c(1,10)]
probe2gene=probe2gene[probe2gene$GENE_SYMBOL != '',]
head(probe2gene)
save(probe2gene,file='probe2gene.Rdata')
}
load(file='probe2gene.Rdata')
head(probe2gene)
ids=probe2gene
colnames(ids)=c('probe_id','symbol')
ids=ids[ids$probe_id %in% rownames(dat),]
ids$probe_id=as.character(ids$probe_id)
dat[1:4,1:4]
dat=dat[ids$probe_id,]
ids$median=apply(dat,1,median) #ids新建median這一列体捏，列名為median男摧，同時(shí)對(duì)dat這個(gè)矩陣按行操作，取每一行的中位數(shù)译打，將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對(duì)ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序耗拓，將對(duì)應(yīng)的行賦值為一個(gè)新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項(xiàng)，'!'為否奏司，即取出不重復(fù)的項(xiàng)乔询，去除重復(fù)的gene ，保留每個(gè)基因最大表達(dá)量結(jié)果s
dat=dat[ids$probe_id,] #新的ids取出probe_id這一列韵洋，將dat按照取出的這一列中的每一行組成一個(gè)新的dat
rownames(dat)=ids$symbol#把ids的symbol這一列中的每一行給dat作為dat的行名
dat[1:4,1:4] #保留每個(gè)基因ID第一次出現(xiàn)的信息
dim(dat)
boxplot(dat)
dat=log(dat+1)
dat[1:4,1:4]
save(dat,group_list,file = 'step1-output.Rdata')
load(file = 'step1-output.Rdata')
檢查數(shù)據(jù)
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
## 下面是畫PCA的必須操作竿刁，需要看說(shuō)明書(shū)。
dat=t(dat)#畫PCA圖時(shí)要求是行名時(shí)樣本名搪缨，列名時(shí)探針名食拜，因此此時(shí)需要轉(zhuǎn)換
dat=as.data.frame(dat)
dat=cbind(dat,group_list)
library("FactoMineR")
library("factoextra")
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#現(xiàn)在dat最后一列是group_list，需要重新賦值給一個(gè)dat.pca,這個(gè)矩陣是不含有分組信息的
fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = dat$group_list,
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
legend.title = "Groups"
)
ggsave('all_samples_PCA_by_pCR.png')
rm(list = ls())
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取副编，'2'是按列雀旱椤）取每一行的方差，從小到大排序，取最大的1000個(gè)
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對(duì)那些提取出來(lái)的1000個(gè)基因所在的每一行取出呻待，組合起來(lái)為一個(gè)新的表達(dá)矩陣
n=t(scale(t(dat[cg,]))) # 'scale'可以對(duì)log-ratio數(shù)值進(jìn)行歸一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,filename = 'heatmap_top1000_sd.png')
差異分析
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
table(group_list)
boxplot(dat[1,]~group_list)
bp=function(g){
library(ggpubr)
df=data.frame(gene=g,stage=group_list)
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
# Add p-value
p + stat_compare_means()
}
bp(dat[1,]) ## 調(diào)用上面定義好的函數(shù)打月，避免同樣的繪圖代碼重復(fù)多次敲。
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
dim(dat)
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
## 但是上面的用法做不到隨心所欲的指定任意兩組進(jìn)行比較
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)
exprSet=dat
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts("tnbc-normal",
levels = design)
contrast.matrix ##這個(gè)矩陣聲明蚕捉，我們要把 Tumor 組跟 Normal 進(jìn)行差異分析比較
colnames(design)
deg = function(exprSet,design,contrast.matrix){
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)
head(deg)
save(deg,file = 'deg.Rdata')
load(file = 'deg.Rdata')
head(deg)
bp(dat[rownames(deg)[1],])
## for volcano
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.01,'stable', #if 判斷：如果這一基因的P.Value>0.01奏篙，則為stable基因
ifelse( df$logFC >1,'up', #接上句else 否則：接下來(lái)開(kāi)始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調(diào)）基因
ifelse( df$logFC < -1,'down','stable') )#接上句else 否則：接下來(lái)開(kāi)始判斷那些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)), #挑選一些基因在圖中顯示出來(lái)
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano.png')
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )
ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
ggsave('MA.png')
}
## for heatmap
if(T){
load(file = 'step1-output.Rdata')
# 每次都要檢測(cè)數(shù)據(jù)
dat[1:4,1:4]
table(group_list)
x=deg$logFC #deg取logFC這列并將其重新賦值給x
names(x)=rownames(deg) #deg取probe_id這列敛熬，并將其作為名字給x
cg=c(names(head(sort(x),100)),#對(duì)x進(jìn)行從小到大排列充易，取前100及后100，并取其對(duì)應(yīng)的探針名荸型，作為向量賦值給cg
names(tail(sort(x),100)))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對(duì)dat按照cg取行盹靴，所得到的矩陣來(lái)畫熱圖
n=t(scale(t(dat[cg,])))#通過(guò)“scale”對(duì)log-ratio數(shù)值進(jìn)行歸一化，現(xiàn)在的dat是行名為探針瑞妇，列名為樣本名稿静，由于scale這個(gè)函數(shù)應(yīng)用在不同組數(shù)據(jù)間存在差異時(shí)，需要行名為樣本辕狰，因此需要用t(dat[cg,])來(lái)轉(zhuǎn)換改备，最后再轉(zhuǎn)換回來(lái)
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #將ac的行名也就分組信息（是‘no TNBC’還是‘TNBC’）給到n的列名，即熱圖中位于上方的分組信息
pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = T,
annotation_col=ac,filename = 'heatmap_top200_DEG.png') #列名注釋信息為ac即分組信息
}
write.csv(deg,file = 'deg.csv')
上面都是獲得了每個(gè)表達(dá)矩陣的差異分析結(jié)果蔓倍，上面的代碼在另外兩個(gè)數(shù)據(jù)集再走一波悬钳，當(dāng)然表達(dá)矩陣和分組信息需要自己制作好。每一個(gè)數(shù)據(jù)集得到的差異分析的“deg.Rdata”都保存在自己的文件夾偶翅，即使都叫“deg.Rdata”也沒(méi)關(guān)系默勾。下面就是做韋恩圖，獲得3個(gè)數(shù)據(jù)集的UP及DOWN基因的交集聚谁，篩選閾值可以按照文章中或自己調(diào)整
- 韋恩圖的代碼如下
rm(list = ls()) ## 魔幻操作母剥，一鍵清空~
# P<0.05 and |logFC|≥2.0,
# GSE38959, 515 upregulated genes and 337 downregulated genes.
# GSE45827, 2,117 genes were upregulated, and 878 genes were downregulated.
# GSE65194, 2,130 upregulated genes and 901 downregulated genes were identified.
load(file = '../GSE38959/deg.Rdata')
head(deg)
## 不同的閾值，篩選到的差異基因數(shù)量就不一樣形导，后面的超幾何分布檢驗(yàn)結(jié)果就大相徑庭环疼。
logFC_t=2
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC > logFC_t,'UP',
ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
up_GSE38959=rownames(deg[deg$g=='UP',])
down_GSE38959=rownames(deg[deg$g=='DOWN',])
load(file = '../GSE45827/deg.Rdata')
head(deg)
## 不同的閾值朵耕，篩選到的差異基因數(shù)量就不一樣伪阶，后面的超幾何分布檢驗(yàn)結(jié)果就大相徑庭形娇。
logFC_t=2
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC > logFC_t,'UP',
ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
up_GSE45827=rownames(deg[deg$g=='UP',])
down_GSE45827=rownames(deg[deg$g=='DOWN',])
load(file = '../GSE65194/deg.Rdata')
head(deg)
## 不同的閾值癣缅，篩選到的差異基因數(shù)量就不一樣友存，后面的超幾何分布檢驗(yàn)結(jié)果就大相徑庭直晨。
logFC_t=2
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC > logFC_t,'UP',
ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
up_GSE65194=rownames(deg[deg$g=='UP',])
down_GSE65194=rownames(deg[deg$g=='DOWN',])
save(up_GSE65194,up_GSE45827,up_GSE38959,
down_GSE65194,down_GSE45827,down_GSE38959,
file = 'all_up_down.rdata')
load(file = 'all_up_down.rdata')
library(VennDiagram)
venn.diagram(list( up_GSE65194=up_GSE65194 ,
up_GSE45827=up_GSE45827,
up_GSE38959=up_GSE38959 ),
fill=c("red","green","blue"),
filename="up_VennDiagram.tiff")
venn.diagram(list( down_GSE65194=down_GSE65194 ,
down_GSE45827=down_GSE45827,
down_GSE38959=down_GSE38959 ),
fill=c("red","green","blue"),
filename="down_VennDiagram.tiff")
image-20191121212838253
