1. #单细胞差异分析后的gsea富集分析
    2. rm(list = ls())
    3. library(Seurat)
    4. library(SeuratData)
    5. library(ggplot2)
    6. library(patchwork)
    7. library(dplyr)
    9. load(file = 'deg-by-MAST-for-mono-2-cluster.Rdata')
    10. colnames(deg)
    11. head(deg)
    13. #得到差异基因列表后取出 ，p值和logFC
    14. nrDEG=deg[,c(2,1)]
    15. colnames(nrDEG)=c('log2FoldChange','pvalue')
    16. head(nrDEG) 
    18. library(org.Hs.eg.db) 
    19. #Y叔的R包，把SYMBOL转换为ENTREZID，后面可以直接做 KEGG 和 GO
    20. library(clusterProfiler)
    21. gene <- bitr(rownames(nrDEG),         #转换的列是nrDEG的列名
    22.              fromType = "SYMBOL",     #需要转换ID类型
    23.              toType =  "ENTREZID",    #转换成的ID类型
    24.              OrgDb = org.Hs.eg.db)    #对应的物种，小鼠的是org.Mm.eg.db
    26. #让基因名、ENTREZID、logFC对应起来
    27. gene$logFC <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
    28. head(gene)
    29. geneList=gene$logFC
    30. names(geneList)=gene$ENTREZID 
    32. #按照logFC的值来排序geneList
    33. geneList=sort(geneList,decreasing = T)
    34. head(geneList)
    36. library(clusterProfiler)
    37. kk_gse <- gseKEGG(geneList     = geneList,
    38.                   organism     = 'hsa',
    39.                   nPerm        = 1000,
    40.                   minGSSize    = 10,
    41.                   pvalueCutoff = 0.9,
    42.                   verbose      = FALSE)
    44. #提取结果
    45. tmp=kk_gse@result
    46. kk=DOSE::setReadable(kk_gse,
    47.                      OrgDb='org.Hs.eg.db',
    48.                      keyType='ENTREZID')
    49. tmp=kk@result
    51. pro='test_gsea'
    52. write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
    55. if(F){
    56.   #按照enrichment score从高到低排序，便于查看富集通路
    57.   kk_gse=kk
    58.   sortkk<-kk_gse[order(kk_gse$enrichmentScore, decreasing = T),]
    59.   head(sortkk)
    60.   dim(sortkk)
    61.   #write.table(sortkk,"gsea_output2.txt",sep = "\\t",quote = F,col.names = T,row.names = F)
    62.   #可以根据自己想要的通路画出需要的图
    63.   library(enrichplot)
    64.   gseaplot(kk_gse, "hsa04510")
    65.   gseaplot2(kk_gse, "hsa04510", color = "firebrick",
    66.             rel_heights=c(1, .2, .6))#改变更多参数，为了美观
    69.   ##同时展示多个pathways的结果
    70.   #画出排名前四的通路
    71.   gseaplot2(kk_gse, row.names(sortkk)[1:4])
    73.   #上图用的是ES排名前4个画图，也可以用你自己感兴趣的通路画图
    74.   # 自己在刚才保存的txt文件里挑选就行。
    75.   paths <- c("hsa04510", "hsa04512", "hsa04974")
    76.   paths <- row.names(sortkk)[5:8]
    77.   paths
    78.   gseaplot2(kk_gse, paths)
    80.   #这里的GSEA分析其实由三个图构成，GSEA分析的runningES折线图
    81.   # 还有下面基因的位置图，还有所谓的蝴蝶图。如果不想同时展示，还可以通过subplots改变。
    82.   gseaplot2(kk_gse, paths, subplots=1)#只要第一个图
    83.   gseaplot2(kk_gse, paths, subplots=1:2)#只要第一和第二个图
    84.   gseaplot2(kk_gse, paths, subplots=c(1,3))#只要第一和第三个图
    86.   #如果想把p值标上去，也是可以的。
    87.   gseaplot2(kk_gse, paths, color = colorspace::rainbow_hcl(4),
    88.             pvalue_table = TRUE)
    90.   #最后的总结代码
    91.   gseaplot2(kk_gse,#数据
    92.             row.names(sortkk)[39],#画那一列的信号通路
    93.             title = "Prion disease",#标题
    94.             base_size = 15,#字体大小
    95.             color = "green",#线条的颜色
    96.             pvalue_table = TRUE,#加不加p值
    97.             ES_geom="line")#是用线，还是用d点
    99. }
    101. ###当然，这里不知道具体需要什么通路，就全部都画出来
    102. # 这里找不到显著下调的通路，可以选择调整阈值，或者其它。
    103. down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.4,];down_kegg$group=-1
    104. up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.4,];up_kegg$group=1
    105. dat=rbind(up_kegg,down_kegg)
    106. colnames(dat)
    107. dat$pvalue = -log10(dat$pvalue)
    108. dat$pvalue=dat$pvalue*dat$group 
    109. dat=dat[order(dat$pvalue,decreasing = F),]
    110. library(ggplot2)
    111. #条形图
    112. g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
    113.   geom_bar(stat="identity") + 
    114.   scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
    115.   scale_x_discrete(name ="Pathway names") +
    116.   scale_y_continuous(name ="log10P-value") +
    117.   coord_flip() + theme_bw(base_size = 15)+
    118.   theme(plot.title = element_text(hjust = 0.5),  axis.text.y = element_text(size = 15))+
    119.   ggtitle("Pathway Enrichment") 
    120. g_kegg
    121. print(g_kegg)
    122. ggsave(g_kegg,filename = paste0(pro,'_kegg_gsea.pdf'))
    124. library(enrichplot)
    125. gesa_res=kk@result
    127. ###画出每张kegg通路
    128. lapply(1:nrow(down_kegg), function(i){ 
    129.   gseaplot2(kk,down_kegg$ID[i],
    130.             title=down_kegg$Description[i],pvalue_table = T)
    131.   ggsave(paste0(pro,'_down_kegg_',
    132.                 gsub('/','-',down_kegg$Description[i])
    133.                 ,'.pdf'))
    134. })
    135. lapply(1:nrow(up_kegg), function(i){ 
    136.   gseaplot2(kk,up_kegg$ID[i],
    137.             title=up_kegg$Description[i],pvalue_table = T)
    138.   ggsave(paste0(pro,'_up_kegg_',
    139.                 gsub('/','-',up_kegg$Description[i]),
    140.                 '.pdf'))
    141. })
    145. ego <- gseGO(geneList     = geneList,
    146.              OrgDb        = org.Hs.eg.db,
    147.              ont          = "ALL",
    148.              nPerm        = 1000,   ## 排列数
    149.              minGSSize    = 100,
    150.              maxGSSize    = 500,
    151.              pvalueCutoff = 0.9,
    152.              verbose      = FALSE)  ## 不输出结果
    153. go=ego@result
    155. write.csv(go,file = 'gse_go.csv')