1. 从基因列表到生物学故事为什么单细胞分析离不开富集分析如果你也和我一样每天和单细胞转录组数据打交道那你肯定对Scanpy这个工具不陌生。它帮我们完成了从原始数据到细胞分群、差异基因分析这一整套标准流程。但每次拿到那一长串差异表达基因列表时我总会遇到一个更挠头的问题这些基因到底意味着什么它们共同指向了哪些生物学过程是细胞在应激、在增殖还是在分化这时候一个强大的工具——gseapy就成为了连接“基因列表”与“生物学解读”的关键桥梁。今天我就结合自己处理神经发育和肿瘤微环境数据的实战经验来聊聊如何用Scanpy结合gseapy把冷冰冰的基因ID变成有血有肉的生物学故事。简单来说富集分析Enrichment Analysis就是一套“翻译”工具。它把我们从单细胞数据中筛选出的、成百上千个基因符号比如TP53, EGFR, CDKN1A与已知的生物学知识数据库比如GO、KEGG、Reactome进行比对然后告诉我们“嘿你这堆基因里参与‘细胞周期调控’和‘DNA损伤应答’这两个通路的基因特别多这很可能就是你样本里正在发生的主要事件。” gseapy就是Python生态里做这件事的佼佼者它封装了多种经典算法如GSEA、ORA并且与Scanpy的数据结构AnnData能很好地配合。接下来的内容我会带你走通从Scanpy差异分析结果出发到用gseapy完成富集分析并生成可发表级图表的所有关键步骤并分享那些在官方文档里不会写的参数调优心得和避坑指南。2. 环境搭建与数据准备构建可复现的分析流水线工欲善其事必先利其器。一个稳定、版本匹配的环境是后续所有分析可复现的基础。很多人直接pip install了事但不同工具包版本间的隐性冲突常常是诡异错误的源头。2.1 创建独立的Conda环境与精准版本控制我的习惯是为每个重要项目创建独立的Conda环境。这不仅避免了包冲突也便于记录和复现。以下是我为单细胞富集分析项目准备的environment.yml文件核心内容name: sc_enrichment channels: - conda-forge - bioconda - defaults dependencies: - python3.9 - scanpy1.9.3 - anndata0.8.0 - gseapy1.0.4 - pandas1.4.0 - numpy1.22.0 - scipy1.8.0 - matplotlib3.5.2 - seaborn0.11.2 - ipykernel - jupyter这里有几个关键点。首先我固定了Python 3.9这是一个在生物信息学工具链中兼容性极好的版本。其次我明确指定了Scanpy和gseapy的版本。gseapy的1.0.x版本与其早期的0.10.x版本在部分API上略有不同固定版本能确保教程代码的准确性。通过conda env create -f environment.yml创建环境后别忘了在Jupyter中注册这个内核python -m ipykernel install --user --namesc_enrichment --display-name单细胞富集分析。2.2 从Scanpy的AnnData对象中提取差异基因列表假设我们已经用Scanpy完成了标准的预处理、降维、聚类和差异表达分析。我们通常会将差异分析结果存储在adata.uns或adata.varm中。这里以使用scanpy.tl.rank_genes_groups函数并针对特定细胞簇比如簇‘0’进行分析为例import scanpy as sc import pandas as pd # 假设 adata 是已经完成差异分析的AnnData对象 # 提取簇‘0’ vs 其余所有细胞的差异分析结果 de_result adata.uns[rank_genes_groups] # 将其转换为一个便于操作的DataFrame de_df pd.DataFrame({ names: de_result[names][0], scores: de_result[scores][0], logfoldchanges: de_result[logfoldchanges][0], pvals: de_result[pvals][0], pvals_adj: de_result[pvals_adj][0] }) # 查看前10个差异基因 print(de_df.head(10))现在de_df这个DataFrame就包含了基因名、统计分数、log2倍变化和p值等信息。这是gseapy所需的原材料。但直接使用前我们必须进行关键的数据清洗。2.3 基因标识符的清洗与映射避免“无效基因名”的坑这是实战中第一个高频踩坑点。Scanpy分析中使用的基因标识符gene identifier可能是Ensembl ID如ENSG00000141510也可能是Gene Symbol如TP53。而gseapy的基因集数据库如MSigDB默认使用人类的Gene Symbol。如果标识符不匹配富集分析将找不到大部分基因导致结果为空或不可靠。解决方案是进行统一的基因标识符映射。我强烈推荐使用mygene这个Python包进行自动化的ID转换它背后整合了多个权威数据库。import mygene # 假设我们的基因名是Ensembl ID存储在 de_df[names] 中 gene_list de_df[names].tolist() mg mygene.MyGeneInfo() # 批量查询将Ensembl ID转换为Gene Symbol # scopes: 指定输入ID的类型这里是ensembl.gene # fields: 指定要返回的字段这里需要symbol # species: 指定物种人类是‘human’ gene_info mg.querymany(gene_list, scopesensembl.gene, fieldssymbol, specieshuman) # 将查询结果整理成一个映射字典 id_to_symbol {} for g in gene_info: # 确保查询成功且返回了symbol字段 if symbol in g: id_to_symbol[g[query]] g[symbol] else: # 记录未映射成功的基因便于后续检查 print(fWarning: Failed to map {g[query]}) # 将映射后的Gene Symbol添加到DataFrame中 de_df[gene_symbol] de_df[names].map(id_to_symbol) # 删除未能成功映射的行这些基因无法用于后续分析 de_df_clean de_df.dropna(subset[gene_symbol]).copy()注意mygene的在线查询可能有速率限制。对于非常大的基因列表可以考虑先保存本地映射文件或使用g:Profiler、biomaRtR包的等效方法。确保映射后的Gene Symbol没有重复如果有需要根据logFC或p值进行去重保留最显著的一个。至此我们得到了一个干净的、包含Gene Symbol和其对应统计显著性如p值、logFC的基因列表de_df_clean。这是进行后续富集分析的“弹药”。3. gseapy核心分析实战ORA与GSEA方法的选择与执行gseapy主要支持两种富集分析方法过表达分析ORA和基因集富集分析GSEA。它们适用场景不同理解其区别是正确解读结果的前提。3.1 方法选择何时用ORA何时用GSEAORA (Over-Representation Analysis)核心思想给定一个“感兴趣”的基因列表通常是差异分析中满足某个阈值如p_val_adj 0.05 \|logFC\| 1的基因看这些基因在某个已知通路/基因集中是否“过表达”。输入一个二分化的基因列表是或不是差异基因。优点计算简单结果直观容易理解。缺点需要人为设定阈值阈值的选择具有主观性且会丢失阈值附近基因的信息忽略了基因表达变化的方向和幅度。适用场景当你有一个明确的、筛选后的差异基因列表时想快速看看这些基因富集在哪些通路上。GSEA (Gene Set Enrichment Analysis)核心思想不设阈值使用全部基因按差异程度如logFC或p值排序。检验一个预先定义的基因集是否倾向于集中在排序列表的顶部或底部。输入一个按某种度量如logFC排序的所有基因列表。优点无需硬性阈值利用了全部数据信息能检测到那些基因表达变化幅度不大但协调一致的微妙效应可以区分基因集在列表顶部上调或底部下调的富集。缺点计算更复杂对基因集的构成和质量更敏感。适用场景当你不想或无法设定一个明确的差异基因阈值时当你不仅关心是否富集还关心富集的方向上调/下调时当你怀疑生物学效应是由许多基因的微弱但一致的变化引起时。在单细胞分析中我个人的经验是两者结合使用。先用ORA对高置信度的差异基因做一个快速、直观的概览再用GSEA对全基因排序列表进行分析捕捉更全面的信号并验证ORA的结果。3.2 执行ORA分析参数详解与结果解读假设我们从清洗后的数据中筛选出了显著上调的基因列表。import gseapy as gp # 1. 定义感兴趣的基因列表显著上调的基因 (示例阈值) sig_up_genes de_df_clean[(de_df_clean[pvals_adj] 0.05) (de_df_clean[logfoldchanges] 1)][gene_symbol].tolist() print(fNumber of significant up-regulated genes: {len(sig_up_genes)}) # 2. 选择基因集数据库。gseapy内置了多个常用的是 # - KEGG_2021_Human: KEGG通路 # - GO_Biological_Process_2021: GO生物过程 # - Reactome_2022: Reactome通路 # - MSigDB_Hallmark_2020: MSigDB Hallmark基因集高度概括非常常用 gene_sets KEGG_2021_Human # 3. 执行ORA分析 ora_results gp.enrichr(gene_listsig_up_genes, gene_setsgene_sets, organismHuman, # 必须指定 outdirNone, # 不输出文件结果保存在变量中 cutoff0.05 # 富集结果的FDR/q-value阈值 ) # enrichr返回一个包含多个DataFrame的字典我们通常需要的是.results ora_df ora_results.results # 查看最显著的前10个富集通路 print(ora_df.sort_values(Adjusted P-value).head(10)[[Term, Overlap, P-value, Adjusted P-value, Odds Ratio]])关键参数与结果解读organism必须正确指定否则会找不到对应的基因集。cutoff输出结果的显著性阈值基于校正后的P值Adjusted P-value通常是FDR。结果列解读Term富集到的通路/基因集名称。Overlap格式为X/Y表示你的基因列表中有X个基因属于该通路该通路总共有Y个基因。P-value富集分析的原始p值。Adjusted P-value经过多重检验校正后的p值如FDR这是我们判断结果是否显著的主要依据通常要求0.05。Odds Ratio比值比大于1表示正富集你的基因在该通路中富集数值越大富集程度越强。Genes你的基因列表中属于该通路的基因具体是哪些。实操心得enrichr函数默认会同时查询多个数据库如果gene_sets是一个列表。但初次使用时我建议先专注于一个数据库如‘KEGG_2021_Human’避免结果过于庞杂。另外outdir参数可以指定路径这样函数会自动将结果保存为文本文件和图表非常方便。3.3 执行GSEA分析排序列表的构建与富集信号挖掘GSEA需要的是一个排好序的基因列表。通常我们使用log2FoldChange作为排序指标因为它包含了变化的方向和幅度。# 1. 为GSEA准备数据一个包含Gene Symbol和排序指标的DataFrame # 我们使用清洗后的全部基因按log2FoldChange从大到小排序正值越大排越前负值越小排越后 gsea_data de_df_clean[[gene_symbol, logfoldchanges]].copy() gsea_data gsea_data.set_index(gene_symbol) # GSEA需要的是一个Series索引是基因名值是排序指标 gsea_rank gsea_data[logfoldchanges].sort_values(ascendingFalse) # 2. 执行GSEA分析 # 这里使用MSigDB的Hallmark基因集它包含了50个精炼的、具有明确生物学状态的基因集 gsea_results gp.gsea(datagsea_rank, gene_setsMSigDB_Hallmark_2020, organismhuman, permutation_num1000, # 置换检验次数默认1000增加次数更稳定但更慢 outdir./gsea_output, # 指定输出目录会生成详细结果和图表 formatpng, seed42 # 设置随机种子保证结果可复现 ) # gsea函数返回一个GSEA结果对象其.res2d属性是主要的结果DataFrame gsea_df gsea_results.res2d # 查看富集最显著NES绝对值大FDR小的前10个基因集 # NES (Normalized Enrichment Score) 标准化富集分数是核心指标 significant_results gsea_df[(gsea_df[FDR q-val] 0.25) (abs(gsea_df[NES]) 1.0)] print(significant_results.sort_values(NES, ascendingFalse).head(10)[[Term, NES, FDR q-val, Lead_genes]])GSEA核心结果解读NES标准化富集分数这是GSEA的灵魂指标。正负号表示方向NES 0表示该基因集在排序列表的顶部富集即在你的差异分析中该基因集对应的基因普遍上调。NES 0表示在底部富集普遍下调。绝对值表示强度|NES|越大富集程度越强。通常认为|NES| 1.0有生物学意义。FDR q-val错误发现率。GSEA官方推荐使用FDR 0.25作为显著性阈值这比ORA的0.05宽松因为GSEA探测的是更微妙的协调性变化。当然在严格的研究中也可以使用0.05。Lead_genes核心贡献基因。即在该基因集的富集信号中贡献最大的那些基因位于排序列表最顶端或最底端且属于该基因集。它们是解释该富集结果的关键。ES富集分数未标准化的原始分数受基因集大小影响一般看NES即可。参数设置经验permutation_num置换检验次数。默认1000次对于大多数分析足够。如果你追求极致的p值精度或者基因集很小/很大可以增加到5000或10000次但计算时间会线性增加。seed务必设置。这能确保每次运行的结果完全一致对于可复现性至关重要。outdirgsea函数的一个巨大优点是会自动生成可视化结果图富集分析图Enrichment Plot这对于理解和展示结果非常有帮助。4. 高级结果可视化与生物学解读让数据自己说话得到富集分析表格只是第一步如何将结果高效、美观地呈现出来并挖掘其背后的生物学故事才是更见功力的地方。gseapy和matplotlib/seaborn的结合可以创造出强大的可视化效果。4.1 绘制经典富集分析气泡图/条形图气泡图是展示富集结果最直观的方式之一它同时展示了通路的显著性-log10(p-value)或FDR、富集强度Odds Ratio或NES和基因重叠数量。import matplotlib.pyplot as plt import seaborn as sns # 假设 ora_df 是我们之前ORA分析的结果DataFrame # 选取FDR最显著的前15个通路进行可视化 top_ora ora_df.sort_values(Adjusted P-value).head(15).copy() # 计算 -log10(Adjusted P-value) 用于作图 top_ora[-log10(FDR)] -np.log10(top_ora[Adjusted P-value]) # 从‘Overlap’列中提取重叠基因数 top_ora[Overlap_Count] top_ora[Overlap].apply(lambda x: int(x.split(/)[0])) plt.figure(figsize(10, 8)) # 绘制气泡图x轴为Odds Ratioy轴为通路名称气泡大小为重叠基因数颜色为 -log10(FDR) scatter plt.scatter(xtop_ora[Odds Ratio], ytop_ora[Term], stop_ora[Overlap_Count]*20, # 气泡大小乘以一个系数放大视觉效果 ctop_ora[-log10(FDR)], cmapviridis_r, # 使用颜色映射_r表示反转使得值越大颜色越深 alpha0.7, edgecolorsk, linewidth0.5) plt.axvline(x1, colorgrey, linestyle--, linewidth0.8) # 添加Odds Ratio1的参考线 plt.xlabel(Odds Ratio, fontsize12) plt.ylabel(Pathway/Term, fontsize12) plt.title(ORA Enrichment Analysis (Top 15 Pathways), fontsize14, pad20) # 添加颜色条 cbar plt.colorbar(scatter) cbar.set_label(-log10(FDR), fontsize11) # 添加图例表示气泡大小需要手动创建 import matplotlib.patches as mpatches # 创建几个示例大小的图例句柄 size_legend_handles [mpatches.Circle((0,0), radiuss, facecolorgray, edgecolork, alpha0.7) for s in [5, 10, 15]] # 半径对应 Overlap_Count * 20 后的视觉大小 size_labels [5, 10, 15] # 对应的重叠基因数 plt.legend(handlessize_legend_handles, labelssize_labels, titleGene Count, loccenter left, bbox_to_anchor(1.05, 0.5)) plt.tight_layout() plt.show()这张图可以一目了然地告诉我们哪些通路最显著颜色深、富集程度最强Odds Ratio远大于1、并且有足够多的基因支持气泡大。4.2 解读GSEA富集分析图gseapy在运行GSEA时生成的PNG图片位于指定的outdir中是标准输出。以“HALLMARK_INFLAMMATORY_RESPONSE”为例生成的图片包含三部分顶部主图Enrichment Score Plot中间的曲线是富集分数ES在排序基因列表上的运行轨迹。峰值出现在哪里就表示该基因集的成员基因在排序列表的那个位置集中出现。峰值在左侧排名靠前表示上调富集NES0在右侧表示下调富集NES0。图中的绿色竖线就是Lead_genes核心基因的位置。中间部分Gene List Metric显示了所有基因按排序指标这里是logFC的分布情况像一条山脉。这直观展示了你输入数据的排序质量。底部热图Heatmap显示了该基因集中每个基因在原始数据集如果你提供了表达矩阵或在你排序列表中的相对位置/表达模式。在单细胞语境下如果我们提供了所有细胞的表达矩阵这里可以显示出该基因集在哪些细胞簇中高表达。解读技巧不要只看FDR和NES的数字。一定要打开这些图观察ES曲线的形状。一个“漂亮”的富集信号通常表现为一个早期达到峰值或谷值的、单峰的、清晰的曲线。如果曲线来回震荡峰值不明显即使FDR勉强显著其生物学意义也可能值得怀疑。4.3 整合多组结果与自定义基因集分析有时我们需要比较不同细胞簇或不同实验条件之间的富集结果。这时可以绘制多组对比气泡图或热图。# 假设我们对三个细胞簇cluster0,1,2都做了ORA分析结果存储在字典中results_dict {cluster0: ora_df0, cluster1: ora_df1, cluster2: ora_df2} # 我们想找出那些在多个簇中都富集的通路 common_terms set() for cluster, df in results_dict.items(): # 获取每个簇中FDR0.05的显著通路 sig_terms set(df[df[Adjusted P-value] 0.05][Term]) if not common_terms: common_terms sig_terms else: common_terms common_terms.intersection(sig_terms) print(fPathways enriched in all three clusters: {common_terms}) # 为了可视化我们可以创建一个矩阵行是通路列是细胞簇值是 -log10(FDR) vis_data [] for term in list(common_terms)[:20]: # 取前20个共同通路示例 row {Term: term} for cluster, df in results_dict.items(): # 在对应簇的结果中查找该通路的FDR值 term_row df[df[Term] term] if not term_row.empty: row[cluster] -np.log10(term_row.iloc[0][Adjusted P-value]) else: row[cluster] 0 # 如果不显著赋值为0 vis_data.append(row) vis_df pd.DataFrame(vis_data).set_index(Term) # 绘制热图 plt.figure(figsize(8, 10)) sns.heatmap(vis_df, cmapYlOrRd, linewidths0.5, annotTrue, fmt.2f, cbar_kws{label: -log10(FDR)}) plt.title(Common Enriched Pathways Across Clusters, fontsize14, pad20) plt.tight_layout() plt.show()此外gseapy支持使用自定义基因集。你可以将自己关注的、来自文献或特定研究的基因列表保存为.gmt格式文件然后通过gene_sets‘/path/to/your_custom.gmt’参数调用。这在进行非常聚焦的、与特定生物学问题相关的研究时极其有用。5. 实战避坑与进阶技巧来自一线踩坑的经验纸上得来终觉浅绝知此事要躬行。下面分享几个我在实际项目中反复遇到并且耗费不少时间才解决的问题。5.1 物种匹配错误与基因同源映射问题分析小鼠单细胞数据时直接使用Gene Symbol进行富集分析结果一片空白或非常奇怪。根因gseapy的许多内置数据库如KEGG, GO默认是面向人类的。虽然它也支持小鼠‘Mouse’、大鼠等但基因标识符必须与数据库匹配。小鼠的基因符号如Trp53和人类的TP53不同。解决方案使用正确的物种参数在执行gp.enrichr或gp.gsea时务必设置organism‘Mouse’。gseapy会自动调用对应物种的数据库。进行严谨的基因同源映射如果你的物种不在gseapy直接支持的范围如斑马鱼、果蝇或者你想跨物种比较就需要进行同源基因映射。可以使用mygene的homologene功能或专业的同源数据库如Ensembl BioMart将你的基因列表映射到人类同源基因上再使用人类数据库进行分析。这是一个复杂但有时不可避免的步骤。# 使用mygene进行小鼠到人类的同源映射示例简化 mg mygene.MyGeneInfo() mouse_genes [Trp53, Myc, Sox2] # 查询小鼠基因的人类同源物 homologs mg.querymany(mouse_genes, scopessymbol, fieldshomologene.human_symbol, speciesmouse) human_gene_list [] for h in homologs: if homologene in h and human_symbol in h[homologene]: human_gene_list.append(h[homologene][human_symbol]) else: print(fNo human homolog found for {h.get(query)}) print(fMapped human genes: {human_gene_list})5.2 背景基因集的正确设置问题ORA分析结果中某些通路的Odds Ratio异常高比如几百或者p值极其显著但重叠基因数很少比如1/1结果不可信。根因默认情况下enrichr使用其服务器上对应物种的“所有基因”作为背景background。但在单细胞RNA-seq中由于技术限制我们只能检测到约10000-20000个基因远少于全基因组背景~20000个蛋白编码基因。使用不匹配的背景集会导致统计偏差。解决方案自定义背景基因集。你应该将本次单细胞实验中所有被检测到的基因即adata.var_names映射为Gene Symbol后的列表作为背景。# 获取本次分析中所有检测到的基因背景集 # 假设我们已经将adata的索引转换为Gene Symbol或有一个映射关系 # 这里假设 adata.var[gene_symbol] 存储了Gene Symbol background_genes adata.var[gene_symbol].dropna().unique().tolist() print(fBackground gene set size: {len(background_genes)}) # 在执行enrichr时通过 background 参数传入 ora_results_custom_bg gp.enrichr(gene_listsig_up_genes, gene_setsKEGG_2021_Human, organismHuman, backgroundbackground_genes, # 关键参数 outdirNone, cutoff0.05)使用自定义背景后Odds Ratio和p值的计算会更加准确能有效过滤掉因背景集过大而造成的假阳性富集信号。5.3 处理大规模分析与结果过滤问题一次分析几十个细胞簇每个簇都做一次富集分析产生海量结果难以管理和聚焦。解决方案自动化脚本 结构化结果过滤。自动化循环将分析过程封装成函数对每个簇进行循环分析。结构化存储将每个簇的显著结果如FDR0.05存储到一个统一的、结构化的DataFrame或字典中包含簇名、通路名、NES、FDR、重叠基因等关键信息。结果过滤与聚合基于FDR和NES阈值这是最基本的一步。去除过于宽泛的通路GO或KEGG中有些通路非常宽泛如“代谢过程”、“信号转导”虽然显著但信息量低。可以手动建立一个“宽泛通路黑名单”进行过滤或者通过通路的基因数量Term大小进行筛选例如只保留基因数在10到500之间的通路。聚类相似通路使用gseapy的plot模块中的dotplot函数或GOplot、clusterProfilerR等工具可以对富集到的通路进行语义相似性聚类将相似的通路归为一类使结果更简洁。# 示例循环分析多个簇并收集显著结果 all_significant_results [] for cluster in adata.obs[leiden].cat.categories: # 1. 提取该簇的差异基因 (假设已有差异分析结果) de_genes get_de_genes_for_cluster(adata, cluster) # 这是一个自定义函数示例 if len(de_genes) 5: # 如果差异基因太少跳过 continue # 2. 执行富集分析 ora_res gp.enrichr(gene_listde_genes, gene_setsGO_Biological_Process_2021, ...) ora_df ora_res.results # 3. 过滤显著结果 sig_df ora_df[ora_df[Adjusted P-value] 0.05].copy() sig_df[cluster] cluster # 标记来源簇 # 4. 收集 all_significant_results.append(sig_df[[cluster, Term, Adjusted P-value, Odds Ratio, Overlap]]) # 合并所有结果 final_results_df pd.concat(all_significant_results, ignore_indexTrue) # 后续可以方便地用pandas进行筛选、排序和可视化5.4 与Scanpy生态的深度整合将通路分数嵌入AnnData富集分析的结果不仅可以用来画图和写文章还可以反向嵌入到Scanpy的AnnData对象中为后续的细胞层面分析提供新维度。例如我们可以计算每个细胞在某个特定通路如“上皮-间质转化EMT”上的“活性分数”。一种常见的方法是使用单样本基因集富集分析ssGSEA或AUCell算法。虽然gseapy本身不直接提供与Scanpy的深度接口但我们可以利用其基因集和Scanpy的sc.tl.score_genes功能进行近似计算。# 假设我们从GSEA结果中找到了一个感兴趣的基因集例如‘HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION’ # 我们需要这个基因集的成员基因列表 emt_gene_set ... # 从gsea_results.results_on_cell或MSigDB官网获取该基因集的基因列表 # 使用Scanpy的score_genes函数计算每个细胞在这个基因集上的“富集分数” sc.tl.score_genes(adata, emt_gene_set, score_nameEMT_score, ctrl_size50, use_rawFalse) # 参数说明 # - adata: AnnData对象 # - emt_gene_set: 基因列表 # - score_name: 存储在adata.obs中的列名 # - ctrl_size: 用于标准化分数的随机背景基因集大小 # - use_raw: 是否使用adata.raw # 现在adata.obs中多了一列‘EMT_score’ # 我们可以可视化这个分数在UMAP图上的分布 sc.pl.umap(adata, color[EMT_score, leiden], cmapRdBu_r, vcenter0) # 或者看这个分数在不同细胞簇间的差异 sc.pl.violin(adata, [EMT_score], groupbyleiden)这样我们就将一个通路水平的富集信息转化为了每个细胞的一个连续型特征EMT_score可以用于后续的聚类验证、差异分析、甚至作为回归模型的输入极大地拓展了分析的深度和灵活性。这个技巧在分析肿瘤异质性、细胞分化轨迹等复杂生物学过程中尤其有用。
1. 从基因列表到生物学故事为什么单细胞分析离不开富集分析如果你也和我一样每天和单细胞转录组数据打交道那你肯定对Scanpy这个工具不陌生。它帮我们完成了从原始数据到细胞分群、差异基因分析这一整套标准流程。但每次拿到那一长串差异表达基因列表时我总会遇到一个更挠头的问题这些基因到底意味着什么它们共同指向了哪些生物学过程是细胞在应激、在增殖还是在分化这时候一个强大的工具——gseapy就成为了连接“基因列表”与“生物学解读”的关键桥梁。今天我就结合自己处理神经发育和肿瘤微环境数据的实战经验来聊聊如何用Scanpy结合gseapy把冷冰冰的基因ID变成有血有肉的生物学故事。简单来说富集分析Enrichment Analysis就是一套“翻译”工具。它把我们从单细胞数据中筛选出的、成百上千个基因符号比如TP53, EGFR, CDKN1A与已知的生物学知识数据库比如GO、KEGG、Reactome进行比对然后告诉我们“嘿你这堆基因里参与‘细胞周期调控’和‘DNA损伤应答’这两个通路的基因特别多这很可能就是你样本里正在发生的主要事件。” gseapy就是Python生态里做这件事的佼佼者它封装了多种经典算法如GSEA、ORA并且与Scanpy的数据结构AnnData能很好地配合。接下来的内容我会带你走通从Scanpy差异分析结果出发到用gseapy完成富集分析并生成可发表级图表的所有关键步骤并分享那些在官方文档里不会写的参数调优心得和避坑指南。2. 环境搭建与数据准备构建可复现的分析流水线工欲善其事必先利其器。一个稳定、版本匹配的环境是后续所有分析可复现的基础。很多人直接pip install了事但不同工具包版本间的隐性冲突常常是诡异错误的源头。2.1 创建独立的Conda环境与精准版本控制我的习惯是为每个重要项目创建独立的Conda环境。这不仅避免了包冲突也便于记录和复现。以下是我为单细胞富集分析项目准备的environment.yml文件核心内容name: sc_enrichment channels: - conda-forge - bioconda - defaults dependencies: - python3.9 - scanpy1.9.3 - anndata0.8.0 - gseapy1.0.4 - pandas1.4.0 - numpy1.22.0 - scipy1.8.0 - matplotlib3.5.2 - seaborn0.11.2 - ipykernel - jupyter这里有几个关键点。首先我固定了Python 3.9这是一个在生物信息学工具链中兼容性极好的版本。其次我明确指定了Scanpy和gseapy的版本。gseapy的1.0.x版本与其早期的0.10.x版本在部分API上略有不同固定版本能确保教程代码的准确性。通过conda env create -f environment.yml创建环境后别忘了在Jupyter中注册这个内核python -m ipykernel install --user --namesc_enrichment --display-name单细胞富集分析。2.2 从Scanpy的AnnData对象中提取差异基因列表假设我们已经用Scanpy完成了标准的预处理、降维、聚类和差异表达分析。我们通常会将差异分析结果存储在adata.uns或adata.varm中。这里以使用scanpy.tl.rank_genes_groups函数并针对特定细胞簇比如簇‘0’进行分析为例import scanpy as sc import pandas as pd # 假设 adata 是已经完成差异分析的AnnData对象 # 提取簇‘0’ vs 其余所有细胞的差异分析结果 de_result adata.uns[rank_genes_groups] # 将其转换为一个便于操作的DataFrame de_df pd.DataFrame({ names: de_result[names][0], scores: de_result[scores][0], logfoldchanges: de_result[logfoldchanges][0], pvals: de_result[pvals][0], pvals_adj: de_result[pvals_adj][0] }) # 查看前10个差异基因 print(de_df.head(10))现在de_df这个DataFrame就包含了基因名、统计分数、log2倍变化和p值等信息。这是gseapy所需的原材料。但直接使用前我们必须进行关键的数据清洗。2.3 基因标识符的清洗与映射避免“无效基因名”的坑这是实战中第一个高频踩坑点。Scanpy分析中使用的基因标识符gene identifier可能是Ensembl ID如ENSG00000141510也可能是Gene Symbol如TP53。而gseapy的基因集数据库如MSigDB默认使用人类的Gene Symbol。如果标识符不匹配富集分析将找不到大部分基因导致结果为空或不可靠。解决方案是进行统一的基因标识符映射。我强烈推荐使用mygene这个Python包进行自动化的ID转换它背后整合了多个权威数据库。import mygene # 假设我们的基因名是Ensembl ID存储在 de_df[names] 中 gene_list de_df[names].tolist() mg mygene.MyGeneInfo() # 批量查询将Ensembl ID转换为Gene Symbol # scopes: 指定输入ID的类型这里是ensembl.gene # fields: 指定要返回的字段这里需要symbol # species: 指定物种人类是‘human’ gene_info mg.querymany(gene_list, scopesensembl.gene, fieldssymbol, specieshuman) # 将查询结果整理成一个映射字典 id_to_symbol {} for g in gene_info: # 确保查询成功且返回了symbol字段 if symbol in g: id_to_symbol[g[query]] g[symbol] else: # 记录未映射成功的基因便于后续检查 print(fWarning: Failed to map {g[query]}) # 将映射后的Gene Symbol添加到DataFrame中 de_df[gene_symbol] de_df[names].map(id_to_symbol) # 删除未能成功映射的行这些基因无法用于后续分析 de_df_clean de_df.dropna(subset[gene_symbol]).copy()注意mygene的在线查询可能有速率限制。对于非常大的基因列表可以考虑先保存本地映射文件或使用g:Profiler、biomaRtR包的等效方法。确保映射后的Gene Symbol没有重复如果有需要根据logFC或p值进行去重保留最显著的一个。至此我们得到了一个干净的、包含Gene Symbol和其对应统计显著性如p值、logFC的基因列表de_df_clean。这是进行后续富集分析的“弹药”。3. gseapy核心分析实战ORA与GSEA方法的选择与执行gseapy主要支持两种富集分析方法过表达分析ORA和基因集富集分析GSEA。它们适用场景不同理解其区别是正确解读结果的前提。3.1 方法选择何时用ORA何时用GSEAORA (Over-Representation Analysis)核心思想给定一个“感兴趣”的基因列表通常是差异分析中满足某个阈值如p_val_adj 0.05 \|logFC\| 1的基因看这些基因在某个已知通路/基因集中是否“过表达”。输入一个二分化的基因列表是或不是差异基因。优点计算简单结果直观容易理解。缺点需要人为设定阈值阈值的选择具有主观性且会丢失阈值附近基因的信息忽略了基因表达变化的方向和幅度。适用场景当你有一个明确的、筛选后的差异基因列表时想快速看看这些基因富集在哪些通路上。GSEA (Gene Set Enrichment Analysis)核心思想不设阈值使用全部基因按差异程度如logFC或p值排序。检验一个预先定义的基因集是否倾向于集中在排序列表的顶部或底部。输入一个按某种度量如logFC排序的所有基因列表。优点无需硬性阈值利用了全部数据信息能检测到那些基因表达变化幅度不大但协调一致的微妙效应可以区分基因集在列表顶部上调或底部下调的富集。缺点计算更复杂对基因集的构成和质量更敏感。适用场景当你不想或无法设定一个明确的差异基因阈值时当你不仅关心是否富集还关心富集的方向上调/下调时当你怀疑生物学效应是由许多基因的微弱但一致的变化引起时。在单细胞分析中我个人的经验是两者结合使用。先用ORA对高置信度的差异基因做一个快速、直观的概览再用GSEA对全基因排序列表进行分析捕捉更全面的信号并验证ORA的结果。3.2 执行ORA分析参数详解与结果解读假设我们从清洗后的数据中筛选出了显著上调的基因列表。import gseapy as gp # 1. 定义感兴趣的基因列表显著上调的基因 (示例阈值) sig_up_genes de_df_clean[(de_df_clean[pvals_adj] 0.05) (de_df_clean[logfoldchanges] 1)][gene_symbol].tolist() print(fNumber of significant up-regulated genes: {len(sig_up_genes)}) # 2. 选择基因集数据库。gseapy内置了多个常用的是 # - KEGG_2021_Human: KEGG通路 # - GO_Biological_Process_2021: GO生物过程 # - Reactome_2022: Reactome通路 # - MSigDB_Hallmark_2020: MSigDB Hallmark基因集高度概括非常常用 gene_sets KEGG_2021_Human # 3. 执行ORA分析 ora_results gp.enrichr(gene_listsig_up_genes, gene_setsgene_sets, organismHuman, # 必须指定 outdirNone, # 不输出文件结果保存在变量中 cutoff0.05 # 富集结果的FDR/q-value阈值 ) # enrichr返回一个包含多个DataFrame的字典我们通常需要的是.results ora_df ora_results.results # 查看最显著的前10个富集通路 print(ora_df.sort_values(Adjusted P-value).head(10)[[Term, Overlap, P-value, Adjusted P-value, Odds Ratio]])关键参数与结果解读organism必须正确指定否则会找不到对应的基因集。cutoff输出结果的显著性阈值基于校正后的P值Adjusted P-value通常是FDR。结果列解读Term富集到的通路/基因集名称。Overlap格式为X/Y表示你的基因列表中有X个基因属于该通路该通路总共有Y个基因。P-value富集分析的原始p值。Adjusted P-value经过多重检验校正后的p值如FDR这是我们判断结果是否显著的主要依据通常要求0.05。Odds Ratio比值比大于1表示正富集你的基因在该通路中富集数值越大富集程度越强。Genes你的基因列表中属于该通路的基因具体是哪些。实操心得enrichr函数默认会同时查询多个数据库如果gene_sets是一个列表。但初次使用时我建议先专注于一个数据库如‘KEGG_2021_Human’避免结果过于庞杂。另外outdir参数可以指定路径这样函数会自动将结果保存为文本文件和图表非常方便。3.3 执行GSEA分析排序列表的构建与富集信号挖掘GSEA需要的是一个排好序的基因列表。通常我们使用log2FoldChange作为排序指标因为它包含了变化的方向和幅度。# 1. 为GSEA准备数据一个包含Gene Symbol和排序指标的DataFrame # 我们使用清洗后的全部基因按log2FoldChange从大到小排序正值越大排越前负值越小排越后 gsea_data de_df_clean[[gene_symbol, logfoldchanges]].copy() gsea_data gsea_data.set_index(gene_symbol) # GSEA需要的是一个Series索引是基因名值是排序指标 gsea_rank gsea_data[logfoldchanges].sort_values(ascendingFalse) # 2. 执行GSEA分析 # 这里使用MSigDB的Hallmark基因集它包含了50个精炼的、具有明确生物学状态的基因集 gsea_results gp.gsea(datagsea_rank, gene_setsMSigDB_Hallmark_2020, organismhuman, permutation_num1000, # 置换检验次数默认1000增加次数更稳定但更慢 outdir./gsea_output, # 指定输出目录会生成详细结果和图表 formatpng, seed42 # 设置随机种子保证结果可复现 ) # gsea函数返回一个GSEA结果对象其.res2d属性是主要的结果DataFrame gsea_df gsea_results.res2d # 查看富集最显著NES绝对值大FDR小的前10个基因集 # NES (Normalized Enrichment Score) 标准化富集分数是核心指标 significant_results gsea_df[(gsea_df[FDR q-val] 0.25) (abs(gsea_df[NES]) 1.0)] print(significant_results.sort_values(NES, ascendingFalse).head(10)[[Term, NES, FDR q-val, Lead_genes]])GSEA核心结果解读NES标准化富集分数这是GSEA的灵魂指标。正负号表示方向NES 0表示该基因集在排序列表的顶部富集即在你的差异分析中该基因集对应的基因普遍上调。NES 0表示在底部富集普遍下调。绝对值表示强度|NES|越大富集程度越强。通常认为|NES| 1.0有生物学意义。FDR q-val错误发现率。GSEA官方推荐使用FDR 0.25作为显著性阈值这比ORA的0.05宽松因为GSEA探测的是更微妙的协调性变化。当然在严格的研究中也可以使用0.05。Lead_genes核心贡献基因。即在该基因集的富集信号中贡献最大的那些基因位于排序列表最顶端或最底端且属于该基因集。它们是解释该富集结果的关键。ES富集分数未标准化的原始分数受基因集大小影响一般看NES即可。参数设置经验permutation_num置换检验次数。默认1000次对于大多数分析足够。如果你追求极致的p值精度或者基因集很小/很大可以增加到5000或10000次但计算时间会线性增加。seed务必设置。这能确保每次运行的结果完全一致对于可复现性至关重要。outdirgsea函数的一个巨大优点是会自动生成可视化结果图富集分析图Enrichment Plot这对于理解和展示结果非常有帮助。4. 高级结果可视化与生物学解读让数据自己说话得到富集分析表格只是第一步如何将结果高效、美观地呈现出来并挖掘其背后的生物学故事才是更见功力的地方。gseapy和matplotlib/seaborn的结合可以创造出强大的可视化效果。4.1 绘制经典富集分析气泡图/条形图气泡图是展示富集结果最直观的方式之一它同时展示了通路的显著性-log10(p-value)或FDR、富集强度Odds Ratio或NES和基因重叠数量。import matplotlib.pyplot as plt import seaborn as sns # 假设 ora_df 是我们之前ORA分析的结果DataFrame # 选取FDR最显著的前15个通路进行可视化 top_ora ora_df.sort_values(Adjusted P-value).head(15).copy() # 计算 -log10(Adjusted P-value) 用于作图 top_ora[-log10(FDR)] -np.log10(top_ora[Adjusted P-value]) # 从‘Overlap’列中提取重叠基因数 top_ora[Overlap_Count] top_ora[Overlap].apply(lambda x: int(x.split(/)[0])) plt.figure(figsize(10, 8)) # 绘制气泡图x轴为Odds Ratioy轴为通路名称气泡大小为重叠基因数颜色为 -log10(FDR) scatter plt.scatter(xtop_ora[Odds Ratio], ytop_ora[Term], stop_ora[Overlap_Count]*20, # 气泡大小乘以一个系数放大视觉效果 ctop_ora[-log10(FDR)], cmapviridis_r, # 使用颜色映射_r表示反转使得值越大颜色越深 alpha0.7, edgecolorsk, linewidth0.5) plt.axvline(x1, colorgrey, linestyle--, linewidth0.8) # 添加Odds Ratio1的参考线 plt.xlabel(Odds Ratio, fontsize12) plt.ylabel(Pathway/Term, fontsize12) plt.title(ORA Enrichment Analysis (Top 15 Pathways), fontsize14, pad20) # 添加颜色条 cbar plt.colorbar(scatter) cbar.set_label(-log10(FDR), fontsize11) # 添加图例表示气泡大小需要手动创建 import matplotlib.patches as mpatches # 创建几个示例大小的图例句柄 size_legend_handles [mpatches.Circle((0,0), radiuss, facecolorgray, edgecolork, alpha0.7) for s in [5, 10, 15]] # 半径对应 Overlap_Count * 20 后的视觉大小 size_labels [5, 10, 15] # 对应的重叠基因数 plt.legend(handlessize_legend_handles, labelssize_labels, titleGene Count, loccenter left, bbox_to_anchor(1.05, 0.5)) plt.tight_layout() plt.show()这张图可以一目了然地告诉我们哪些通路最显著颜色深、富集程度最强Odds Ratio远大于1、并且有足够多的基因支持气泡大。4.2 解读GSEA富集分析图gseapy在运行GSEA时生成的PNG图片位于指定的outdir中是标准输出。以“HALLMARK_INFLAMMATORY_RESPONSE”为例生成的图片包含三部分顶部主图Enrichment Score Plot中间的曲线是富集分数ES在排序基因列表上的运行轨迹。峰值出现在哪里就表示该基因集的成员基因在排序列表的那个位置集中出现。峰值在左侧排名靠前表示上调富集NES0在右侧表示下调富集NES0。图中的绿色竖线就是Lead_genes核心基因的位置。中间部分Gene List Metric显示了所有基因按排序指标这里是logFC的分布情况像一条山脉。这直观展示了你输入数据的排序质量。底部热图Heatmap显示了该基因集中每个基因在原始数据集如果你提供了表达矩阵或在你排序列表中的相对位置/表达模式。在单细胞语境下如果我们提供了所有细胞的表达矩阵这里可以显示出该基因集在哪些细胞簇中高表达。解读技巧不要只看FDR和NES的数字。一定要打开这些图观察ES曲线的形状。一个“漂亮”的富集信号通常表现为一个早期达到峰值或谷值的、单峰的、清晰的曲线。如果曲线来回震荡峰值不明显即使FDR勉强显著其生物学意义也可能值得怀疑。4.3 整合多组结果与自定义基因集分析有时我们需要比较不同细胞簇或不同实验条件之间的富集结果。这时可以绘制多组对比气泡图或热图。# 假设我们对三个细胞簇cluster0,1,2都做了ORA分析结果存储在字典中results_dict {cluster0: ora_df0, cluster1: ora_df1, cluster2: ora_df2} # 我们想找出那些在多个簇中都富集的通路 common_terms set() for cluster, df in results_dict.items(): # 获取每个簇中FDR0.05的显著通路 sig_terms set(df[df[Adjusted P-value] 0.05][Term]) if not common_terms: common_terms sig_terms else: common_terms common_terms.intersection(sig_terms) print(fPathways enriched in all three clusters: {common_terms}) # 为了可视化我们可以创建一个矩阵行是通路列是细胞簇值是 -log10(FDR) vis_data [] for term in list(common_terms)[:20]: # 取前20个共同通路示例 row {Term: term} for cluster, df in results_dict.items(): # 在对应簇的结果中查找该通路的FDR值 term_row df[df[Term] term] if not term_row.empty: row[cluster] -np.log10(term_row.iloc[0][Adjusted P-value]) else: row[cluster] 0 # 如果不显著赋值为0 vis_data.append(row) vis_df pd.DataFrame(vis_data).set_index(Term) # 绘制热图 plt.figure(figsize(8, 10)) sns.heatmap(vis_df, cmapYlOrRd, linewidths0.5, annotTrue, fmt.2f, cbar_kws{label: -log10(FDR)}) plt.title(Common Enriched Pathways Across Clusters, fontsize14, pad20) plt.tight_layout() plt.show()此外gseapy支持使用自定义基因集。你可以将自己关注的、来自文献或特定研究的基因列表保存为.gmt格式文件然后通过gene_sets‘/path/to/your_custom.gmt’参数调用。这在进行非常聚焦的、与特定生物学问题相关的研究时极其有用。5. 实战避坑与进阶技巧来自一线踩坑的经验纸上得来终觉浅绝知此事要躬行。下面分享几个我在实际项目中反复遇到并且耗费不少时间才解决的问题。5.1 物种匹配错误与基因同源映射问题分析小鼠单细胞数据时直接使用Gene Symbol进行富集分析结果一片空白或非常奇怪。根因gseapy的许多内置数据库如KEGG, GO默认是面向人类的。虽然它也支持小鼠‘Mouse’、大鼠等但基因标识符必须与数据库匹配。小鼠的基因符号如Trp53和人类的TP53不同。解决方案使用正确的物种参数在执行gp.enrichr或gp.gsea时务必设置organism‘Mouse’。gseapy会自动调用对应物种的数据库。进行严谨的基因同源映射如果你的物种不在gseapy直接支持的范围如斑马鱼、果蝇或者你想跨物种比较就需要进行同源基因映射。可以使用mygene的homologene功能或专业的同源数据库如Ensembl BioMart将你的基因列表映射到人类同源基因上再使用人类数据库进行分析。这是一个复杂但有时不可避免的步骤。# 使用mygene进行小鼠到人类的同源映射示例简化 mg mygene.MyGeneInfo() mouse_genes [Trp53, Myc, Sox2] # 查询小鼠基因的人类同源物 homologs mg.querymany(mouse_genes, scopessymbol, fieldshomologene.human_symbol, speciesmouse) human_gene_list [] for h in homologs: if homologene in h and human_symbol in h[homologene]: human_gene_list.append(h[homologene][human_symbol]) else: print(fNo human homolog found for {h.get(query)}) print(fMapped human genes: {human_gene_list})5.2 背景基因集的正确设置问题ORA分析结果中某些通路的Odds Ratio异常高比如几百或者p值极其显著但重叠基因数很少比如1/1结果不可信。根因默认情况下enrichr使用其服务器上对应物种的“所有基因”作为背景background。但在单细胞RNA-seq中由于技术限制我们只能检测到约10000-20000个基因远少于全基因组背景~20000个蛋白编码基因。使用不匹配的背景集会导致统计偏差。解决方案自定义背景基因集。你应该将本次单细胞实验中所有被检测到的基因即adata.var_names映射为Gene Symbol后的列表作为背景。# 获取本次分析中所有检测到的基因背景集 # 假设我们已经将adata的索引转换为Gene Symbol或有一个映射关系 # 这里假设 adata.var[gene_symbol] 存储了Gene Symbol background_genes adata.var[gene_symbol].dropna().unique().tolist() print(fBackground gene set size: {len(background_genes)}) # 在执行enrichr时通过 background 参数传入 ora_results_custom_bg gp.enrichr(gene_listsig_up_genes, gene_setsKEGG_2021_Human, organismHuman, backgroundbackground_genes, # 关键参数 outdirNone, cutoff0.05)使用自定义背景后Odds Ratio和p值的计算会更加准确能有效过滤掉因背景集过大而造成的假阳性富集信号。5.3 处理大规模分析与结果过滤问题一次分析几十个细胞簇每个簇都做一次富集分析产生海量结果难以管理和聚焦。解决方案自动化脚本 结构化结果过滤。自动化循环将分析过程封装成函数对每个簇进行循环分析。结构化存储将每个簇的显著结果如FDR0.05存储到一个统一的、结构化的DataFrame或字典中包含簇名、通路名、NES、FDR、重叠基因等关键信息。结果过滤与聚合基于FDR和NES阈值这是最基本的一步。去除过于宽泛的通路GO或KEGG中有些通路非常宽泛如“代谢过程”、“信号转导”虽然显著但信息量低。可以手动建立一个“宽泛通路黑名单”进行过滤或者通过通路的基因数量Term大小进行筛选例如只保留基因数在10到500之间的通路。聚类相似通路使用gseapy的plot模块中的dotplot函数或GOplot、clusterProfilerR等工具可以对富集到的通路进行语义相似性聚类将相似的通路归为一类使结果更简洁。# 示例循环分析多个簇并收集显著结果 all_significant_results [] for cluster in adata.obs[leiden].cat.categories: # 1. 提取该簇的差异基因 (假设已有差异分析结果) de_genes get_de_genes_for_cluster(adata, cluster) # 这是一个自定义函数示例 if len(de_genes) 5: # 如果差异基因太少跳过 continue # 2. 执行富集分析 ora_res gp.enrichr(gene_listde_genes, gene_setsGO_Biological_Process_2021, ...) ora_df ora_res.results # 3. 过滤显著结果 sig_df ora_df[ora_df[Adjusted P-value] 0.05].copy() sig_df[cluster] cluster # 标记来源簇 # 4. 收集 all_significant_results.append(sig_df[[cluster, Term, Adjusted P-value, Odds Ratio, Overlap]]) # 合并所有结果 final_results_df pd.concat(all_significant_results, ignore_indexTrue) # 后续可以方便地用pandas进行筛选、排序和可视化5.4 与Scanpy生态的深度整合将通路分数嵌入AnnData富集分析的结果不仅可以用来画图和写文章还可以反向嵌入到Scanpy的AnnData对象中为后续的细胞层面分析提供新维度。例如我们可以计算每个细胞在某个特定通路如“上皮-间质转化EMT”上的“活性分数”。一种常见的方法是使用单样本基因集富集分析ssGSEA或AUCell算法。虽然gseapy本身不直接提供与Scanpy的深度接口但我们可以利用其基因集和Scanpy的sc.tl.score_genes功能进行近似计算。# 假设我们从GSEA结果中找到了一个感兴趣的基因集例如‘HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION’ # 我们需要这个基因集的成员基因列表 emt_gene_set ... # 从gsea_results.results_on_cell或MSigDB官网获取该基因集的基因列表 # 使用Scanpy的score_genes函数计算每个细胞在这个基因集上的“富集分数” sc.tl.score_genes(adata, emt_gene_set, score_nameEMT_score, ctrl_size50, use_rawFalse) # 参数说明 # - adata: AnnData对象 # - emt_gene_set: 基因列表 # - score_name: 存储在adata.obs中的列名 # - ctrl_size: 用于标准化分数的随机背景基因集大小 # - use_raw: 是否使用adata.raw # 现在adata.obs中多了一列‘EMT_score’ # 我们可以可视化这个分数在UMAP图上的分布 sc.pl.umap(adata, color[EMT_score, leiden], cmapRdBu_r, vcenter0) # 或者看这个分数在不同细胞簇间的差异 sc.pl.violin(adata, [EMT_score], groupbyleiden)这样我们就将一个通路水平的富集信息转化为了每个细胞的一个连续型特征EMT_score可以用于后续的聚类验证、差异分析、甚至作为回归模型的输入极大地拓展了分析的深度和灵活性。这个技巧在分析肿瘤异质性、细胞分化轨迹等复杂生物学过程中尤其有用。
