GSEA — 基因集富集分析（Gene Set Enrichment Analysis）¶

一句话说明¶

GSEA 不依赖任意 p 值截断，而是对全部基因按表达变化排序，检验预定义基因集（通路、功能模块）是否整体向排序列表的顶端或底端富集（当前 GSEA v4.3+）。

安装与配置¶

# 方法1：下载 GSEA 图形界面（Java 应用）
# 从官网下载：https://www.gsea-msigdb.org/gsea/downloads.jsp
wget https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp \
    -O GSEAPreranked.jar              # 命令行版本

# 方法2：使用 R 包 fgsea（更适合脚本化分析）
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("fgsea")        # 快速实现 GSEA 的 R 包

# 方法3：通过 clusterProfiler 包（推荐，下一章详细介绍）
BiocManager::install("clusterProfiler")

# 安装 MSigDB 数据包（提供 GSEA 基因集数据库）
install.packages("msigdbr")          # 包含 MSigDB 全部基因集

# 验证
library(fgsea)
library(msigdbr)
packageVersion("fgsea")

核心用法（R 包 fgsea 实现）¶

第一步：准备有序基因列表¶

library(fgsea)
library(msigdbr)
library(dplyr)

# 从 DESeq2 或 edgeR 结果准备有序基因列表
# 方法：用 -log10(pvalue) * sign(log2FC) 作排序统计量
deseq_res <- read.csv("DESeq2_results.csv", row.names=1)

# 创建命名向量：基因名 → 排序统计量（越大越上调）
gene_ranks <- deseq_res$stat           # DESeq2 的 Wald 统计量（推荐）
# 或手动计算：
# gene_ranks <- -log10(deseq_res$pvalue) * sign(deseq_res$log2FoldChange)

names(gene_ranks) <- rownames(deseq_res)  # 基因名作为名称

# 去除 NA 并排序
gene_ranks <- na.omit(gene_ranks)     # 去除 NA（pvalue 为 NA 的基因）
gene_ranks <- sort(gene_ranks, decreasing=TRUE)  # 降序排列

cat("总基因数:", length(gene_ranks), "\n")

第二步：获取基因集数据库¶

# 从 msigdbr 获取 MSigDB 基因集
# 常用基因集类别：
# H  — Hallmark（50个典型生物学状态，最常用）
# C2 — 经典通路（KEGG/Reactome/WikiPathways）
# C5 — GO terms（BP/MF/CC）
# C7 — 免疫相关

# 获取 Hallmark 基因集（人类）
hallmark_gene_sets <- msigdbr(species = "Homo sapiens",
                               category = "H")    # Hallmark

# 转换为 fgsea 需要的列表格式（基因集名 → 基因向量）
pathways <- split(hallmark_gene_sets$gene_symbol,  # 按基因集名分割
                  hallmark_gene_sets$gs_name)       # 基因符号作为基因 ID

cat("基因集数量:", length(pathways), "\n")

# 获取 KEGG 通路（C2 类别）
kegg_gene_sets <- msigdbr(species="Homo sapiens",
                           category="C2",
                           subcategory="CP:KEGG")
kegg_pathways <- split(kegg_gene_sets$gene_symbol,
                        kegg_gene_sets$gs_name)

第三步：运行 GSEA 分析¶

# 运行 fgsea（快速 GSEA 实现）
# pathways：基因集列表
# stats：有序基因列表（命名数值向量）
# minSize/maxSize：基因集大小过滤
# nPerm/nPermSimple：置换次数（越多越精确，越慢）
set.seed(42)                           # 设置随机种子保证可重现
gsea_result <- fgsea(
    pathways = pathways,               # 基因集列表
    stats = gene_ranks,                # 排序统计量
    minSize = 15,                      # 最小基因集大小（过滤太小的基因集）
    maxSize = 500,                     # 最大基因集大小（过滤太大的基因集）
    nPermSimple = 10000                # 简单置换次数（快速）
)

# 按 NES（标准化富集得分）排序
gsea_result <- gsea_result[order(gsea_result$NES, decreasing=TRUE), ]

# 筛选显著结果（padj < 0.25 是 GSEA 惯例）
sig_pathways <- subset(gsea_result, padj < 0.25)
cat("显著富集通路数:", nrow(sig_pathways), "\n")

# 保存结果（leadingEdge 列需要特殊处理）
gsea_export <- gsea_result
gsea_export$leadingEdge <- sapply(gsea_result$leadingEdge,
                                   paste, collapse=",")  # 将列表转为字符串
write.csv(gsea_export, "GSEA_results.csv", row.names=FALSE)

参数详解¶

参数	说明	默认值
`minSize`	基因集最小基因数	1
`maxSize`	基因集最大基因数	Inf
`nPermSimple`	简单置换次数	1000
`nproc`	并行线程数	0（全部）
`eps`	极小 p 值的下界	1e-10

可视化¶

# 1. 富集图（单个通路）
plotEnrichment(pathways[["HALLMARK_GLYCOLYSIS"]],  # 指定通路名
               gene_ranks) +                        # 排序统计量
    ggplot2::labs(title = "HALLMARK_GLYCOLYSIS")    # 标题

# 2. 气泡图（多通路汇总）
library(ggplot2)
top_pathways <- head(sig_pathways[order(sig_pathways$padj),], 20)

ggplot(top_pathways, aes(x=NES, y=reorder(pathway, NES),
                          size=size, color=padj)) +
    geom_point() +                                  # 气泡图
    scale_color_gradient(low="red", high="blue") +  # 颜色按 padj
    labs(x="NES", y="Pathway", title="GSEA Results") +
    theme_bw()

# 3. 条形图（上调/下调通路分开）
ggplot(top_pathways, aes(x=NES, y=reorder(pathway, NES),
                          fill=NES>0)) +
    geom_bar(stat="identity") +
    scale_fill_manual(values=c("FALSE"="blue","TRUE"="red")) +
    labs(x="NES", y="Pathway") + theme_bw()

常见报错与解决¶

报错信息	原因	解决方法
`No genes in pathways`	基因 ID 格式不匹配	统一用 gene symbol 或 Entrez ID
`All p-values are 1`	置换次数太少	增大 `nPermSimple`
`Infinite values in stats`	统计量有 Inf	用 `replace(x, is.infinite(x), max(x[is.finite(x)]))`
`pathway not found`	通路名拼写错误	`grep("GLYCOLYSIS", names(pathways))` 搜索正确名
结果 NES 极端	基因集与数据相关性强	正常现象，检查 leading edge 基因

速查表¶

# 完整最简流程
library(fgsea); library(msigdbr)
# 准备基因排序向量（来自 DESeq2 结果）
ranks <- setNames(deseq_res$stat, rownames(deseq_res))
ranks <- sort(na.omit(ranks), decreasing=TRUE)
# 获取基因集
pathways <- split(msigdbr(species="Homo sapiens", category="H")$gene_symbol,
                  msigdbr(species="Homo sapiens", category="H")$gs_name)
# 运行 GSEA
res <- fgsea(pathways, ranks, minSize=15, maxSize=500, nPermSimple=10000)
write.csv(res[order(res$padj),], "gsea_results.csv", row.names=FALSE)