前面獲得的來自于不同celltype的marker，或者不同cluster的marker，或者每個celltype中不同分組的差異基因，上述marker構(gòu)成的列表形式，可以直接輸入，獲得富集結(jié)果的list以及圖片。

1. 加載包

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

2. 富集分析,org指定物種"mmu"或者"hsa"，markers_list為單細胞數(shù)據(jù)中獲得的maker的list

EnrichmentForMarkerlist <- function(markers_list,org,item.size=8,tag,path.out){
  if (org=="mmu") {
  db <- "org.Mm.eg.db"
}else if(org=="hsa"){
  db <- "org.Hs.eg.db"
}
GOenrich_result_list <- list()
KEGGenrich_result_list <- list()
for (i in names(markers_list)){
    print(paste("cluster",i,"enrichment",sep=" "))
    if (nrow(markers_list[[i]]) == 0){
        print(paste("nrow cluster",i,"is",nrow(markers_list[[i]]),sep=" "))
        next
    }
  #GO富集分析
  print(paste("cluster",i,"GO entichment",sep=" "))
    GOenrich_result <- enrichGO(rownames(markers_list[[i]]),db,keyType="SYMBOL", ont="all", pvalueCutoff=0.05)
#畫三種分類的圖
print(paste("cluster",i,"GO all entichment",sep=" "))
if (nrow(as.data.frame(GOenrich_result)) > 0){
    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("cluster_",i,"_GO_bar_all_plot.png"),GO_bar_split,path=path.out,width=15,height=10))
#只畫BP的圖
 print(paste("cluster",i,"GO BP entichment",sep=" "))
    GOenrich_result_BP <- enrichGO(rownames(markers_list[[i]]),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("cluster_",i,"_GO_bar_BP_plot.png"),GO_bar_BP,path=path.out,limitsize = FALSE))
}
#KEGG富集分析
    ID_transform_table <- bitr(rownames(markers_list[[i]]),fromType="SYMBOL",toType="ENTREZID",OrgDb=db)
    ID_Entrez_gene_set <- ID_transform_table$ENTREZID
    print("KEGGenrich start")
    KEGGenrich_result <- enrichKEGG(ID_Entrez_gene_set, organism = org, keyType = 'kegg', pvalueCutoff = 0.05,pAdjustMethod = 'BH', minGSSize = 10,maxGSSize = 500,use_internal_data = FALSE)
#繪制氣泡圖
print("KEGGenrich start map")
    try(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)))
    print("map saved")
    try(ggsave(paste0("cluster_",i,"_KEGG_dot_plot.png"),KEGG_dot,path=path.out,width=10,height=6))

try(GOenrich_result_list[[i]] <- GOenrich_result_BP)
try(KEGGenrich_result_list[[i]] <- KEGGenrich_result)
print("result saved")
}
print("start output")
  write.xlsx(GOenrich_result_list,file = paste0(path.out,"/",tag,'_GOenrich_result_list.xlsx'),row.names=T,sep='\t',overwrite=T)
  write.xlsx(KEGGenrich_result_list,file = paste0(path.out,"/",tag,'_KEGGenrich_result_list.xlsx'),row.names=T,sep='\t',overwrite=T)
  saveRDS(GOenrich_result_list,file=paste0(path.out,"/",tag,'_GO_EnrichResult.rds'))
  saveRDS(KEGGenrich_result_list,file=paste0(path.out,"/",tag,'_KEGG_EnrichResult.rds'))
}

3.加載數(shù)據(jù)

EnrichmentForMarkerlist(markers_list=all_cluster_markers,org="hsa",item.size=8,tag="clusters",path.out=path_out)