1. 项目概述：从单细胞数据到生物学洞见的桥梁

如果你正在处理单细胞转录组数据，用Scanpy做完差异表达分析，拿到一长串差异基因列表后，是不是经常有种“老虎吃天，无从下口”的感觉？几百上千个基因名字摆在那里，它们到底意味着什么生物学过程被激活或抑制了？这时候，富集分析就是你不可或缺的“翻译官”。这个项目要解决的，就是如何将Scanpy分析得到的“基因列表”这份“原材料”，通过gseapy这个强大的工具，烹饪成一份易于理解的“通路解读”报告。这不仅仅是跑个代码，更是连接高通量数据与具体生物学意义的关键一步。

Scanpy作为单细胞分析领域的瑞士军刀，其差异分析结果（比如sc.tl.rank_genes_groups的输出）为我们筛选出了在特定细胞群或条件下表达发生显著变化的基因。但这些基因符号本身是沉默的。富集分析的核心思想是“物以类聚”，它基于一个基本假设：功能相关的基因往往会协同变化。通过将我们的基因列表与已知的基因功能数据库（如GO、KEGG、Reactome、MSigDB）进行比对，我们可以找出哪些生物学通路、分子功能或细胞组分在我们的数据中呈现出统计学上的显著富集。gseapy（Gene Set Enrichment Analysis in Python）正是执行这一任务的利器，它封装了多种经典算法（如GSEA、ORA），并提供了友好的Python接口，能与Scanpy生态无缝衔接。

这个实战指南适合所有正在或即将使用Scanpy分析单细胞数据的研究者、生物信息学入门者以及对功能注释感兴趣的实验生物学家。你将学到的不只是几行代码，而是从数据到生物学故事的一整套可复现、可解释的分析流程。我们将避开那些只讲函数调用的浅显教程，深入到参数选择背后的逻辑、结果解读的陷阱以及提升分析效率的实战技巧中。

2. 核心思路与工具选型：为什么是gseapy？

面对一个基因列表，进行富集分析的路子有很多。在R语言生态里，clusterProfiler几乎是标配，功能强大且社区活跃。那为什么在Python环境里，我们要选择gseapy？这背后是一系列针对单细胞数据分析场景的针对性考量。

2.1 与Scanpy的无缝集成：工作流连贯性

单细胞分析工作流通常较长，从原始数据质控、归一化、降维聚类到差异分析，Scanpy提供了一条龙服务。如果在差异分析后，为了做富集分析而被迫切换到R环境，不仅需要数据格式的转换（增加出错风险），也打断了分析思维的连续性。gseapy作为纯Python包，可以直接读取Scanpy的AnnData对象中存储的差异基因结果，或者直接处理Python列表、Pandas DataFrame，保证了从数据处理到生物学解释都在同一个Jupyter Notebook或Python脚本中完成。这种连贯性对于构建可复现的分析管道至关重要。

2.2 算法完备性与数据库支持

gseapy并非一个功能单一的玩具包。它实现了最常用的两种富集分析策略：

1. 过代表分析（Over-Representation Analysis, ORA）：这是最直观的方法。它需要一个预先设定的“差异基因”列表（例如，logFC > 1 & p_val_adj < 0.05），然后检验这个列表中的基因在某个通路基因集中是否显著过多。方法简单粗暴，适用于有明显阈值的情况。
2. 基因集富集分析（Gene Set Enrichment Analysis, GSEA）：这是一种更精细、更强大的方法。它不需要预先设定阈值来筛选差异基因，而是利用所有基因的排序信息（例如，按logFC从大到小排序）。GSEA会检验一个通路中的基因是否倾向于集中在排序列表的顶部或底部。这种方法能发现那些基因表达变化幅度虽不大但协调一致的通路，避免了ORA因阈值选择而丢失信息的问题。

在数据库方面，gseapy内置了连接到MSigDB（分子特征数据库）的接口，这是目前最全面、最权威的基因集资源之一，包含Hallmark、C2（KEGG、Reactome等）、C5（GO）等多个精选集合。同时，它也支持用户自定义的基因集文件（.gmt格式）。这意味着，你既可以使用KEGG、GO这些经典通路，也可以使用针对特定疾病、细胞类型或实验条件定制的基因集，灵活性极高。

2.3 性能与易用性平衡

