2026版单细胞测序scRNA-seq入门实战（十五）：irGSEA基因集富集分析——拒绝单一算法偏好，多算法联合打分

本节概览：

irGSEA简介：了解什么是irGSEA，以及它为什么要“集成”多种算法。
核心优势：一键调用6种主流算法 + 稳健秩聚合（RRA），消除单一算法偏好。
安装配置：irGSEA及其依赖包的安装。
数据准备：从Seurat对象出发，准备基因集与表达数据。
核心分析：使用irGSEA.score()一键完成多算法打分。
结果可视化：气泡图、山脊图、热图等多维度展示富集结果。
结果解读：从多算法一致性中锁定真实生物学信号。

传统GSEA是针对样本间比较设计的——它需要两个分组（如处理组vs对照组）来计算差异。

然而，在单细胞数据分析中，我们面临的是一个不同的问题：我们想知道某个基因集（如“炎症反应通路”、“缺氧标志基因”）在哪些细胞亚群中富集，而不是在两个预定义的组之间比较。

为此，单细胞领域发展出了多种基于排名的单样本基因集打分方法，如AUCell、UCell、ssGSEA、singscore等。它们只依赖于单个细胞内部的基因表达排名，不需要样本分组。

但问题来了：这么多方法，我该选哪个？

不同的算法在原理上各有侧重，对同一数据集的打分结果可能不一致甚至相互矛盾。单一算法的偏好性可能导致假阳性或遗漏真实信号。

irGSEA（integrated ranking Gene Set Enrichment Analysis）正是为解决这一问题而生的。它集成了6种主流算法，并通过稳健秩聚合（Robust Rank Aggregation, RRA） 整合多算法结果，帮助你从多算法的一致性中锁定真正的生物学信号。

一、irGSEA的核心概念

1.1 什么是“基于排名的单样本基因集打分”？

与需要两个分组比较的传统GSEA不同，基于排名的方法对每个细胞独立计算：

将该细胞中的所有基因按表达量从高到低排序
检查目标基因集（如“凋亡通路”）的基因在排序列表中的位置
如果这些基因集中在前列（高表达区域），则说明该基因集在该细胞中富集

这种方法的优势在于：不依赖样本分组、对数据规模和组成变化具有鲁棒性。

1.2 irGSEA集成了哪6种方法？

方法	核心原理	特点
AUCell	计算基因集在表达排序中的曲线下面积	最常用，稳健性好
UCell	AUCell的改进版，计算效率更高	适合大规模数据
singscore	基于基因排名的简单评分	计算极快，易于解释
ssGSEA	单样本GSEA的扩展	经典方法，广泛认可
JASMINE	基于排名的富集分析方法	适用于稀疏数据
Viper	基于蛋白质活性推断的富集分析	可整合调控网络信息

这些方法的共同特点是：仅依赖于单个细胞的相对基因表达水平，而非绝对值。

1.3 稳健秩聚合（RRA）：irGSEA的“投票机制”

irGSEA最核心的创新在于集成策略：

分别用6种方法对每个细胞计算基因集富集分数
对每种方法的结果进行归一化
使用稳健秩聚合（RRA）算法整合6套结果
输出一个综合富集分数和统计显著性

二、准备工作：安装irGSEA

2.1 安装依赖包

irGSEA的依赖包较多，但大部分可以通过CRAN和Bioconductor一键安装：

cran.packages <- c("aplot",  "doParallel", "doRNG",  "ggfun", "gghalves", 
                   "ggplotify", "ggridges", "ggsci", "irlba", "magrittr", "pagoda2", "plyr", "pointr", "purrr", 
                   "RcppML", "readr", "reticulate", "rlang", "RMTstat", 
                   "RobustRankAggreg", "roxygen2", "stringr", "tibble", "tidyr", "tidyselect", "tidytree", "VAM")

for (i in cran.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    install.packages(i, ask = FALSE, update = FALSE)
  }
}
bioconductor.packages <- c("fgsea", "ggtree", "GSEABase", "Nebulosa", "scde", 
                           "singscore", "SummarizedExperiment", "sparseMatrixStats")

