REACTOME 是一個(gè)開源、放開、人工整理并經(jīng)過同行評審的通路數(shù)據(jù)庫。目標(biāo)是為各個(gè)通路提供可視化、可解釋的分析，以支持基礎(chǔ)研究、臨床研究、基因組分析、建模、系統(tǒng)生物學(xué)等。
和我們所熟知的GO/KEGG類似，Reactome這些年也越來越的的在文章中出現(xiàn)，截止至筆者寫本文的時(shí)間，數(shù)據(jù)收錄信息如下：
不同于KEGG這種收費(fèi)的數(shù)據(jù)庫，Reactome的使用非常的簡單，現(xiàn)在讓我們一步步的開始搭建分析工具吧！
首選需要安裝以下包：clusterProfiler、ReactomePA、stringr、ggplot2，這些包均可以使用 BiocManager::install()或install.package()命令安裝這里就不再贅述。下面直接貼出代碼供大家使用。
suppressMessages(library(clusterProfiler))
suppressMessages(library(ReactomePA))
suppressMessages(library(stringr))
suppressMessages(library(ggplot2))
args <- commandArgs(T)
ref <- args[1]
genelist <- args[2]
outdir <- args[3]
top <- args[4]
GetGeneType<-function(df,filename){
if(str_detect(df[1,1],"EN")){
if( str_detect(df[1,1],"T") )
return("ENSEMBLTRANS")
else
return("ENSEMBL")
}else if(is.na(as.integer(df[1,1])) == F){
return("ENTREZID")
}else if(str_detect(df[1,1],"AT")){
return("TAIR")
}
else{
return("SYMBOL")
}
}
sample_name <- strsplit( basename(genelist) ,split=".",fixed=TRUE)[[1]][1]
df <- read.table(file=genelist,sep='\t',header=F)
setwd(outdir)
input_type = GetGeneType(df,genelist)
if(ref=='hs'){
suppressMessages(library(org.Hs.eg.db))
Db <- org.Hs.eg.db
organism <- 'hsa'
ref <- 'human'
}else if(ref=='mm'){
suppressMessages(library(org.Mm.eg.db))
Db <- org.Mm.eg.db
organism <- 'mmu'
ref <- 'mouse'
}else if( ref =='rn'){
suppressMessages(library(org.Rn.eg.db))
Db <- org.Rn.eg.db
organism <- 'rno'
ref <- 'rat'
}
bitrdf <- bitr(df[,1],fromType=input_type,toType="ENTREZID",OrgDb = Db)
Reactome <- as.data.frame(enrichPathway(gene=bitrdf[,'ENTREZID'] , pvalueCutoff = 1, readable=TRUE, organism = ref))
write.table(Reactome,paste0(sample_name,".Reactome.xls"),sep="\t",quote=F,row.names=F)
sigdf <- subset(Reactome, pvalue < 0.05)
write.table(sigdf,paste0(sample_name,".SigReactome.xls"),sep="\t",quote=F,row.names=F)
library(tidyverse)
plotdf <- separate(data = sigdf, col = BgRatio, into = c("BgCount", "BgAll"), sep = "/")
plotdf$BgCount = as.integer(plotdf$BgCount)
plotdf$RichFactor <- plotdf$Count / plotdf$BgCount
plotdf <- plotdf[order(plotdf$pvalue),]
plotdf <- head(plotdf,top)
p <- ggplot(data = plotdf,mapping = aes(x = RichFactor,y = reorder(Description,RichFactor)))+
geom_point(aes(color= pvalue,size = Count)) +
scale_colour_gradient(low = "red", high = "blue") +
theme_bw()+
labs(title = paste('Top',nrow(plotdf),'of Pathway Enrichment'),
x = 'Rich factor',
y = 'Pathway')+ scale_size("GeneNumber")+theme(plot.title = element_text(hjust = 0.5))
ggsave(paste0(sample_name,".bubble.png"), width = 8.5)
ggsave(paste0(sample_name,".bubble.pdf"), width = 8.5)
安裝好對應(yīng)的包后，將上述代碼保存為Reactome.r使用命令：
Rscript Reactome.r hs gene.glist ./ 20 即可進(jìn)行分析，其中
hs：物種信息，如人hs，小鼠mm，大鼠rn
gene.glist：基因列表一行一個(gè)基因，并去重，可以是ENSEMBL ID，ENTREZID, Symbol等，這里使用了bitrdf()將所有輸入全部轉(zhuǎn)為ENTREZID
./：輸出路徑，請先建立好該目錄，./代表當(dāng)前目錄
20：繪制TopN的氣泡圖
下面是實(shí)際演示：
輸入文件如下：
gene.glist文件
運(yùn)行
Rscript Reactome.r hs gene.glist ./ 20
輸出
輸出文件
*.bubble.png/pdf是氣泡圖
氣泡圖
*.Reactome.xls和 *,SigReactome.xls是輸出結(jié)果，前者是所有的富集結(jié)果后者是P<0.05的結(jié)果。
Reactome富集結(jié)果
這樣我們就完成了整個(gè)Reactome富集分析了，心動(dòng)了嗎？那快加入你的文章當(dāng)中吧！
