0.清空環(huán)境，加載R包

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

1.GO 富集分析

(1)輸入數(shù)據(jù)

#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)富集

#以下步驟耗時(shí)很長(zhǎng)，設(shè)置了存在即跳過
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)#gene ID自動(dòng)轉(zhuǎn)換成gene symbol
  #ont參數(shù)：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))#當(dāng)x軸太長(zhǎng)時(shí)設(shè)置了折疊

#geneList 用于設(shè)置下面圖的顏色
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,這個(gè)函數(shù)最近更新過，版本不同代碼會(huì)不同
Biobase::package.version("enrichplot")

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

(4)展示通路關(guān)系

https://zhuanlan.zhihu.com/p/99789859
goplot(ego)可以實(shí)現(xiàn)

(5)Heatmap-like functional classification

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

2.KEGG pathway analysis

（1）輸入數(shù)據(jù) 上調(diào)、下調(diào)、差異、所有基因

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）對(duì)上調(diào)/下調(diào)/所有差異基因進(jìn)行富集分析

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)#結(jié)果為FAUSE即沒有富集的結(jié)果，不要懷疑自己

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是不打開腳本的情況下運(yùn)行代碼
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4))#改橫坐標(biāo)軸，應(yīng)全為正值
ggsave(g_kegg,filename = 'kegg_up_down.png')

3.GSEA作kegg和GO富集分析

http://www.itdecent.cn/p/c5b7b7dbf29b
GSEA是把全部基因進(jìn)行富集，在GO和KEGG無(wú)法富集到的時(shí)候可以選擇

(1)查看示例數(shù)據(jù)

data(geneList, package="DOSE")

(2)將我們的數(shù)據(jù)轉(zhuǎn)換成示例數(shù)據(jù)的格式

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.能看懂的資料越來(lái)越多

GSEA學(xué)習(xí)更多：http://www.itdecent.cn/p/baf85b51752e
富集分析學(xué)習(xí)更多：http://yulab-smu.top/clusterProfiler-book/index.html
弦圖：http://www.itdecent.cn/p/e4bb41865b7f