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

连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中几个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！

他前面的分享是：

下面是他对我们b站转录组视频课程的详细笔记

承接上节：RNA-seq入门实战（四）：差异分析前的准备——数据检查，以及 RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较 

本节概览：

1.GSEA简单介绍

2.创建GSEA分析所需的geneList，包含log2FoldChange和ENTREZID信息

3.利用clusterProfiler进行GSEA富集GO与KEGG通路

4.GSEA富集结果可视化：GSEA结果图、 gsearank plot 、ridgeplot

1. GSEA简单介绍

以下对GSEA涉及的一些重要概念进行了简单介绍，详细介绍见：
一文掌握GSEA，超详细教程 - 云+社区 - 腾讯云 (tencent.com)史上最全GSEA可视化教程，今天让你彻底搞懂GSEA！ - 知乎 (zhihu.com)

1.1 GSEA定义与基本原理：

• 定义：
  基因集富集分析(Gene Set Enrichment Analysis, GSEA)是一种计算方法，用来确定一组先验定义的基因集是否在两种生物状态之间显示出统计学上显著的、一致的差异。
  官网地址：GSEA (gsea-msigdb.org)

• 基本原理：
  使用预定义的基因集（通常来自功能注释或先前实验的结果），将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合富集分析检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果

• 与GO\KEGG差异基因富集分析区别：
  差异基因富集分析是先筛选差异基因，再判断差异基因在哪些注释的通路存在富集；这涉及到阈值的设定，存在一定主观性并且只能用于表达变化较大的基因，即我们定义的显著差异基因。而GSEA则不局限于差异基因，从基因集的富集角度出发，理论上更容易囊括细微但协调性的变化对生物通路的影响。

1.2 MSigDB（Molecular Signatures Database）:

分子特征数据库。一般进行GSEA或GSVA使用的就是该数据库中的基因集，我们也可以自定义基因集。MSigDB所包含的基因集如下所示:

• KEGG信息包含在C2中，GO信息包含在C5中。

1.3 GSEA中关键概念

• ES（Enrichment Score）：富集得分ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是，从基因集L的第一个基因开始，计算一个累计统计值。当遇到一个落在s里面的基因，则增加统计值。遇到一个不在s里面的基因，则降低统计值。 每一步统计值增加或减少的幅度与基因的表达变化程度（fold-change值）是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集，负值ES表示基因集在列表的底部富集。
  p-value用来评估富集得分(ES)的显著性，通过排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。

• NES (Normalized Enrichment Score)：标准化富集得分每个基因子集s计算得到的ES根据基因集的大小进行标准化得到标准化富集得分Normalized Enrichment Score (NES)。随后会针对NES计算假阳性率FDR。

• Leading-edge subset：领头基因亚集对富集贡献最大的基因成员

• 一般认为|NES|>1，p-value<0.05，FDR<0.25的通路是显著富集的。
  |NES|值越大，FDR值就越小，说明分析的结果可信度越高。

---

2. 创建GSEA分析所需的geneList

在了解了GSEA基本概念后就可以正式开始实操了，首先需要将基因按照在两类样本中的差异表达程度排序。
下面我们构建包含了geneList，里面含有从大到小排序的log2FoldChange和对应的ENTREZID信息：

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)   #从大到小排序

---

3. 利用clusterProfiler包进行GSEA富集

clusterProfiler包内的gseGO()和gseKEGG()函数可以很方便地对GO与KEGG通路进行GSEA， 再使用DOSE::setReadable转化id 。

##### 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")

---

4. GSEA富集结果可视化

GSEA的可视化主要是GSEA结果图、 gsearank plot和ridgeplot山脊图。
同样也可以进行其他可视化如barplot、dotplot、cnetplot等等，详见RNA-seq入门的简单实战（六）：GO、KEGG富集分析与超全可视化攻略或者参阅说明书Chapter 15 Visualization of functional enrichment result | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)，这里就不再进行展示啦

4.1 gseaplot GSEA结果图

下面选取KEGG通路的富集结果进行gseaplot绘图示范

• 首先对富集结果进行条件筛选，一般认为|NES|>1，NOM pvalue<0.05，FDR（padj）<0.25的通路是显著富集的；还可以从结果中细分出上下调通路单独绘图，以下代码仅展示KEGG通路富集结果的上调通路。

• gseaplot2(）函数既可以对单独的通路绘图，也可以合并几个通路一起绘图；各类详细参数设置见以下代码处

##选取富集结果
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)

下面解释一下GSEA图的含义：

• 第1部分是ES折线图，离垂直距离x=0轴最远的峰值便是基因集的ES值，峰出现在排序基因集的前端（ES值大于0）则说明通路上调，出现在后端（ES值小于0）则说明通路下调。

• 第二部分为基因集成员位置图，用竖线标记了基因集中各成员出现在基因排序列表中的位置。若竖线集中分布在基因排序列表的前端或后端，说明该基因集通路上调或下调；若竖线较均匀分布在基因排序列表中，则说明该基因集通路在比较的两个数据中无明显变化。
  红色部分对应的基因在实验组中高表达，蓝色部分对应的基因在对照组中高表达，
  leading edge subset 是（0,0）到曲线峰值ES出现对应的这部分基因成员。

• 第三部分是排序后所有基因rank值（由log2FoldChang值计算得出）的分布，以灰色面积图显展示。

4.2 gsearank plot 绘制特定基因集的基因排序列表

gsearank()展示特定基因集的排序，横坐标为基因排序，纵坐标为ES值，利用cowplot和ggplot2包可以批量出图。

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

4.3 ridgeplot山脊图

展示富集通路的核心富集基因的表达分布，x轴为富集通路的核心富集基因表达变化的log2倍，值为正值表示表达上调，值为负值表示表达下调。

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

（之前运行报错解决方法见ridgeplot报错：Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : replacement has ... ）

4.4 其他富集结果可视化图

dotplot cnetplot emapplot treeplot heatplot upsetplot
详见RNA-seq入门的简单实战（六）：GO、KEGG富集分析与超全可视化攻略

---

GSEA分析和可视化到这就结束啦，下一节介绍GSVA的使用

参考资料

📖 Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)

一文掌握GSEA，超详细教程 - 云+社区 - 腾讯云 (tencent.com)

史上最全GSEA可视化教程，今天让你彻底搞懂GSEA！ - 知乎 (zhihu.com)

GitHub - jmzeng1314/GEO

【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili

【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年5~6月场，你的生物信息学第一课