RNA-seq入门实战（七）：GSEA——基因集富集分析

rm(list = ls())  
options(stringsAsFactors = F)
# library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(ggstatsplot)

setwd("C:/Users/Lenovo/Desktop/test")
source('FUNCTIONS.R')
load(list.files(path = "./3.DEG",pattern = 'DEG_results.Rdata',full.names = T))
dir.create("5.GSEA_kegg_go")
setwd("5.GSEA_kegg_go")

## 物种设置
organism = 'mmu'    #  人类'hsa' 小鼠'mmu'   
OrgDb = 'org.Mm.eg.db'#人类"org.Hs.eg.db" 小鼠"org.Mm.eg.db"

#### 按照需要可选择不同的DEG方法数据集 ####
need_DEG <- DEG_DESeq2
need_DEG <- need_DEG[,c(2,5)] #选择log2FoldChange和pvalue（凑成数据框）

colnames(need_DEG) <- c('log2FoldChange','pvalue')
need_DEG$SYMBOL <- rownames(need_DEG)

##### 创建gsea分析的geneList（包含从大到小排列的log2FoldChange和ENTREZID信息）####
#转化id  
df <- bitr(rownames(need_DEG), 
           fromType = "SYMBOL",
           toType =  "ENTREZID",
           OrgDb = OrgDb) #人数据库org.Hs.eg.db 小鼠org.Mm.eg.db
need_DEG <- merge(need_DEG, df, by='SYMBOL')  #按照SYMBOL合并注释信息
geneList <- need_DEG$log2FoldChange
names(geneList) <- need_DEG$ENTREZID
geneList <- sort(geneList, decreasing = T)   #从大到小排序

##### gsea富集 ####
KEGG_kk_entrez <- gseKEGG(geneList     = geneList,
                   organism     = organism, #人hsa 鼠mmu
                   pvalueCutoff = 0.25)  #实际为padj阈值可调整 
KEGG_kk <- DOSE::setReadable(KEGG_kk_entrez, 
                             OrgDb=OrgDb,
                             keyType='ENTREZID')#转化id             
  
GO_kk_entrez <- gseGO(geneList     = geneList,
               ont          = "ALL",  # "BP"、"MF"和"CC"或"ALL"
               OrgDb        = OrgDb,#人类org.Hs.eg.db 鼠org.Mm.eg.db
               keyType      = "ENTREZID",
               pvalueCutoff = 0.25)   #实际为padj阈值可调整
GO_kk <- DOSE::setReadable(GO_kk_entrez, 
                           OrgDb=OrgDb,
                           keyType='ENTREZID')#转化id 

save(KEGG_kk_entrez, GO_kk_entrez, file = "GSEA_result.RData")

##选取富集结果
kk_gse <- KEGG_kk
kk_gse_entrez <- KEGG_kk_entrez

###条件筛选 
#一般认为|NES|>1，NOM pvalue<0.05，FDR（padj）<0.25的通路是显著富集的
kk_gse_cut <- kk_gse[kk_gse$pvalue<0.05 & kk_gse$p.adjust<0.25 & abs(kk_gse$NES)>1]
kk_gse_cut_down <- kk_gse_cut[kk_gse_cut$NES < 0,]
kk_gse_cut_up <- kk_gse_cut[kk_gse_cut$NES > 0,]

#选择展现NES前几个通路 
down_gsea <- kk_gse_cut_down[tail(order(kk_gse_cut_down$NES,decreasing = T),10),]
up_gsea <- kk_gse_cut_up[head(order(kk_gse_cut_up$NES,decreasing = T),10),]
diff_gsea <- kk_gse_cut[head(order(abs(kk_gse_cut$NES),decreasing = T),10),]


#### 经典的GSEA图 
up_gsea$Description
i=2
gseap1 <- gseaplot2(kk_gse,
                    up_gsea$ID[i],#富集的ID编号
                    title = up_gsea$Description[i],#标题
                    color = "red", #GSEA线条颜色
                    base_size = 20,#基础字体大小
                    rel_heights = c(1.5, 0.5, 1),#副图的相对高度
                    subplots = 1:3,   #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
                    ES_geom = "line", #enrichment score用线还是用点"dot"
                    pvalue_table = T) #显示pvalue等信息
ggsave(gseap1, filename = 'GSEA_up_1.pdf', width =10, height =8)
 
#### 合并 GSEA通路 
gseap2 <- gseaplot2(kk_gse,
                    up_gsea$ID,#富集的ID编号
                    title = "UP_GSEA_all",#标题
                    color = "red",#GSEA线条颜色
                    base_size = 20,#基础字体大小
                    rel_heights = c(1.5, 0.5, 1),#副图的相对高度
                    subplots = 1:3, #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
                    ES_geom = "line",#enrichment score用线还是用点"dot"
                    pvalue_table = T) #显示pvalue等信息
ggsave(gseap2, filename = "GSEA_up_all.pdf",width =12,height =12)

## gsearank plot 绘制出属于特定基因集的基因排序列表
##绘制up_gsea前3个富集通路
library(cowplot)
library(ggplot2)
pp <- lapply(1:3, function(i) {
  anno <- up_gsea[i, c("NES", "pvalue", "p.adjust")]
  lab <- paste0(names(anno), "=",  round(anno, 3), collapse="\n")
  
  gsearank(kk_gse,
           up_gsea$ID[1], 
           title = up_gsea$Description[i]) + 
    xlab(NULL) +
    # ylab(NULL) +
    annotate("text", 10000,
             up_gsea[i, "enrichmentScore"] * .75, 
             label = lab, 
             hjust=0, vjust=0)
})
rankp <- plot_grid(plotlist=pp, ncol=1)
ggsave(rankp, filename = "gsearank_up.pdf",width=8,height=10)

## ridgeplot
ridgep <- ridgeplot(kk_gse_entrez,
                    showCategory = 15,
                    fill = "p.adjust",
                    core_enrichment = TRUE,
                    label_format = 30, #设置轴标签文字的每行字符数长度，过长则会自动换行。
                    orderBy = "NES",
                    decreasing = F) 
ggsave(ridgep,filename = 'ridgeplot.pdf',width =10,height =10)