对于单细胞数据，我们常常需要对多个细胞簇（cluster）或对比组分别进行富集分析。gseapy的API设计允许进行批量操作，例如，可以一次性对多个基因列表进行ORA分析，并以结构化的DataFrame返回结果，方便后续整理和可视化。其输出结果直接是Pandas DataFrame，与Python数据科学生态（如Matplotlib, Seaborn, Plotly）的整合天衣无缝，制作发表级的图表非常方便。

注意：虽然R的clusterProfiler在基因集资源和某些高级功能上可能更丰富，但gseapy在满足单细胞数据分析核心需求（集成、批量、可视化）方面已经做得足够出色。对于绝大多数应用场景，gseapy是Python环境下的最优解，避免了跨语言调用的复杂度。

3. 实战准备：从Scanpy结果到gseapy输入

理论说再多，不如动手做。我们假设你已经用Scanpy完成了一个标准的分析流程，得到了聚类结果，并针对某个感兴趣的细胞簇（比如Cluster 0）进行了差异表达分析。现在，我们要把Scanpy的输出，转化成gseapy能“吃”下去的格式。

3.1 提取差异基因列表

差异分析后，Scanpy将结果存储在adata.uns[‘rank_genes_groups’]中。我们需要从中提取特定分组的基因名和统计量。这里有两种主流输入准备方式，对应gseapy的两种主要分析模式。

方式一：用于ORA的基因列表（Symbol列表）ORA需要一个明确的“差异基因”列表。通常我们根据调整后p值（p_val_adj）和log2折叠变化（logfoldchanges）来筛选。

import scanpy as sc import pandas as pd # 假设你已经有了包含差异分析结果的adata对象 # 将差异分析结果转换为便于操作的DataFrame dea_result = sc.get.rank_genes_groups_df(adata, group='0') # 提取cluster 0 vs rest的结果 # 设定阈值筛选差异表达基因 signif_genes_df = dea_result[(dea_result['pvals_adj'] < 0.05) & (dea_result['logfoldchanges'].abs() > 1)] # 获取基因符号列表 ora_gene_list = signif_genes_df['names'].tolist() print(f“筛选得到 {len(ora_gene_list)} 个差异表达基因用于ORA分析。”)

这段代码的核心是sc.get.rank_genes_groups_df函数，它把Scanpy内部存储的差异结果变成了一个规整的DataFrame，后续的筛选和提取就变得非常直观。阈值（p_val_adj < 0.05, |logFC| > 1）是常用起点，但并非金科玉律，需要根据数据实际情况（如测序深度、细胞数）调整。

方式二：用于GSEA的基因排序列表（带分数的DataFrame）GSEA需要所有被检测基因的排序信息。通常我们按log2折叠变化（logFC）降序排列，得到一个从上调最显著到下调最显著的基因列表。

# 为GSEA准备数据：包含基因名和排序指标（如logFC）的DataFrame # 我们使用完整的差异分析结果，并按logFC排序 gsea_data_df = dea_result[['names', 'logfoldchanges']].copy() gsea_data_df.columns = ['gene_name', 'logFC'] # 重命名列以符合gseapy习惯 # 按logFC降序排列 gsea_data_df = gsea_data_df.sort_values(by='logFC', ascending=False) # 注意：GSEA也可以使用其他统计量（如t值、p值衍生值）进行排序，logFC是最直观的之一。

这里的关键是提供一个包含两列的DataFrame：基因名和用于排序的数值型指标。排序决定了GSEA算法检验的方向。

3.2 安装与配置gseapy环境

gseapy可以通过pip直接安装。强烈建议在虚拟环境（如conda环境）中进行。

pip install gseapy

安装后，首次使用可能需要下载基因集数据库。gseapy提供了在线和离线两种方式。对于国内用户，网络连接MSigDB官网可能不稳定，提前下载好数据库文件是更稳妥的做法。

import gseapy as gp # 查看可用的内置基因集库 print(gp.get_library_name()) # 在线方式（需稳定网络）：直接使用库名，gseapy会自动下载 # 离线方式：提前从MSigDB官网（https://www.gsea-msigdb.org/gsea/msigdb）下载.gmt文件 # 例如，下载了 c2.cp.kegg.v2023.1.Hs.symbols.gmt gene_sets = ‘./path/to/your/c2.cp.kegg.v2023.1.Hs.symbols.gmt’