for (i in bioconductor.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    BiocManager::install(i, ask = FALSE, update = FALSE)
  }
}

remotes::install_github("erocoar/gghalves")
pak::pkg_install("chuiqin/irGSEA")

⏱️ 安装提示：依赖包较多，安装过程可能需要10-20分钟，请耐心等待。

2.2 加载R包

# 加载R包
rm(list = ls())
library(Seurat)
library(tidyverse)
library(irGSEA)
library(msigdbr)        # 用于获取MSigDB基因集
library(ComplexHeatmap)
library(ggplot2)
library(patchwork)

set.seed(12345)

# 读取已注释的Seurat对象
seurat_obj <- readRDS("4_scRNA_celltype.rds")
dir.create("irGSEA", showWarnings = FALSE)
setwd("irGSEA")

# 查看细胞类型
table(seurat_obj$cell_type)

三、核心分析：irGSEA一键打分

3.1 运行irGSEA.score()

这是irGSEA最核心的函数——一键完成6种算法的打分：

# 运行irGSEA多算法联合打分
# 注意：irGSEA.score() 会自动识别基因集并进行多算法打分
# 对于大规模数据，建议先对细胞进行子集化（如随机抽取500-1000个细胞）
# 或按细胞类型分组分析

# 设定子集细胞数（例如 500 或 1000）
n_cells <- 500

# 按细胞类型分层抽样
cells_by_type <- split(colnames(seurat_obj), seurat_obj$cell_type)
cells_sub <- unlist(lapply(cells_by_type, function(x) {
  sample(x, size = min(length(x), round(n_cells * length(x) / ncol(seurat_obj))))
}))

seurat_sub <- seurat_obj[, cells_sub]
table(seurat_sub$cell_type)  # 检查比例
seurat_obj<-seurat_sub
seurat_obj <- irGSEA.score(object = seurat_obj, assay = "RNA",
                           slot = "data", seeds = 123, 
                           ncores = 1,
                           custom = F, geneset = NULL, msigdb = T,
                           species = "Mus musculus",
                           category = "H",
                           subcategory = NULL, 
                           geneid = "symbol",
                           min.cells = 3, min.feature = 0,
                           method = c("AUCell","UCell","singscore",
                                      "ssgsea", "JASMINE", "viper"),
                           aucell.MaxRank = NULL, 
                           ucell.MaxRank = NULL,
                           kcdf = 'Gaussian')
Assays(seurat_obj)

参数说明：

参数	说明
`method`	选择要运行的算法，可选 `"AUCell"`, `"UCell"`, `"singscore"`, `"ssgsea"`, `"JASMINE"`, `"Viper"`
`min.cells`	基因至少在多少个细胞中表达才纳入分析
`ncores`	并行核心数（Windows建议设为1）

⏱️ 运行时间：取决于细胞数和基因集数量。

3.2 整合差异基因集

result.dge <- irGSEA.integrate(object = sc_dataset2,
                               group.by = "celltype",
                               method = c("AUCell","UCell","singscore",
                                          "ssgsea", "JASMINE", "viper"))
class(result.dge)

四、结果可视化

4.1 气泡图：各细胞类型中富集的基因集

这是irGSEA最具代表性的可视化方式，展示哪些基因集在哪些细胞类型中显著富集：

# 气泡图展示基因集在细胞类型中的富集情况

# irGSEA提供了专门的气泡图函数
# 需要指定分组信息（如细胞类型）和感兴趣的方法
p_bubble <- irGSEA.bubble(
  object = seurat_obj,
  group.by = "cell_type",          # 分组列名
  method = "AUCell",               # 使用哪种方法的结果
  threshold = 0.05,                # 显著性阈值
  top = 15                         # 显示Top基因集数量
)

# 保存
ggsave("irGSEA_bubble_AUCell.pdf", p_bubble, width = 12, height = 8)

解读：

横轴：细胞类型
纵轴：基因集名称

4.2 山脊图（Ridge Plot）：展示基因集在细胞间的分布

# 山脊图展示特定基因集的活性分布
# 选择感兴趣的基因集
target_geneset <- "HALLMARK−BILE−ACID−METABOLISM"

