1. rm(list = ls())  
    2. load(file = 'step4output.Rdata')
    3. library(clusterProfiler)
    4. library(ggthemes)
    5. library(org.Hs.eg.db)
    6. library(dplyr)
    7. library(ggplot2)
    8. library(stringr)
    9. library(enrichplot)
    1. # 1.GO 富集分析----
    3. #(1)输入数据
    4. gene_up = deg$ENTREZID[deg$change == 'up'] 
    5. gene_down = deg$ENTREZID[deg$change == 'down'] 
    6. gene_diff = c(gene_up,gene_down)
    8. #(2)富集
    9. #以下步骤耗时很长，设置了存在即跳过
    10. if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
    11.   ego <- enrichGO(gene = gene_diff,
    12.                   OrgDb= org.Hs.eg.db,
    13.                   ont = "ALL",
    14.                   readable = TRUE)
    15.   ego_BP <- enrichGO(gene = gene_diff,
    16.                   OrgDb= org.Hs.eg.db,
    17.                   ont = "BP",
    18.                   readable = TRUE)
    19.   #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
    20.   save(ego,ego_BP,file = paste0(gse_number,"_GO.Rdata"))
    21. }
    22. load(paste0(gse_number,"_GO.Rdata"))

    可视化

    1. #条带图
    2. barplot(ego)
    3. #气泡图
    4. dotplot(ego)
    6. dotplot(ego, split = "ONTOLOGY", font.size = 10, 
    7.         showCategory = 5) + 
    8.   facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") + 
    9.   scale_y_discrete(labels = function(x) str_wrap(x, width = 45))
    11. #geneList 用于设置下面图的颜色
    12. geneList = deg$logFC
    13. names(geneList)=deg$ENTREZID
    14. #(3)展示top通路的共同基因，要放大看。
    15. #Gene-Concept Network
    16. cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
    17. cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
    18. #Enrichment Map,这个函数最近更新过，版本不同代码会不同
    20. Biobase::package.version("enrichplot")
    22. if(T){
    23.   emapplot(pairwise_termsim(ego)) #新版本
    24. }else{
    25.   emapplot(ego)#旧版本
    26. }
    27. #(4)展示通路关系 https://zhuanlan.zhihu.com/p/99789859
    28. #goplot(ego)
    29. goplot(ego_BP)
    31. #(5)Heatmap-like functional classification
    32. heatplot(ego,foldChange = geneList,showCategory = 8)

    KEGG pathway analysis

    1. # 2.KEGG pathway analysis----
    2. #上调、下调、差异、所有基因
    3. #（1）输入数据
    4. gene_up = deg[deg$change == 'up','ENTREZID'] 
    5. gene_down = deg[deg$change == 'down','ENTREZID'] 
    6. gene_diff = c(gene_up,gene_down)
    8. #（2）对上调/下调/所有差异基因进行富集分析
    10. if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
    11.   kk.up <- enrichKEGG(gene         = gene_up,
    12.                       organism     = 'hsa')
    13.   kk.down <- enrichKEGG(gene         =  gene_down,
    14.                         organism     = 'hsa')
    15.   kk.diff <- enrichKEGG(gene         = gene_diff,
    16.                         organism     = 'hsa')
    17.   save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
    18. }
    19. load(paste0(gse_number,"_KEGG.Rdata"))
    21. #(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
    22. table(kk.diff@result$p.adjust<0.05)
    23. table(kk.up@result$p.adjust<0.05)
    24. table(kk.down@result$p.adjust<0.05)
    26. #(4)双向图
    27. # 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值，在文中说明即可
    28. down_kegg <- kk.down@result %>%
    29.   filter(pvalue<0.05) %>% #筛选行
    30.   mutate(group=-1) #新增列
    32. up_kegg <- kk.up@result %>%
    33.   filter(pvalue<0.05) %>%
    34.   mutate(group=1)
    36. source("kegg_plot_function.R")
    37. g_kegg <- kegg_plot(up_kegg,down_kegg)
    38. g_kegg
    39. #g_kegg +scale_y_continuous(labels = c(2,0,2,4,6))
    40. ggsave(g_kegg,filename = 'kegg_up_down.png')
    42. # 3.GSEA作kegg和GO富集分析----
    43. #https://www.yuque.com/xiaojiewanglezenmofenshen/dbwkg1/ytawgg
    45. #(1)查看示例数据
    46. data(geneList, package="DOSE")
    47. #(2)将我们的数据转换成示例数据的格式
    48. geneList=deg$logFC
    49. names(geneList)=deg$ENTREZID
    50. geneList=sort(geneList,decreasing = T)
    51. #(3)GSEA富集
    52. kk_gse <- gseKEGG(geneList     = geneList,
    53.                   organism     = 'hsa',
    54.                   verbose      = FALSE)
    55. down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
    56. up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
    57. #(4)可视化
    58. g2 = kegg_plot(up_kegg,down_kegg)
    59. g2