实操心得：对于常用数据库（如KEGG、GO），建议在项目开始前就下载好对应的.gmt文件。这不仅能避免每次分析时的网络延迟和潜在失败，也保证了分析环境的可复现性。将.gmt文件存放在项目目录的data/或resources/子文件夹下是个好习惯。

4. 核心分析执行：ORA与GSEA详解

万事俱备，只欠东风。接下来我们分别进行ORA和GSEA分析，并解读核心输出。

4.1 ORA分析实战与结果解读

我们使用之前准备好的ora_gene_list和KEGG基因集进行ORA分析。

# 执行ORA分析 ora_res = gp.enrichr(gene_list=ora_gene_list, gene_sets=[‘KEGG_2021_Human’], # 可以同时指定多个库，如 [‘KEGG_2021_Human’, ‘GO_Biological_Process_2021’] organism=‘Human’, # 物种，必须与基因集匹配 outdir=None, # 设为None则不生成输出文件，结果只保存在变量中 cutoff=0.05 # 显著性截断值，通常看‘Adjusted P-value’ ) # 获取结果DataFrame ora_results_df = ora_res.results # 查看显著富集的前10条通路 print(ora_results_df.head(10)[[‘Term’, ‘Overlap’, ‘P-value’, ‘Adjusted P-value’, ‘Odds Ratio’, ‘Combined Score’]])

gp.enrichr是gseapy中执行ORA分析的函数。这里有几个关键参数：

• gene_sets: 可以传入库名（在线）或本地.gmt文件路径。传入列表可以一次性分析多个数据库。
• organism: 至关重要！必须与你的基因标识符（通常是Gene Symbol）和基因集数据库的物种一致，否则匹配不上。
• cutoff: 用于筛选最终展示结果的阈值，基于Adjusted P-value（经过多重检验校正的p值）。

结果解读要点：

• Term: 富集到的通路名称。
• Overlap: 格式如“15/200”，表示你的基因列表中有15个基因属于该通路，而该通路总共有200个基因。这个比例是富集的基础。
• P-value&Adjusted P-value: 富集分析的显著性p值。一定要看调整后的p值（Adjusted P-value），它控制了假阳性率。通常认为Adj. P-value < 0.05是显著的。
• Odds Ratio（比值比）: 表示你的基因列表中基因属于该通路的几率与背景基因相比的倍数。OR > 1表示正富集（即该通路在你列表中过代表），数值越大富集程度越强。
• Combined Score: 一个综合了p值和OR值的评分，用于对富集结果进行排序，分数越高通常意味着该结果越可信、越显著。

4.2 GSEA分析实战与结果解读

接下来，使用为GSEA准备的数据框gsea_data_df进行分析。

# 执行GSEA分析 gsea_res = gp.gsea(data=gsea_data_df, # 包含基因名和排序指标的DataFrame gene_sets=‘./data/c2.cp.kegg.v7.5.1.symbols.gmt’, # 使用本地KEGG基因集 cls=None, # 对于单列表排序GSEA，cls设为None permutation_num=1000, # 置换检验次数，默认1000，增加次数更稳定但更慢 outdir=‘./gsea_results’, # 输出目录，gseapy会生成一系列报告和图表 method=‘signal_to_noise’, # 排名度量方法。对于我们的logFC数据，用‘signal_to_noise’或‘t_test’均可 permutation_type=‘gene_set’, # 置换类型，‘gene_set’是标准做法 seed=42, # 随机种子，保证结果可重复 ) # GSEA的结果对象更复杂，核心富集结果在.results中 gsea_results_df = gsea_res.results # 查看富集分数（ES）最正（上调）和最负（下调）的前几条通路 print(“Top positively enriched pathways:”) print(gsea_results_df.sort_values(by=‘NES’, ascending=False).head(5)[[‘Term’, ‘NES’, ‘NOM p-val’, ‘FDR q-val’]]) print(“\nTop negatively enriched pathways:”) print(gsea_results_df.sort_values(by=‘NES’, ascending=True).head(5)[[‘Term’, ‘NES’, ‘NOM p-val’, ‘FDR q-val’]])

GSEA的参数更多，核心在于：

