GO and KEGG富集分析

GO

基因本体涉及的基因和基因产物词汇分为三大类，涵盖生物学的三个方面:
细胞组分(cellular component)CC:细胞的每个部分和细胞外环境。
分子功能(molecular function)MF:可以描述为分子水平的活性(activity)，如催化(catalytic)或结合(binding)活性。
生物过程(biological process)BP:生物学过程系指由一个或多个分子功能有序组合而产生的系列事件。其定义有广义和狭义之分，在词义上可以区分为泛指和特指。一般规律是，一个过程是由多个不同的步骤组成。需要指出的是，生物学过程与途径或通路(pathway)不是同一回事。

KEGG

GO分析好比是将基因分门别类放入一个个功能类群的篮子，而pathway则是将基因一个个具体放到代谢网络中的指定位置。

kegg物种支持查询

http://www.genome.jp/kegg/catalog/org_list.html
ssc #sus
hsa #human
mmu #Mus musculus (mouse)

#示例数据
gene = "TP53,BRCA1,KRAS,PTEN,RUNDC3A,MAPK8IP2,PSD,JPH3,PELI3,STMN3,CDK5R2,AP3B2,SVOP,GRIN1,SNCB,RLTPR,STX1B,RP11-1263C18.1,KCNC1,LRRC73,KCNT1,PRRT1,SEPT5,CTD-2228A4.1,RP11-482M8.3,CPLX1,PNMA3,CKMT1A,PPP1R3F,CABP1,CKMT1B,SCAMP5,PPFIA3,PDZD4,FBXO41,CHRM1,AC091878.1,SPRN,UBE2QL1,CACNA1I,FAIM2,SYP,ARHGDIG,SYT4,MADD,CEND1,RP11-95O2.1,DLG4,CNNM1,EIF4E1B,MAPRE3,PRKAG2-AS1"
gene = unlist(strsplit(gene,','))

#安装clusterProfiler包
#if (!requireNamespace("BiocManager", quietly = TRUE){install.packages("BiocManager")}
#BiocManager::install("clusterProfiler")
#BiocManager::install("org.Hs.eg.db")

library(clusterProfiler)
library("org.Hs.eg.db")

#gene symbol id to entrez id
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
gene = gene$ENTREZID
#下面这个也可以进行id转换
#entrezid <- select(org.Hs.eg.db, keys=gene, columns = 'ENTREZID', keytype = 'SYMBOL')

#GO分析
enrich_go <- enrichGO(gene=gene, 
                      OrgDb = 'org.Hs.eg.db', 
                      keyType = 'ENTREZID', 
                      ont = 'ALL', #BP,CC,MF
                      pvalueCutoff = 0.05, #pvalue cutoff
                      qvalueCutoff = 0.05, #qvalue cutoff
                      readable = TRUE #whether mapping gene ID to gene Name
                      )

write.csv(summary(enrich_go),"GO_enrich.csv",row.names =F)
barplot(enrich_go, showCategory = 30,title="Enrichment GO")
dotplot(enrich_go, showCategory = 30,orderBy="x",title="Enrichment GO")

#KEGG分析
enrich_kegg <- enrichKEGG(gene = gene, 
                          organism = 'hsa', # 对应 KEGG 物种缩写
                          keyType = 'kegg',
                          pvalueCutoff = 0.05,
                          qvalueCutoff = 0.05
                          )

barplot(enrich_kegg, showCategory = 30,title="Enrichment KEGG")
dotplot(enrich_kegg, showCategory = 30,title="Enrichment KEGG",orderBy="x")

#enrichKEGG没有readable这个参数?,自己写个函数实现类似功能的转换
#entrezid to symbol
entrezid_2_symbol <- function(gene_list){
  gene_list=bitr(unlist(strsplit(gene_list, '/')), fromType="ENTREZID", toType="SYMBOL", OrgDb="org.Hs.eg.db")
  gene_list=paste(gene_list$SYMBOL,collapse = '/')
  gene_list
}

kegg_result=summary(enrich_kegg)
kegg_result$geneID=unlist(lapply(kegg_result$geneID, entrezid_2_symbol))
write.csv(kegg_result,"KEGG_enrich.csv",row.names =F)