使用R语言进行GSEA富集分析

        GSEA（Gene Set Enrichment Analysis，基因集富集分析），可以评估整个基因表达谱，以确定一个基因集的显著性排名。在这种情况下，一组基因可能是与特定的生化途径相关的基因集合，也可能是与特定的生理功能、疾病过程或药理反应相关的基因集合。GSEA识别的是基因集合的行为。

        GSEA官网有官方推荐使用的基于JAVA的软件（下面称之为Broad GSEA 软件）。当然，GSEA在R语言中也是可以进行的。R语言实现的是 Broad GSEA 软件的 Pre‑Ranked GSEA 模式（后续会提及），等效于 gene‑set permutation 方法。

置换类型	英文	用途	适用场景
表型/样本置换	phenotype permutation	在样本标签之间随机置换，用来构建富集得分的随机分布	样本数量足够多（一般 ≥ 7 个/组）
基因集置换	gene set permutation	在排序向量内部随机打乱基因的顺序，用于样本数少、无法稳定表型置换时	小样本（< 7 个每组）或只有一个对比结果时

        因此，在大样本分析下两种方法结果的p 值和 FDR 会稍有差异。

一、准备数据

1. 准备好你要跑的基因集

        MSigDB 是一个广泛使用的基因集合注释数据库，它包含了大量关于基因集的注释信息，这些信息可以用于各种基因表达分析，尤其是在癌症生物学、免疫学和其他基因组学研究领域。

• https://www.gsea-msigdb.org/gsea/msigdb/index.jsp

        如上图，官网中有整合人和鼠的基因集，此外我们也可以按一定逻辑自行构建基因集（基因集本质就是若干个“通路+对应的基因”的集合）。根据自己的分析需求选择需要的基因集即可。

           Tips: 人源还是鼠源注意要分清楚，请选择对应物种的基因集进行分析！

2. 准备好基因排序（gene_list）

        GSEA输入的gene_list是一个将差异分析结果转化为基因排序的向量。GSEA 最理想的数据是全基因排序向量，过滤掉低表达噪音即可。
        对于 DESeq2 分析 （利用R语言的DESeq2包进行RNA-seq数据处理），最推荐使用 stat 列作为排序依据，因为它同时考虑了变化方向和显著程度。

         Broad GSEA 软件中也可以利用DESeq2获得的gene_list进行分析， 等效于Broad GSEA 软件的 Pre‑Ranked GSEA 模式，大家可以试试：

二、R语言操作

分析目标

        利用 DESeq2 差异分析得到的有序基因列表（gene_list），在 MSigDB 的 H (Hallmark) 基因集中进行 GSEA 富集分析，并可视化结果。

1. 加载R包

# 加载所有包，可利用BiocManager下载
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(patchwork)
library(dplyr)

GSEA 分析依赖几个常用包：

• clusterProfiler：核心 GSEA 分析；
• enrichplot：绘制 GSEA 曲线；
• ggplot2 / patchwork：美化和拼接图像。

2. 导入基因集（以H集为例）

        读取本地的基因集或你自定义的集：

# 读取本地的基因集文件，返回 TERM2GENE 格式
msig_h <- read.gmt("E:/R/h.all.v2025.1.Hs.symbols.txt") #改成你的集

3. GSEA富集分析

# 以 DESeq2 结果 res 为例，构建 gene_list
gene_list <- res$stat
names(gene_list) <- rownames(res)
gene_list <- gene_list[!is.na(gene_list)]
gene_list <- sort(gene_list, decreasing = TRUE)

set.seed(123)  # 保证结果可复现

# 运行 GSEA（等价于 Broad GSEA 的 Pre-ranked 模式）
gsea_h <- GSEA(
  geneList      = gene_list,
  TERM2GENE     = msig_h,
  minGSSize     = 10,
  pvalueCutoff  = 1,        # 不提前过滤 p 值
  pAdjustMethod = "fdr",
  nPerm         = 1000,
  verbose       = FALSE,
  seed          = TRUE
)

4. 保存结果

gsea_df_h <- as.data.frame(gsea_h)

write.csv(
  gsea_df_h,
  file = "E:/R/GSEA_Hallmark_Result.csv",
  row.names = FALSE
)

       把分析结果转为数据框格式并导出 CSV，包含了NES，p值，q值等信息，方便后续查看或在 Excel 里筛选。

5. 绘制GSEA富集图

        定义一个函数方便绘制任意通路的 GSEA 曲线，并保存为 PNG。这是本人目前使用下来最美观且兼容的绘制函数，可直接复制。

plot_gsea_enrichment <- function(gsea_result, pathway, title, filename) {
  # 绘制并组合上/下两部分图
  p <- gseaplot2(
    gsea_result,
    geneSetID = pathway,
    title     = title,
    color     = "darkgreen",
    base_size = 10,
    subplots  = 1:2,
    ES_geom   = "line"
  )

  p[[1]] <- p[[1]] + 
    theme(panel.grid.major.y = element_line(color = "gray90", size = 0.4))

  p_combined <- p[[1]] + p[[2]] +
    patchwork::plot_layout(ncol = 1, heights = c(3, 2)) +
    patchwork::plot_annotation(title = title)

  # 保存
  ggsave(
    filename  = filename,
    plot      = p_combined,
    device    = "png",
    width     = 6,
    height    = 4,
    dpi       = 600
  )
}

        示例：绘制 Hallmark 基因集中的 TNFA_SIGNALING 通路：

target_pathway <- "HALLMARK_TNFA_SIGNALING_VIA_NFKB" 

plot_gsea_enrichment(
  gsea_result = gsea_h,
  pathway     = target_pathway,
  title       = "Hallmark: TNFA Signaling via NFKB",
  filename    = "E:/R/GSEA_TNFA_SIGNALING.png"
)

Tips: 所需要绘制的通路必须存在于gsea结果文件中。

参考文献

Xu S, Hu E, Cai Y, et al. Using clusterProfiler to characterize multiomics data. Nat Protoc. 2024;19(11):3292-3320.