2020-04-20 Go analysis (clusterp

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