• permutation_num: 用于计算p值的置换检验次数。1000次是常用起点，对于非常小的基因集或需要极高精度时，可以增加到10000次，但计算时间会线性增加。
• method: 排名度量。因为我们直接提供了logFC，所以选择‘signal_to_noise’是合适的。如果你提供的是包含表达矩阵和表型标签的完整数据，gseapy可以自己计算排名。
• permutation_type: ‘gene_set’是默认且最常用的，它通过打乱基因标签来构建零分布。

GSEA结果解读要点：

• NES（Normalized Enrichment Score，标准化富集分数）: 这是GSEA的核心结果。NES > 0表示该通路基因在排序列表的顶部（上调端）富集；NES < 0表示在底部（下调端）富集。NES的绝对值越大，富集程度越强。
• NOM p-val（Nominal p-value）: 置换检验得到的原始p值。
• FDR q-val（False Discovery Rate q-value）: 经过多重检验校正后的q值。这是判断通路是否显著的主要依据，通常要求FDR q-val < 0.25。注意，GSEA的FDR阈值通常比ORA的0.05宽松，这是由算法本身的特点决定的，0.25是GSEA原始论文推荐的常用阈值。
• Lead_Edge（Leading Edge）: 结果中还会包含一个“Leading Edge”列，它列出了对该通路富集分数贡献最大的核心基因。这是后续进行深入机制研究的关键线索。

注意事项：ORA和GSEA的结果可能不完全一致。这是正常的，因为它们回答的是略有不同的问题。ORA问：“我的差异基因列表里，哪些通路特别多？”GSEA问：“在我的所有基因排序中，哪些通路的基因倾向于聚集在顶部或底部？”通常，GSEA能发现更细微、更协调的变化，而ORA对强差异基因的响应更直接。建议两者结合看，互相佐证。

5. 结果可视化：让洞见一目了然

再好的数据，如果不能直观呈现，其影响力也会大打折扣。gseapy内置了实用的绘图函数，但我们也完全可以利用其结果DataFrame，用Seaborn或Matplotlib制作更定制化的图表。

5.1 富集分析标准图

gseapy的dotplot和barplot是快速查看结果的利器。

# 绘制ORA结果的点图 gp.dotplot(ora_res.results, column=‘Adjusted P-value’, # 颜色映射的列 x=‘Gene_set’, # 分组，这里我们只用了KEGG，所以显示为一点 size=10, # 点的大小，可以映射到‘Odds Ratio’或‘Count’ title=‘ORA Enrichment Analysis (KEGG)’, cmap=‘viridis_r’, # 颜色映射，_r表示反转 ofname=‘./figures/ora_dotplot.png’ # 保存文件 ) # 绘制GSEA结果的条形图，展示Top N通路 # 首先筛选显著通路 (FDR < 0.25) gsea_sig = gsea_results_df[gsea_results_df[‘FDR q-val’] < 0.25] # 取NES绝对值最大的前10条通路 top_pathways = gsea_sig.reindex(gsea_sig[‘NES’].abs().sort_values(ascending=False).index).head(10) import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize=(10, 8)) # 根据NES正负赋予不同颜色 colors = [‘firebrick’ if x > 0 else ‘navy’ for x in top_pathways[‘NES’]] sns.barplot(data=top_pathways, y=‘Term’, x=‘NES’, palette=colors) plt.axvline(0, color=‘k’, linestyle=‘-’, linewidth=0.5) plt.xlabel(‘Normalized Enrichment Score (NES)’) plt.title(‘Top 10 Significantly Enriched Pathways (GSEA)’) plt.tight_layout() plt.savefig(‘./figures/gsea_top10_bar.png’, dpi=300) plt.show()

点图能同时展示通路的显著性（颜色）和富集基因数量（点大小），信息密度高。条形图则能清晰对比不同通路的NES大小和方向。

5.2 GSEA富集图谱解读

GSEA最经典的输出是富集图谱（Enrichment Plot）。gseapy在运行时会为每个显著富集的通路自动生成该图。你可以在指定的输出目录（如./gsea_results）下找到名为KEGG_XXX等的子文件夹，里面的KEGG_XXX.png就是该通路的富集图谱。 这张图包含三部分：