# 绘制山脊图
p_ridge <- irGSEA.ridgeplot(object = pbmc3k.final,
                              method = "UCell",
                              show.geneset = target_geneset)

ggsave(paste0("ridge_", target_geneset, ".pdf"), p_ridge, width = 8, height = 6)

4.3 热图：基因集 × 细胞类型活性矩阵

p_heatmap <- irGSEA.heatmap(object = result.dge,
                                      method = "RRA",
                                      top = 50,
                                      show.geneset = NULL)
ggsave("irGSEA_heatmap.pdf", p_heatmap, width = 12, height = 8)

4.4 upset_plot

p_upset <- irGSEA.upset(object = result.dge,
                                  method = "RRA")
ggsave("irGSEA_upset.pdf", p_upset, width = 12, height = 8)

4.5 barplot

p_barplot <- irGSEA.barplot(object = result.dge,
                            method = c("AUCell","UCell","singscore",
                                       "ssgsea", "JASMINE", "viper"))
ggsave("irGSEA_barplot.pdf", p_barplot, width = 15, height = 8)

4.6 vlnplot

p_vlnplot <- irGSEA.halfvlnplot(
  object = seurat_obj,                      # 传入 Seurat 对象
  method = "UCell",                         # 方法名
  show.geneset = "HALLMARK-ADIPOGENESIS",
  group.by = "cell_type"                    # 按细胞类型分组
)
ggsave("irGSEA_vlnplot.pdf", p_vlnplot, width = 12, height = 6)

4.7 密度热图

p_densityheatmap <- irGSEA.densityheatmap(
  object = seurat_obj,                      # 使用 Seurat 对象
  method = "UCell",                         # 方法名
  show.geneset = "HALLMARK-ADIPOGENESIS",   # 基因集名称
  group.by = "cell_type"                    # 按细胞类型分组
)
ggsave("irGSEA_densityheatmap.pdf", p_densityheatmap, width = 8, height = 8)

4.7 密度散点图

p_scatterplot <- irGSEA.density.scatterplot(object =seurat_obj  ,
                             method = "UCell",
                             show.geneset = "HALLMARK-ADIPOGENESIS",
                             reduction = "umap")
ggsave("irGSEA_scatterplot.pdf", p_scatterplot, width = 8, height = 8)

五、稳健秩聚合（RRA）整合分析

irGSEA最强大的功能是使用稳健秩聚合（RRA）算法整合多方法结果：

# ============================================================
# RRA整合分析：从多算法一致性中锁定真实信号
# ============================================================

# irGSEA内置了RRA整合功能
# 它会综合所有方法的结果，输出一个整合的富集排名

# 提取整合结果
integrated_results <- seurat_obj@misc$irGSEA$integrated

# 查看整合后的基因集排名
head(integrated_results)

# 可视化整合结果
# 提取各细胞类型中Top基因集的平均秩
integrated_rank <- integrated_results %>%
  group_by(cell_type, geneset) %>%
  summarise(mean_rank = mean(rank, na.rm = TRUE)) %>%
  arrange(cell_type, mean_rank)

# 绘制Top基因集在各细胞类型中的排名热图
rank_matrix <- integrated_rank %>%
  spread(cell_type, mean_rank) %>%
  column_to_rownames("geneset")

# 取每个细胞类型中排名前10的基因集
top_genesets <- integrated_rank %>%
  group_by(cell_type) %>%
  slice_min(mean_rank, n = 10) %>%
  pull(geneset) %>%
  unique()

rank_matrix_top <- rank_matrix[top_genesets, ]

pdf("RRA_integrated_heatmap.pdf", width = 10, height = 12)
Heatmap(
  rank_matrix_top,
  name = "Mean Rank",
  col = colorRamp2(c(min(rank_matrix_top, na.rm = TRUE), 
                     median(rank_matrix_top, na.rm = TRUE), 
                     max(rank_matrix_top, na.rm = TRUE)), 
                   c("red", "white", "blue")),
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  show_row_names = TRUE,
  show_column_names = TRUE,
  row_names_gp = gpar(fontsize = 7)
)
dev.off()

