GEO数据挖掘 富集分析

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot) #南方医某个大佬写的R包 

# 初阶版本:手动修改if后面的F，需要运行的时候改成T
if(F){
  a = 1 #假装是限速步骤
  save(a,file = "a.Rdata")  #保存运行结果，下次运行到这里时直接加载结果
}
load("a.Rdata")

# 高阶版本
f = "a.Rdata"
if(!file.exists(f)){
  #只有f文件(a.Rdata)在工作目录下不存在时才运行，否则跳过这段代码
  a = 1 #假装是限速步骤
  print("bye") #如果运行了这个代码，就输出一个“bye”，起到提示作用
  save(a,file = f)
}
load(f)

# 完美版本：只修改q的内容就可，不同的数据对应不同的q
q = "124"
f = paste0(q,"a.Rdata")
if(!file.exists(f)){
  #只有f文件(a.Rdata)在工作目录下不存在时才运行，否则跳过这段代码
  a = 1 #假装是限速步骤
  print("bye")
  #保存运行结果，下次运行到这里时直接加载结果
  save(a,file = f)
}
load(f)

#1.输入数据
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)
#2.富集分析
f = paste0(gse_number,"_GO.Rdata")
if(!file.exists(f)){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db, #这里是人类物种的数据，如果是别的物种 需要修改
                  ont = "ALL",
                  readable = TRUE) #输出结果里含有gene symbol
  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) 
#3.画图
#条带图
barplot(ego)
barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
#气泡图
dotplot(ego)
dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
#4.展示top通路的共同基因
dotplot(ego)
dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

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

#2.对上调/下调/所有差异基因进行富集分析
f2 = paste0(gse_number,"_KEGG.Rdata")
if(!file.exists(f2)){
  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 = f2)
}
load(f2)

#3.看看富集到了吗？因为KEGG数据框包含的基因不多，所以是有可能没有很多差异基因的
table(kk.diff@result$p.adjust<0.05)
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)

#4.双向图
# 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值，在文中说明即可

source("kegg_plot_function.R") #在不打开脚本时，运行该脚本下的全部代码。
g_kegg <- kegg_plot(kk.up,kk.down) #这个函数可以调节FC、p值的界限
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4,6)) #修改坐标数字使其不带有负号，c里面的内容要根据富集的结果修改

limma::removeBatchEffect() 
sva::ConBat()