1. 顶部：富集分数（ES）曲线。曲线在横轴（基因按排序列表排列）上行走，当遇到属于该通路的基因时向上走，否则向下走。曲线的最终峰值就是ES值。一个在顶部出现高峰的曲线（峰在左侧）表示该通路基因在排序列表顶部富集（上调）；在底部出现低谷（峰在右侧）表示在底部富集（下调）。
2. 中部：基因排序列表的“命中”条（黑色竖线），标记了属于该通路的基因在排序中的位置。
3. 底部：基因排序指标（如logFC）沿排序列表的分布热图或线条图。

解读时，要结合NES值、FDR q-val和图谱形态。一个典型的显著上调通路图谱，其ES曲线应在左侧快速攀升至一个高峰，并且“命中”条密集地集中在排序列表的最前端。

5.3 自定义高级可视化：通路网络与聚类

当富集到的通路很多时，它们之间可能存在功能重叠。我们可以通过通路相似性聚类来简化解读。这需要利用GO或KEGG通路的层级结构或基因重叠信息。

# 示例：基于通路间基因重叠度进行聚类可视化（需要scipy, scikit-learn） from sklearn.metrics.pairwise import pairwise_distances from scipy.cluster.hierarchy import linkage, dendrogram, fcluster import numpy as np # 选取显著ORA结果 sig_ora = ora_results_df[ora_results_df[‘Adjusted P-value’] < 0.05].head(20) # 这里需要一个函数来获取每条通路的基因列表，可能需要从原始enrichr结果或数据库中解析 # 假设我们有一个字典 path_genes，键为通路名，值为基因集合 # 计算Jaccard相似度矩阵 pathway_names = sig_ora[‘Term’].tolist() n = len(pathway_names) jac_matrix = np.zeros((n, n)) for i in range(n): for j in range(n): set_i = path_genes[pathway_names[i]] set_j = path_genes[pathway_names[j]] jac_matrix[i, j] = len(set_i & set_j) / len(set_i | set_j) if (set_i | set_j) else 0 # 转换为距离矩阵 dist_matrix = 1 - jac_matrix # 层次聚类 linkage_matrix = linkage(dist_matrix, method=‘average’) # 绘制树状图 plt.figure(figsize=(12, 8)) dendrogram(linkage_matrix, labels=pathway_names, orientation=‘left’, leaf_font_size=10) plt.title(‘Hierarchical Clustering of Enriched Pathways (based on gene overlap)’) plt.xlabel(‘Distance (1 - Jaccard Similarity)’) plt.tight_layout() plt.show()

这种可视化能帮你发现哪些通路是高度相关的，可能指向同一个核心生物学过程，从而在撰写报告时进行归纳合并，使故事线更清晰。

6. 避坑指南与高级技巧

在实际操作中，你会遇到各种预料之外的问题。下面是我从多次实战中总结出的常见“坑”和应对技巧。

6.1 基因标识符匹配失败

这是新手遇到最多的问题。症状是富集分析结果为空，或者富集到的通路极少。

• 问题根源：你的基因列表中的基因标识符（如TP53）与基因集数据库中的标识符不匹配。常见原因有：
  1. 物种错误：人类数据用了小鼠的基因集。
  2. 标识符类型错误：数据库使用Gene Symbol（TP53），而你的列表是Ensembl ID（ENSG00000141510）或Entrez ID（7157）。
  3. 基因符号过时：你使用的基因符号是旧版本，而数据库是最新的。
• 解决方案：
  1. 统一物种：确保organism参数与数据一致。
  2. 标识符转换：在进行分析前，使用专业的ID转换工具。推荐使用mygene包，它在Python中非常方便。

  import mygene mg = mygene.MyGeneInfo() # 假设你的基因列表是Ensembl ID ensembl_ids = [‘ENSG00000141510’, ‘ENSG00000169083’] # 批量查询，转换为Gene Symbol result = mg.querymany(ensembl_ids, scopes=‘ensembl.gene’, fields=‘symbol’, species=‘human’) # 提取转换后的Symbol symbol_list = [hit[‘symbol’] for hit in result if ‘symbol’ in hit]

  3. 检查并清洗列表：去除重复项、空值和无法识别的基因名。

6.2 背景基因集的选择

ORA分析需要一个“背景”基因集，即所有可能被考虑到的基因集合。默认情况下，gseapy的enrichr会使用该数据库定义的全部基因作为背景。但在单细胞分析中，这有时并不合适。

