rm(list = ls())  
load(file = 'step4output.Rdata')

富集分析考验网速，因此给大家保存了Rdata

上课运行示例数据无需修改，在做自己的数据时请注意把本行之后的load()去掉

library(clusterProfiler)
library(dplyr)
library(ggplot2)
source("kegg_plot_function.R")
#source表示运行整个kegg_plot_function.R脚本，里面是一个function
#以up_kegg和down_kegg为输入数据做图

1.GO database 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)GO分析，分三部分

#以下步骤耗时很长，实际运行时注意把if后面的括号里F改成T，
#自己操作时需要改成T，时间久，课堂上未运行
if(F){
  #细胞组分
  ego_CC <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "CC",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  #生物过程
  ego_BP <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "BP",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  #分子功能：
  ego_MF <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "MF",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  save(ego_CC,ego_BP,ego_MF,file = "ego_GSE42872.Rdata")
}
load(file = "ego_GSE42872.Rdata")

(3)可视化

#条带图
barplot(ego_CC,showCategory=20)
#气泡图
dotplot(ego_CC)
#下面的图需要映射颜色，设置和示例数据一样的geneList
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)

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

#Gene-Concept Network
cnetplot(ego_CC, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego_CC, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map
emapplot(ego_CC)

(4)展示通路关系

goplot(ego_CC)

(5)Heatmap-like functional classification

heatplot(ego_CC,foldChange = geneList)
#太多基因就会糊。可通过调整比例或者减少基因来控制。
pdf("heatplot.pdf",width = 14,height = 5)
heatplot(ego_CC,foldChange = geneList)
dev.off()

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）对上调/下调/所有差异基因进行富集分析

注意这里又有个F

if(F){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      universe     = gene_all,
                      pvalueCutoff = 0.9,
                      qvalueCutoff = 0.9)
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        universe     = gene_all,
                        pvalueCutoff = 0.9,#可以自己设置，富集的结果条件
                        qvalueCutoff =0.9)
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.9)
  save(kk.diff,kk.down,kk.up,file = "GSE42872kegg.Rdata")
}
load("GSE42872kegg.Rdata")

(3)从富集结果中提取出结果数据框

kegg_diff_dt <- kk.diff@result

(4)按照pvalue筛选通路

#在enrichkegg时没有设置pvaluecutoff，在此处筛选
down_kegg <- kk.down@result %>%
  filter(pvalue<0.05) %>% #筛选行
  mutate(group=-1) #新增列
#插播，uniqe返回的不是逻辑值，duplicated去重复才是返回逻辑值，加！表示否定
#多个探针对应一个基因，筛选保留一个的时候，可以取均值或最大值，无正确答应
#table统计差异表达的基因值会重复统计
up_kegg <- kk.up@result %>%
  filter(pvalue<0.05) %>%
  mutate(group=1)

(5)可视化

g_kegg <- kegg_plot(up_kegg,down_kegg)

#g_kegg +scale_y_continuous(labels = c(20,15,10,5,0,5))
#设置叠加用+注意每个研究不一样，自己研究课题时，需要自己清楚
ggsave(g_kegg,filename = 'kegg_up_down.png')

#gsea作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',
                  nPerm        = 1000,
                  minGSSize    = 120,
                  pvalueCutoff = 0.9,#可以修改
                  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)可视化
kegg_plot(up_kegg,down_kegg)
ggsave('kegg_up_down_gsea.png')