富集分析-KEGG/GO

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

#(1)输入数据 组织输入数据
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)

#以下步骤耗时很长，设置了存在即跳过
f = paste0(gse_number,"_GO.Rdata") #如果存在文件ff，就会跳过，直接load(f)
if(!file.exists(f)){
  ego <- enrichGO(gene = gene_diff, 
                  OrgDb= org.Hs.eg.db,  #转换ID的数据库
                  ont = "ALL", #取子集
                  readable = TRUE) #可以直接讲结果转化为simple
  ego_BP <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "BP", #单独做某一个分支，取子集
                  readable = TRUE)
  #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,ego_BP,file = f)
}
load(f)

barplot(ego) #直接拿富集结果来画图
barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")  #结果看图a

dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

#gl 用于设置下图的颜色
gl = deg$logFC
names(gl)=deg$ENTREZID
#Gene-Concept Network,要放大看
cnetplot(ego,
         #layout = "star", #用layout()就换了一种视图圆形的(图f)
         color.params = list(foldChange = gl),
         showCategory = 3)
#结果图e

goplot(ego_BP) # 用于展示通路-通路之间的联系(图g)

#上调、下调、差异、所有基因
#（1）输入数据 与GO不需要作区分是完全一样的
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

#（2）对上调/下调/所有差异基因进行富集分析 （包含存在即跳过if...)
f2 = paste0(gse_number,"_KEGG.Rdata")
if(!file.exists(f2)){
  kk.up <- enrichKEGG(gene         = gene_up, #富集-enrichKEGG()
                      organism     = 'hsa') #指定物种的名称，人类-hsa，Go中是ORGDB
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')  #一共做了三次富集
  save(kk.diff,kk.down,kk.up,file = f2)
}
load(f2)
#KEGG 经常存在一个问题就是没有富集到，所以要看看富集到了没
#(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
> table(kk.diff@result$p.adjust<0.05)
FALSE  TRUE 
  320     2 
> table(kk.up@result$p.adjust<0.05)
FALSE  TRUE 
  292     3 
> table(kk.down@result$p.adjust<0.05) 
FALSE
  290 
  #没有富集到 不满足p.adjust值，是正常的(这是Y叔的包，只能用p.adjust富集，但有的满足p值，可以通过p值挽救一下。或者换一种方式富集，或者调整差异基因的阈值后面会有)
# 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值，在文中说明即可

source("kegg_plot_function.R") #source 是在不打开脚本的情况下运行，也是画图函数，是自定义函数（自己写的）
g_kegg <- kegg_plot(kk.up,kk.down) #可以根据不同阈值画图看图h，默认p值为0.05
g_kegg
#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6)) #强行修改横坐标上的标签，删掉#可把图h中的-2改成2

illma::removeBatchEffect()
sva::ComBat()