• 问题：单细胞RNA-seq并非检测所有基因，许多低表达或未检测到的基因不应被纳入背景。使用全基因组背景可能导致富集分析灵敏度下降或假阳性。
• 解决方案：使用检测到的基因集合作为背景。你可以从Scanpy的adata.var_names中获取所有在数据集中被检测到的基因。

  # 获取所有检测到的基因作为背景 background_genes = adata.var_names.tolist() # 在enrichr中，可以通过自定义基因集库的方式间接实现，但enrichr函数本身不直接接受背景参数。 # 一个更直接的方法是使用gseapy的`prerank`函数（类似GSEA）或使用其他支持自定义背景的ORA工具（如scipy.stats.fisher_exact手动计算）。 # 对于gseapy，更常见的做法是确保你的基因列表是从这个检测到的基因集合中筛选出来的，这本身已经隐含了背景信息。

  重要提示：严格来说，enrichr在线版本使用的是其预设背景。对于更精确的控制，可以考虑使用gseapy的enrich函数（如果支持）或转向R的clusterProfiler（其enricher函数支持自定义背景）。在Python中，如果背景问题影响重大，手动实现基于超几何检验的ORA也是一个选择。

6.3 结果太多或太少

• 结果太多（数百条显著通路）：这通常意味着差异基因筛选阈值太宽松（如只用了p值<0.05，没看logFC），导致输入基因列表过长、噪声大。收紧筛选条件（如p_val_adj < 0.01 & |logFC| > 1.5）。另外，在解读时不要只看p值，要结合Odds Ratio或Combined Score，关注富集程度强且生物学意义明确的通路。
• 结果太少（没有或只有几条显著通路）：
  1. 检查基因标识符匹配（见6.1）。
  2. 放宽差异基因筛选阈值，适当增加输入基因数量。ORA需要一定的基因数量才能有统计效力。
  3. 尝试GSEA。GSEA不依赖硬阈值，可能能发现ORA漏掉的、基因表达变化温和但一致的通路。
  4. 考虑使用更广泛的基因集数据库，如GO_Biological_Process_2021比KEGG_2021_Human包含的通路更多、更细。

6.4 提升分析与解读效率的技巧

1. 批量处理多个细胞簇：写一个循环，对每个感兴趣的细胞簇进行差异分析和富集分析，并将结果保存到字典或列表中，最后统一汇总比较。

   clusters_of_interest = [‘0’, ‘1’, ‘2’] enrichment_results = {} for cluster in clusters_of_interest: dea_df = sc.get.rank_genes_groups_df(adata, group=cluster) sig_genes = dea_df[(dea_df[‘pvals_adj’] < 0.05) & (dea_df[‘logfoldchanges’].abs() > 1)][‘names’].tolist() if len(sig_genes) > 5: # 避免基因数太少 ora_res = gp.enrichr(gene_list=sig_genes, gene_sets=[‘KEGG_2021_Human’], organism=‘Human’) enrichment_results[cluster] = ora_res.results

2. 结果自动化报告：使用Python的Jinja2或WeasyPrint库，将每个簇的Top富集通路、关键基因和图表自动整合成HTML或PDF报告，极大节省时间。
3. 生物学解读不是罗列结果：不要简单地把Top 10通路扔进文章。要归纳。例如，如果富集到的通路大量涉及“细胞周期”、“DNA复制”，那么该细胞簇可能处于活跃增殖状态；如果涉及“炎症反应”、“TNF信号通路”，则可能提示免疫激活。结合你研究的生物学背景，将通路归类，讲一个连贯的故事。
4. 利用Leading Edge Genes：对于GSEA显著的通路，仔细查看其Lead_Edge基因。这些是驱动该通路富集的核心基因。对它们进行额外的表达模式检查（如绘制热图），能为你的机制假设提供最直接的证据。

富集分析是单细胞数据分析从描述性统计迈向生物学解释的关键一跃。掌握Scanpy+gseapy这套组合拳，意味着你不仅能告诉别人“这些细胞不同”，还能深入地阐述“它们为什么不同，以及这种不同可能意味着什么”。记住，工具是死的，生物学问题是活的。始终带着你的研究问题去审视富集分析的结果，让数据为你讲述一个可信的生物学故事。