2020-04-20 Go analysis (clusterp
2020-04-20 本文已影响0人
xiaoguolaile
row.names=FALSE 行名不输出
或者
col.names=NA
readable= TRUE #gene ID 转换为gene symbol
############################################################################################
#若第一次使用请先安装clusterprofile包,若已安装请忽略此条
# install.packages("BiocManager") ##首先安装BiocManager，若已安装忽略此条
#
# library(BiocManager) ##加载安装包
#
# BiocManager::install("clusterProfiler") ##用BiocManager安装clusterprofiler
#
# BiocManager::install("org.Mm.eg.db") ##安装注释文件，此为小鼠的，其他物种见后文
#
# BiocManager::install("Rgraphviz") ## 安装作图包
#
# BiocManager::install("topGO") ##安装topGO
#
# BiocManager::install("pathview") ##安装pathview
#除注释外，所有符号必须是半角，即英文输入，包括括号、引号等，要不然会报错
#流程：设置工作路径----->读取数据------->作图-------->调整参数------>输出图片/数据
############################################################################################
setwd("C:/Users/zxg/Desktop/for_zyc") #设置工作目录，可在顶部Session--set working directory---choose..
library(clusterProfiler) #加载R包
library(org.Rn.eg.db) #加载注释包，根据物种选择安装并加载
library(Rgraphviz) #加载作图包
library(topGO) #加载GO包
library(pathview) #加载pathway包
x<- read.table("brain.txt",head=FALSE) #读取基因列表，基因ID暂定ENSEMBLE_ID,若更换类型需修改后续
genelist <- bitr(x[,1],fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Rn.eg.db") # id coversion
entrezid <- genelist$ENTREZID #extract entrezid
go <- enrichGO(entrezid , #gene list（entrez_id）
OrgDb = org.Rn.eg.db, #注释信息
ont='ALL', #注释方式
pAdjustMethod = 'BH', #富集方式
pvalueCutoff = 0.05, #p值设定
qvalueCutoff = 0.2, #q值设定
keyType = 'ENTREZID', #基因ID类型
readable= TRUE) #gene ID 转换为gene symbol
write.table(go,file="go_AD_up.xls",sep="\t",row.names = FALSE,col.names = )
pdf("AD_up.pdf")
barplot(go,showCategory=20,drop=T,main="whole") #整体作图，柱图，显示前20，前方数字可修改
dotplot(go,showCategory=20,title="whole") #整体作图，点图，显示前10，前方数字可修改
go.BP <- enrichGO(entrezid, OrgDb = org.Rn.eg.db, ont='BP', #仅BP-Biological process 富集
pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
go.MF <- enrichGO(entrezid, OrgDb = org.Rn.eg.db, ont='MF', #仅MF-Molecular function富集
pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
go.CC <- enrichGO(entrezid, OrgDb = org.Rn.eg.db, ont='CC', #仅CC-cell component富集
pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
write.table(go.BP,file="BP.xls",sep="\t") #BP富集信息写入表格，用制表符分隔，文件名字file=“”可修改
write.table(go.MF,file="MF.xls",sep="\t") #MF富集信息写入表格，用制表符分隔，文件名字file=“”可修改
write.table(go.CC,file="CC.xls",sep="\t") #CC富集信息写入表格，用制表符分隔，文件名字file=“”可修改
barplot(go.BP,showCategory=50,drop=T,title="BP") #BP 柱图
dotplot(go.BP,showCategory=10,title="BP") #BP 点图
plotGOgraph(go.BP) #BP 无向环图
barplot(go.MF,showCategory=10,drop=T,title="MF") #MF 柱图
dotplot(go.MF,showCategory=10,title="MF") #MF 点图
plotGOgraph(go.MF) #MF 无向环图
barplot(go.CC,showCategory=10,drop=T,title="CC") #CC 柱图
dotplot(go.CC,showCategory=10,title="CC") #CC 点图
plotGOgraph(go.CC) #CC 无向环图
kegg <- enrichKEGG(entrezid , #基因 entrezID
organism = 'mmu', #物种名称,去KEGG网站查询，人hsa，大鼠rno，小鼠mmu
#https://www.kegg.jp/kegg/catalog/org_list.html
keyType = 'kegg', #富集类型one of "kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'
pvalueCutoff= 0.05, #设定p值pvalue cutoff on enrichment tests to report
qvalueCutoff = 0.05, #设定q值qvalue cutoff on enrichment tests to report as significant.
pAdjustMethod = 'BH', #富集方式one of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
minGSSize = 10, #minimal size of genes annotated by Ontology term for testing
maxGSSize = 500, #maximal size of genes annotated for testing
use_internal_data = FALSE) #logical, use KEGG.db or latest online KEGG data
barplot(kegg,showCategory=10,drop=T) #KEGG 柱图
dotplot(kegg,showCategory=10) #KEGG 点图
plotGOgraph(kegg) #KEGG 无向环图
dev.off()
gsecc <- gseGO(geneList=geneList, ont="CC", OrgDb=org.Mm.eg.db, verbose=F)
gseaplot(gsecc, geneSetID="GO:0000779")
pathview(gene.data=entrezid, #either vector (single sample) or a matrix-like data (multiple sample).
pathway.id = 'rno03030', #the KEGG pathway ID(s), usually 5 digit, include the 3 letter KEGG species code
species="rno", # kegg code, scientific name or the common name
limit=list(gene=max(abs(entrezid)), #a list of two numeric elements with "gene" and "cpd" as the names.
cpd=1))
#######################################################################
#物种及对应注释包
#packages organism
#org.Ag.eg.db Anopheles
#org.At.tair.db Arabidopsis
#org.Bt.eg.db Bovine
#org.Ce.eg.db Worm
#org.Cf.eg.db Canine
#org.Dm.eg.db Fly
#org.Dr.eg.db Zebrafish
#org.EcK12.eg.db E coli strain K12
#org.EcSakai.eg.db E coli strain Sakai
#org.Gg.eg.db Chicken
#org.Hs.eg.db Human
#org.Mm.eg.db Mouse
#org.Mmu.eg.db Rhesus
#org.Pf.plasmo.db Malaria
#org.Pt.eg.db Chimp
#org.Rn.eg.db Rat
#org.Sc.sgd.db Yeast
#org.Ss.eg.db Pig
#org.Xl.eg.db Xenopus
