0.清空环境，加载R包

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

1.GO 富集分析

(1)输入数据

#head(deg)
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']

(2)富集

#以下步骤耗时很长，设置了存在即跳过
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)#gene ID自动转换成gene symbol
  #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,file = paste0(gse_number,"_GO.Rdata"))
}
load(paste0(gse_number,"_GO.Rdata"))
#class(ego)
#z=ego@result;z

(3)可视化

条带图

barplot(ego)

气泡图

dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + facet_grid(ONTOLOGY ~ ., scale = "free") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))#当x轴太长时设置了折叠

#geneList 用于设置下面图的颜色
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)

(3)展示top通路的共同基因，要放大看。

Gene-Concept Network

cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map,这个函数最近更新过，版本不同代码会不同
Biobase::package.version("enrichplot")

if(F){
  emapplot(pairwise_termsim(ego)) #新版本
}else{
  emapplot(ego)#老版本
}

(4)展示通路关系

https://zhuanlan.zhihu.com/p/99789859
goplot(ego)可以实现

(5)Heatmap-like functional classification

heatplot(ego,foldChange = geneList,showCategory = 8)

2.KEGG pathway analysis

（1）输入数据 上调、下调、差异、所有基因

gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']

（2）对上调/下调/所有差异基因进行富集分析

if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa')
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')
  save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
}
load(paste0(gse_number,"_KEGG.Rdata"))

(3)看看富集到了吗？

https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)#结果为FAUSE即没有富集的结果，不要怀疑自己

table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)

(4)p.adjust不适用可以按照pvalue筛选通路

down_kegg <- kk.down@result %>%
  filter(pvalue<0.05) %>% #筛选行
  mutate(group=-1) #新增列

up_kegg <- kk.up@result %>%
  filter(pvalue<0.05) %>%
  mutate(group=1)

(5)可视化

source("kegg_plot_function.R")#source是不打开脚本的情况下运行代码
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4))#改横坐标轴，应全为正值
ggsave(g_kegg,filename = 'kegg_up_down.png')

3.GSEA作kegg和GO富集分析

https://www.jianshu.com/p/c5b7b7dbf29b
GSEA是把全部基因进行富集，在GO和KEGG无法富集到的时候可以选择

(1)查看示例数据

data(geneList, package="DOSE")

(2)将我们的数据转换成示例数据的格式

geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)

(3)富集分析

kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  verbose      = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1

(4)可视化

g2 = kegg_plot(up_kegg,down_kegg)
g2

4.能看懂的资料越来越多

GSEA学习更多：https://www.jianshu.com/p/baf85b51752e
富集分析学习更多：http://yulab-smu.top/clusterProfiler-book/index.html
弦图：https://www.jianshu.com/p/e4bb41865b7f