##Time:2017-10-8
##Author:Feng Shengyu
#----------------------------------------------------
#一、安裝必須的R包（推薦使用的R版本3.2.2）
#必須要安裝的包：
#? 1、clusterprofilter
#? source("http://bioconductor.org/biocLite.R")
#? biocLite("clusterprofilter")
#? 2、org.Mm.eg.db/org.Hs.eg.db(對(duì)應(yīng)需要研究的物種-小鼠/人)
#? biocLite("org.Mm.eg.db")/biocLite("org.Hs.eg.db")
#? 3、DOSE
#? biocLite(DOSE)
library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
#二 change the type of gene
#使用的上游數(shù)據(jù)是RNA-seq做完的差異表達(dá)的基因列表
#example:
# 15431
# 244091
# 15430
# 319158
# 13871
# 109663
# 735269
# 378431
# 21384
# 105247262
#讀取gene list
gene <- read.table("C:\\Users\\Feng\\Desktop\\up_regulate_symbolGene.txt")
geneSymbol <- gene[,1]
geneSymbol
#轉(zhuǎn)化基因類型，一般用cufflinks做的結(jié)果是symbol,此時(shí)需要轉(zhuǎn)化為entrzid
geneEntrezID <- bitr(geneSymbol, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Mm.eg.db")
#可以同時(shí)轉(zhuǎn)為多個(gè)類型的基因
#geneEntrezID <- bitr(geneSymbol, fromType="SYMBOL", toType=c("ENTREZID","UNIPROT"), OrgDb="org.Mm.eg.db",)
#三、enrichment analysis
#GO富集分析
ego_cc <- enrichGO(gene = geneEntrezID[,2], #使用entrezID作為輸入
OrgDb=org.Mm.eg.db,
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE
)
setwd("F:\\生信工具大全\\R")
write.table(as.data.frame(ego_cc@result),file="test_CC.txt",sep="\t")
#KEGG富集分析
kk <- enrichKEGG(gene = geneEntrezID[,2],
organism ="mouse",
pvalueCutoff = 0.05,
qvalueCutoff = 0.01,
minGSSize = 1,
use_internal_data =FALSE
)
write.table(as.data.frame(kk@result), file="test_kk.txt",sep="\t")
#作圖展示結(jié)果
barplot(ego_cc, showCategory=15, title="EnrichmentGO_CC") #條狀圖，按p從小到大排的
dotplot(ego_BP,title="EnrichmentGO_CC_dot") #點(diǎn)圖，按富集的數(shù)從大到小的
#--------------------核心代碼-----------------------
setwd("F:\\碩士生\\GO和KEGG富集分析")
library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
gene <- read.table("C:\\Users\\Feng\\Desktop\\up_regulate.gene")
geneSymbol <- gene[,1]
geneEntrezID <- bitr(geneSymbol, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Mm.eg.db")
ego_cc <- enrichGO(gene = geneEntrezID[,2], #使用entrezID作為輸入
OrgDb=org.Mm.eg.db,
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE
)
write.table(as.data.frame(ego_cc@result),file="haimati_M_up_enrich_GO.txt",sep="\t")
barplot(ego_cc, showCategory=15, title="GO_Enrichment") #條狀圖，按p從小到大排的
ego_BP <- enrichKEGG(gene = geneEntrezID[,2],
organism ="mouse",? #http://www.genome.jp/kegg/catalog/org_list.html(species names)
pvalueCutoff = 0.05,
qvalueCutoff = 0.01,
minGSSize = 1,
use_internal_data =FALSE
)
write.table(as.data.frame(ego_BP@result), file="haimati_M_up_enrich_KEGG.txt",sep="\t")
dotplot(ego_BP,title="EnrichmentGO_CC_dot") #點(diǎn)圖，按富集的數(shù)從大到小的