六、结果解读与生物学意义

6.1 从气泡图看细胞类型功能

横轴扫描：每个细胞类型对应的气泡，颜色越红表示富集越显著
纵轴扫描：每个基因集对应的气泡分布，集中在某列说明该基因集在该细胞类型中特异富集
生物学解读：例如“炎症反应”基因集在小胶质细胞中富集，“突触传递”在神经元中富集——这与已知生物学知识一致，验证了分析的可靠性

6.2 从多算法一致性看信号可靠性

如果某个基因集在所有6种方法中都显示富集，说明该信号非常稳健
如果只在1-2种方法中显著，而在其他方法中不显著，可能是算法偏好性导致的假阳性
irGSEA的RRA整合结果自动帮你完成了这个“一致性筛选”

6.3 发现的基因集如何进一步分析？

回到差异基因：查看富集基因集中具体包含哪些基因
与GO/KEGG对比：irGSEA发现的富集基因集是否与之前GO/KEGG结果一致？
与文献对照：发现的基因集是否符合该细胞类型的已知功能？
新发现：如果某个基因集在特定细胞类型中富集但与已知功能不符，可能是新的生物学发现

七、常见问题与注意事项

7.1 irGSEA运行很慢怎么办？

减少细胞数：随机抽样500-1000个细胞进行探索性分析
减少基因集数量：先只分析Hallmark（50个）而非全部GO（数千个）
减少方法数量：先用2-3种方法（如AUCell + UCell + singscore），而非全部6种
增加并行核心：在Linux/Mac上设置ncores = 4或更高

7.2 不同方法结果差异很大怎么办？

这是正常的——不同算法的数学原理不同，对数据的敏感度也不同。这正是irGSEA设计的意义所在：通过RRA整合找出多方法一致的信号。如果某个基因集只在一种方法中显著，建议谨慎解读。

7.3 人类 vs 小鼠基因集

# 人类
msig_hallmark <- msigdbr(species = "Homo sapiens", category = "H")

# 小鼠
msig_hallmark <- msigdbr(species = "Mus musculus", category = "H")

⚠️ 务必根据你的数据物种选择正确的数据库！

7.4 irGSEA vs 传统GSEA

对比项	传统GSEA	irGSEA
适用场景	两组样本比较	单细胞亚群比较
输入	排序基因列表	表达矩阵 + 基因集
输出	ES/NES	每个细胞的富集分数
算法	单一算法	6种算法集成
分组依赖	需要分组	不需要

7.5 内存不足

对于大规模数据（>10000细胞），建议先按细胞类型子集化
或使用method = c("UCell", "singscore")（这两种方法内存占用较小）

八、总结

至此，我们完成了从Seurat对象到irGSEA多算法联合打分的完整流程：

步骤	函数/操作	输出
获取基因集	`msigdbr()` / 自定义	基因集列表
多算法打分	`irGSEA.score()`	6种方法的富集分数矩阵
气泡图	`irGSEA.bubble()`	基因集×细胞类型富集图
山脊图	`ggplot2` + `ggridges`	特定基因集的活性分布
UMAP	`FeaturePlot()`	基因集活性的空间分布
RRA整合	内置算法	多方法一致性排名

irGSEA的核心价值：

一站式解决：无需学习和安装6种不同的R包，一个函数搞定所有
消除算法偏好：通过多算法集成 + RRA整合，让真实生物学信号浮出水面
结果可视化丰富：气泡图、山脊图、热图、UMAP全覆盖
与Seurat无缝衔接：直接输入Seurat对象，输出也存入其中

在单细胞研究中，irGSEA填补了“如何评估基因集在细胞亚群中的富集” 这一分析空白。它既不像传统GSEA那样需要样本分组，也不像GO/KEGG那样只关注差异基因列表——它直接回答“这个基因集在哪些细胞中活跃” 这一核心生物学问题。

这里是两栖生物手册，中科院生物医学博士，持续记录生物医学实验和生物信息学笔记，干湿结合两不误～正努力成为最贴心负责的生信数据分析者和实验技术分享者，欢迎大家关注～