加载需要的包

library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(ggplot2)
library(openxlsx)

自定义函数

整体分析

• 输入为基因的symbol向量
• org是"mmu"或者"hsa"

EnrichmentForSymbol <- function(markers,org,item.size=8,tag,path.out){
if (org=="mmu") {
    db <- "org.Mm.eg.db"
}else if(org=="hsa"){
    db <- "org.Hs.eg.db"
}
#GO富集分析
    GOenrich_result <- enrichGO(markers,db,keyType="SYMBOL", ont="all", pvalueCutoff=0.05)
    write.xlsx(GOenrich_result,file = paste0(path.out,"/",tag,'_GOenrich_result.xlsx'),rowNames=T,sep='\t',overwrite=T)
    saveRDS(GOenrich_result,file=paste0(path.out,"/",tag,'_GO_EnrichResult.rds'))
#画三种分类的图
    GO_bar_split <- barplot(GOenrich_result, split = "ONTOLOGY",showCategory=20) + #GOenrich_result不能是矩阵
    facet_grid(ONTOLOGY~., scale = "free") + 
    scale_y_discrete(labels=function(x) str_wrap(x, width=100)) +  #取消了iterm强制换行，需要加载stringr包
    theme(axis.text.y=element_text(size=7)) 
    try(ggsave(paste0(tag,"_GO_bar_all_plot.png"),GO_bar_split,path=path.out,width=15,height=10))
#只画BP的图
    GOenrich_result_BP <- enrichGO(markers,db,keyType="SYMBOL", ont="BP", pvalueCutoff=0.05)
    GO_bar_BP <- barplot(GOenrich_result_BP,showCategory=30,drop=T) + 
    theme(axis.text.y=element_text(size=item.size))  + 
    scale_y_discrete(labels=function(x) str_wrap(x, width=100)) #取消了iterm强制换行，需要加载stringr包
    try(ggsave(paste0(tag,"_GO_bar_BP_plot.png"),GO_bar_BP,path=path.out,limitsize = FALSE))

#KEGG富集分析
    ID_transform_table <- bitr(markers,fromType="SYMBOL",toType="ENTREZID",OrgDb=db)
    ID_Entrez_gene_set <- ID_transform_table$ENTREZID
    KEGGenrich_result <- enrichKEGG(ID_Entrez_gene_set, organism = org, keyType = 'kegg', pvalueCutoff = 0.1,pAdjustMethod = 'BH', minGSSize = 10,maxGSSize = 500,use_internal_data = FALSE)
    write.xlsx(KEGGenrich_result,file = paste0(path.out,"/",tag,'_KEGGenrich_result.xlsx'),rowNames=T,sep='\t',overwrite=T)
    saveRDS(KEGGenrich_result,file=paste0(path.out,"/",tag,'_KEGG_EnrichResult.rds'))
#绘制气泡图
    KEGG_dot <- dotplot(KEGGenrich_result, showCategory=20) + theme(axis.text.y=element_text(size=item.size)) +
    scale_y_discrete(labels=function(x) str_wrap(x, width=100))
    try(ggsave(paste0(tag,"_KEGG_dot_plot.png"),KEGG_dot,path=path.out,width=10,height=6))
enrich_reult_list <- list(GO=GOenrich_result,KEGG=KEGGenrich_result)
return(enrich_reult_list)  
}

• 提取GO结果中的差异基因

ExportGOGene <- function(go.result){
fram <- data.frame()
for (i in seq(length(go.result$Description))){
    try(fram_1 <- data.frame(Description=go.result$Description[i],symbol=strsplit(go.result$geneID[i],split="/")[[1]],pvalue=go.result$pvalue[i]))
    try(fram <- rbind(fram,fram_1))
}
return(fram)
}

• 提取KEGG中的差异基因

ExportPathwayGene <- function(pathway_result){
fram <- data.frame()
for (i in seq(length(pathway_result$Description))){
    try(fram_1 <- data.frame(Description=pathway_result$Description[i],ENTREZID=strsplit(pathway_result$geneID[i],split="/")[[1]],pvalue=pathway_result$pvalue[i]))
    try(fram <- rbind(fram,fram_1))
}
transformID <- bitr(fram$ENTREZID,toType="SYMBOL",fromType="ENTREZID",OrgDb="org.Hs.eg.db")
symbol <- transformID$SYMBOL[match(fram$ENTREZID,transformID$ENTREZID)]
fram <- cbind(fram,symbol)#通路中基因的datafram
return(fram)
}