GPU加速的11种单细胞通路富集活性打分

摘要LLM 归库摘要

核心内容： 这是一篇 OmicVerse ov.es 模块教程，用同一份小鼠胰腺发育单细胞数据演示 11 种通路/基因集活性打分方法，并比较 CPU 与 GPU 后端结果。

适用场景： 单细胞或空间转录组中，需要把 GO/MSigDB/regulon 等基因集转换为每个细胞的通路活性分数，再用于 UMAP、分组差异或疾病-正常对比。

输入： AnnData 对象、基因集字典 signatures、方法名、engine='cpu'/'gpu'、最小基因数阈值 tmin 等参数。

输出： 各方法分数写入 adata.obsm['score_方法名']，部分方法还会写入 adata.obsm['padj_方法名']，并可同步写入 adata.obs 方便可视化。

关键参数/步骤： 推荐先用 AUCell、GSVA、Z-score 快速初筛；带权重或 regulon 分析可用 ULM/MLM/VIPER；GPU 后端适合大规模矩阵化计算。

常见错误： GSEA 的 CPU/GPU 结果可能因 NES 与 raw ES 尺度不同而数值不完全可比；ORA 需要明确显著基因集合，否则单细胞分辨率有限；MDT/UDT 的 CPU 逐细胞拟合成本较高。

与口腔生信连接： 可用于牙周炎、OSCC、口腔免疫微环境或空间转录组项目中，对炎症、免疫浸润、EMT、增殖、缺氧等通路进行细胞级活性评分。

关联页面建议： 后续可整合到 03 生信教程与工具库/02_wiki/workflows/单细胞通路活性打分或 03 生信教程与工具库/02_wiki/tools/OmicVerse。

下一步处理： 复现时优先确认 OmicVerse 版本、GPU 环境、基因集文件路径，并把常用方法封装为项目模板函数。

归库说明： 分类到 03 生信教程与工具库，置信度 0.92；已标记为暂未试用，便于后续集中复现。

编者按

我们将推出一个新的栏目，重构算法的介绍，在这个栏目中，我们将介绍omicverse开源重构生态最新的成果以及完整的教程，欢迎关注。

第二期教程里，我们给大家带来的是单细胞通路活性打分算法的实战比较：这是omicverse开源重构生态中 ov.es 模块的内容。

在单细胞分析里，我们经常不会只关心某一个基因的表达，而是更想知道一组基因，也就是某条通路、某个调控子、某个功能模块在每个细胞里的活性。过去这些方法分散在不同包里，参数也不完全一样，所以我一直觉得这块需要一个更统一的入口。

OmicVerse 在 ov.es 中重构了 decoupler 里的 11 种 enrichment score 方法，并且提供了 CPU/GPU 两套计算后端。也就是说，如果你想把 AUCell 换成 GSVA，或者把 CPU 换成 GPU，我们只需要改一个参数。

本文AI率20%，阅读大约需要5min。

OmicVerse团队

1. 引言#

做单细胞分析时，我们最常见的操作是找 marker 基因。

但 marker 基因有一个问题：它很零散。

比如你看胰腺发育数据，某个细胞群里可能有几十个基因一起变化。你当然可以一个个基因去看，但如果想回答“这个细胞是不是更偏向某个生物学过程”，那么用单个基因就不太够了。

这时候通路活性打分就很有用。

我们先准备一个基因集，比如 response to vitamin (GO:0033273) ，我们知道这个基因集都有哪些基因，然后计算这个基因集在每个细胞里的活性分数。后面就可以像画基因表达一样，把这个分数画到 UMAP 上。

这类分析在很多场景里都能用：

细胞类型注释后，看某类细胞是否激活特定通路。
做疾病和正常对比时，看通路活性是否发生改变。
做转录因子调控分析时，把 regulon 当成基因集去打分。
做空间转录组时，把通路分数投到空间坐标上。

过去大家常用的方法包括 AUCell、GSEA、GSVA、ORA、VIPER 等。它们背后的统计思想不一样，适合的场景也不一样。OmicVerse 现在把这些方法统一到了 ov.es 下面，并且提供了同一个调用方式。

ov.es.decoupler(
    adata,
    signatures=signatures,
    method='aucell',
    engine='gpu',
    tmin=3,
)

如果你更习惯直接调用具体函数，也可以这样写：

ov.es.aucell(adata, signatures=signatures, engine='gpu', tmin=3)

二者是等价的。

2. 11种方法怎么选#

在decoupler中，作者的benchmark认为ulm是最优的算法，我分析发现ulm确实有可取之处，比AUCell有时候还要好。

除此之外，如果你只是想快速看某个基因集在单细胞里的分布，我会优先从 AUCell、GSVA、Z-score 这几个方法开始。它们比较直观，也比较适合画图。

如果你手里的基因集带有权重，或者想做调控子活性分析，可以考虑 ULM、MLM 和 VIPER。

如果你只是想做一个最简单的加权平均，WAGGR 就够了。

如果你需要更复杂的非线性打分，可以看 MDT 和 UDT。不过这两个方法如果用传统 CPU 逐细胞拟合树模型，会非常慢，所以这个教程里只展示 GPU 版本。

方法	直观理解	适合场景
`AUCell`	看基因集成员在每个细胞的表达排序里是否靠前	单细胞通路活性、regulon 活性
`GSEA`	沿着排序基因列表做累积得分	经典富集分析迁移到单细胞
`GSVA`	GSEA 的单样本版本	样本/细胞级通路打分
`ORA`	显著基因与通路基因的 2×2 表检验	有明确显著基因集合时
`ULM`	一个基因集一个线性模型	带权重的 signature
`MLM`	多个基因集一起进入线性模型	基因集之间有重叠时
`WAGGR`	加权平均或加权求和	想要最直接的分数
`Z-score`	基因集平均表达相对背景的偏移	快速、稳定、容易解释
`VIPER`	把 signature 当成 regulon 计算活性	转录因子/调控网络
`MDT`	多变量决策树打分	GPU 下的非线性 signature 评分
`UDT`	单变量决策树打分	GPU 下的非线性 signature 评分

我个人的使用习惯是，先用一两个简单方法扫一遍，比如 AUCell 和 Z-score；如果结果很稳定，再根据具体问题补充更复杂的方法。

3. 实战开始#

这里使用的是小鼠胰腺内分泌发生数据，来自 Bastidas-Ponce et al.，OmicVerse 已经提供了下载接口。这个数据自带 UMAP，所以我们不用重新跑 PCA、邻居图和 UMAP。

import numpy as np
import omicverse as ov

ov.style()

读取数据：

adata = ov.datasets.pancreatic_endocrinogenesis()
adata

AnnData object with n_obs × n_vars = 3696 × 27998
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
    var: 'highly_variable_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'distances', 'connectivities'

为了让 CPU 版本的对比跑得更快，这里随机抽取 300 个细胞：

rng = np.random.default_rng(0)
adata = adata[rng.choice(adata.n_obs, size=300, replace=False)].copy()

adata.var['robust'] = True
adata = ov.pp.preprocess(adata, target_sum=1e4)

预处理以后保留了 300 个细胞和 14,182 个基因。

然后读取 GO Biological Process 2021 的基因集：

pathway_dict = ov.utils.geneset_prepare(
    'genesets/GO_Biological_Process_2021.txt',
    organism='Mouse',
)

这里我们只选一条通路做演示：

geneset_name = 'response to vitamin (GO:0033273)'
signatures = {geneset_name: pathway_dict[geneset_name]}

注意，这里为了比较 11 种方法，我们只用了一个基因集。实际项目里你当然可以同时传入很多通路，OmicVerse 会给每条通路都写入一个分数。

4. 统一接口#

ov.es 的一个好处是接口比较统一。

ov.es.decoupler(
    adata,
    signatures=signatures,
    method='aucell',
    engine='gpu',
    tmin=3,
)

计算完成后，分数会写入：

adata.obsm['score_aucell']

如果这个方法会计算校正后的 P 值，也会写入类似的位置：

adata.obsm['padj_aucell']

所以我们后面切换方法时，只需要把 method='aucell' 改成 method='gsva' 、 method='zscore' 等即可。

5. AUCell#

AUCell 的思路很直接：对每个细胞里的基因表达做排序，然后看基因集里的基因是不是更容易出现在排序靠前的位置。

ov.es.aucell(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['aucell'] = adata.obsm['score_aucell'][geneset_name].astype(float).values
adata.obs['aucell (CPU)'] = scores_cpu['aucell']

ov.es.aucell(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['aucell'] = adata.obsm['score_aucell'][geneset_name].astype(float).values
adata.obs['aucell (GPU)'] = scores_gpu['aucell']

可视化：

ov.pl.umap(
    adata,
    color=['aucell (CPU)', 'aucell (GPU)'],
    cmap='Reds',
    frameon=False,
)

AUCell CPU 与 GPU 对比

可以看到 CPU 和 GPU 的空间分布基本一致。对于这种基于排序的打分，GPU 的意义主要在于把大规模数据上的等待时间压下来。

6. GSEA、GSVA 和 ORA#

GSEA 和 GSVA 是大家很熟悉的富集分析方法。放到单细胞数据里以后，可以理解成对每个细胞都做一次基因集打分。

先跑 GSEA：

ov.es.gsea(adata, signatures=signatures, tmin=3, engine='cpu', times=50)
scores_cpu['gsea'] = adata.obsm['score_gsea'][geneset_name].astype(float).values
adata.obs['gsea (CPU)'] = scores_cpu['gsea']

ov.es.gsea(adata, signatures=signatures, tmin=3, engine='gpu', times=1)
scores_gpu['gsea'] = adata.obsm['score_gsea'][geneset_name].astype(float).values
adata.obs['gsea (GPU)'] = scores_gpu['gsea']

GSEA CPU 与 GPU 对比

这里要注意一个细节：CPU 版本在 times>1 时返回的是标准化后的 NES，而 GPU 版本在 times=1 时返回的是原始 ES。所以它们的空间模式相近，但数值尺度不完全一样，后面的散点图里也能看到这一点。

GSVA：

ov.es.gsva(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['gsva'] = adata.obsm['score_gsva'][geneset_name].astype(float).values
adata.obs['gsva (CPU)'] = scores_cpu['gsva']

ov.es.gsva(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['gsva'] = adata.obsm['score_gsva'][geneset_name].astype(float).values
adata.obs['gsva (GPU)'] = scores_gpu['gsva']

GSVA CPU 与 GPU 对比

GSVA 在这个数据上的 CPU 用时是 5.4377 秒，GPU 用时是 0.1089 秒。这个差距就很明显了。

ORA：

ov.es.ora(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['ora'] = adata.obsm['score_ora'][geneset_name].astype(float).values
adata.obs['ora (CPU)'] = scores_cpu['ora']

ov.es.ora(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['ora'] = adata.obsm['score_ora'][geneset_name].astype(float).values
adata.obs['ora (GPU)'] = scores_gpu['ora']

ORA CPU 与 GPU 对比

ORA 更依赖“哪些基因被认为显著”这个前提。如果你的输入本身没有一个很明确的显著基因集合，那么它的单细胞分辨率通常不如 AUCell 或 GSVA 那么细。

7. ULM、MLM、WAGGR 和 Z-score#

接下来是几种更适合带权重 signature 的方法。

ULM 是一元线性模型，也就是一个 signature 单独拟合一次：

ov.es.ulm(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['ulm'] = adata.obsm['score_ulm'][geneset_name].astype(float).values
adata.obs['ulm (CPU)'] = scores_cpu['ulm']

ov.es.ulm(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['ulm'] = adata.obsm['score_ulm'][geneset_name].astype(float).values
adata.obs['ulm (GPU)'] = scores_gpu['ulm']

ULM CPU 与 GPU 对比

MLM 是多元线性模型，也就是把多个 signature 一起放进模型里。如果你的基因集之间重叠很多，MLM 会更适合一点。

ov.es.mlm(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['mlm'] = adata.obsm['score_mlm'][geneset_name].astype(float).values
adata.obs['mlm (CPU)'] = scores_cpu['mlm']

ov.es.mlm(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['mlm'] = adata.obsm['score_mlm'][geneset_name].astype(float).values
adata.obs['mlm (GPU)'] = scores_gpu['mlm']

MLM CPU 与 GPU 对比

WAGGR 就更直接了，其实就是基因表达的加权平均或加权求和。

ov.es.waggr(adata, signatures=signatures, tmin=3, engine='cpu', times=1)
scores_cpu['waggr'] = adata.obsm['score_waggr'][geneset_name].astype(float).values
adata.obs['waggr (CPU)'] = scores_cpu['waggr']

ov.es.waggr(adata, signatures=signatures, tmin=3, engine='gpu', times=1)
scores_gpu['waggr'] = adata.obsm['score_waggr'][geneset_name].astype(float).values
adata.obs['waggr (GPU)'] = scores_gpu['waggr']

WAGGR CPU 与 GPU 对比

Z-score 是我觉得很适合初筛的方法。它计算的是基因集平均表达相对于每个细胞背景分布的偏移。

ov.es.zscore(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['zscore'] = adata.obsm['score_zscore'][geneset_name].astype(float).values
adata.obs['zscore (CPU)'] = scores_cpu['zscore']

ov.es.zscore(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['zscore'] = adata.obsm['score_zscore'][geneset_name].astype(float).values
adata.obs['zscore (GPU)'] = scores_gpu['zscore']

Z-score CPU 与 GPU 对比

这几个方法的 CPU/GPU 结果都几乎完全一致。

8. VIPER、MDT 和 UDT#

VIPER 更常用于转录因子活性分析。它会把 signature 当成 regulon，然后计算调控子在每个细胞里的活性。

ov.es.viper(adata, signatures=signatures, tmin=3, engine='cpu')
scores_cpu['viper'] = adata.obsm['score_viper'][geneset_name].astype(float).values
adata.obs['viper (CPU)'] = scores_cpu['viper']

ov.es.viper(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['viper'] = adata.obsm['score_viper'][geneset_name].astype(float).values
adata.obs['viper (GPU)'] = scores_gpu['viper']

VIPER CPU 与 GPU 对比

MDT 和 UDT 是基于决策树的打分方法。它们在 CPU 上逐细胞拟合树模型的成本比较高，所以这个教程只跑 GPU。

MDT：

ov.es.mdt(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['mdt'] = adata.obsm['score_mdt'][geneset_name].astype(float).values
adata.obs['mdt (GPU)'] = scores_gpu['mdt']

MDT GPU 结果

UDT：

ov.es.udt(adata, signatures=signatures, tmin=3, engine='gpu')
scores_gpu['udt'] = adata.obsm['score_udt'][geneset_name].astype(float).values
adata.obs['udt (GPU)'] = scores_gpu['udt']

UDT GPU 结果

我觉得这两个方法更适合在你已经有明确需求时使用。对于大多数通路活性可视化，AUCell、GSVA、Z-score 已经足够了。

9. CPU 和 GPU 的一致性#

这里把所有有 CPU 基线的方法放到一张散点图里。

每个点是一个细胞，横坐标是 CPU 分数，纵坐标是 GPU 分数。为了让不同方法能放到同一张图中，这里对每个方法的分数做了标准化。

import matplotlib.pyplot as plt

_methods = [m for m in scores_gpu if m in scores_cpu]
_cmap = plt.get_cmap('tab20')

fig, ax = plt.subplots(figsize=(5.5, 5.5))
for _i, _m in enumerate(_methods):
    _x = scores_cpu[_m]
    _y = scores_gpu[_m]
    _xz = (_x - _x.mean()) / (_x.std() + 1e-12)
    _yz = (_y - _y.mean()) / (_y.std() + 1e-12)
    ax.scatter(_xz, _yz, s=4, alpha=0.4, label=_m, color=_cmap(_i % _cmap.N))

ax.plot([-6, 6], [-6, 6], '--', color='grey', lw=0.8)

CPU 与 GPU 分数一致性

从图上看，绝大部分方法都落在对角线附近。GSEA 的点会有一点拉伸，原因前面也说过：CPU 版本这里返回的是 NES，GPU 版本返回的是 raw ES。

再看数值表：

方法	Pearson r	Spearman r	RMSE / range
AUCell	1.0000	1.0000	0.0000
GSEA	0.9747	0.8019	0.2224
GSVA	1.0000	1.0000	0.0001
ORA	NaN	NaN	0.0000
ULM	1.0000	1.0000	0.0000
MLM	1.0000	1.0000	0.0000
WAGGR	1.0000	1.0000	0.0000
Z-score	1.0000	1.0000	0.0000
VIPER	1.0000	1.0000	0.0000

这里 ORA 的相关性是 NaN。原因很简单：这个例子里 ORA 的分数没有足够的变化，所以相关系数无法计算。

如果只看这个 300 细胞的小例子，各方法的运行时间如下：

方法	CPU	GPU
AUCell	0.1586s	0.0264s
GSEA	0.8218s	0.1066s
GSVA	5.4377s	0.1089s
ORA	0.6613s	0.0827s
ULM	0.1155s	0.0459s
MLM	0.1668s	0.0716s
WAGGR	1.1193s	0.0299s
Z-score	0.0749s	0.0368s
VIPER	0.3006s	0.0580s
MDT	-	1.6364s
UDT	-	0.1838s

这个表只对应当前教程、当前机器、当前数据规模，不要把它当成正式基准测试。真正的大规模数据里，GPU 的优势会更明显，尤其是 GSVA、WAGGR 这类可以大规模矩阵化的计算。

10. 总结#

在本期教程中，我们用同一份小鼠胰腺发育单细胞数据，对比了 OmicVerse ov.es 里的 11 种 enrichment score 方法：

AUCell、GSEA、GSVA、ORA、ULM、MLM、WAGGR、Z-score、VIPER 都提供 CPU/GPU 两套后端。
MDT 和 UDT 在这个教程里只展示 GPU，因为 CPU 逐细胞拟合树模型的成本不太适合单细胞规模。
除 GSEA 的 NES/raw ES 尺度差异外，绝大多数方法的 CPU/GPU 结果几乎完全一致。
如果只是做通路活性可视化，我会优先用 AUCell、GSVA 或 Z-score；如果做调控子活性，可以考虑 VIPER、ULM 或 MLM。

实际上我觉得 ov.es 这个模块最大的价值在于接口统一。过去分散的 enrichment score 方法，现在可以在同一个对象、同一套参数里完成切换。

你不需要每次换算法就重新学一套输入输出，也不需要在不同包之间来回转换对象。对于一个真实项目来说，这会省掉很多非常无聊但又很容易出错的时间。

作为 omicverse 开源重构生态的一部分，我们会继续把更多常用算法迁移到统一接口里。如果你有想重构的 R 包，或者希望某个算法接入 omicverse，也欢迎联系我们。

11. 交流群#

如果你也很好奇我们omicverse生态的最新进展，想第一时间体验到新功能，欢迎加入我们的交流群～

OmicVerse相关教程 · 目录

阅读原文

继续滑动看下一个

OmicAI

向上滑动看下一个

闪记

复制 LaTeX 公式

12. 完整代码整理｜下次复用版#

说明：下面代码根据原文流程整理为一个相对完整的复现模板。第一次运行前需要确认 omicverse、numpy、pandas、matplotlib 等依赖已安装，并且 genesets/GO_Biological_Process_2021.txt 路径存在。若没有 GPU，可把 ENGINE_GPU = "gpu" 改为 "cpu" 或只运行 CPU 部分。

# ============================================================
# 0. 环境与参数设置
# ============================================================
# 作用：集中管理依赖、随机种子、基因集名称、CPU/GPU 后端等复现参数。

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import omicverse as ov

# 设置 OmicVerse 绘图风格
ov.style()

# 固定随机种子，保证抽样结果可复现
RANDOM_SEED = 0
N_CELLS = 300

# 基因集文件路径：需要按本地实际路径修改
GENESET_PATH = "genesets/GO_Biological_Process_2021.txt"
ORGANISM = "Mouse"
GENESET_NAME = "response to vitamin (GO:0033273)"

# 统一控制 CPU/GPU 后端
ENGINE_CPU = "cpu"
ENGINE_GPU = "gpu"

# 基因集最少命中基因数阈值；太小容易不稳定，太大可能过滤掉目标通路
TMIN = 3


# ============================================================
# 1. 读取并预处理示例单细胞数据
# ============================================================
# 作用：下载 OmicVerse 内置小鼠胰腺内分泌发生数据，随机抽取细胞并做基础预处理。

adata = ov.datasets.pancreatic_endocrinogenesis()

# 为了快速演示，只随机抽取 300 个细胞；正式项目可删除这一步，直接使用全量数据
rng = np.random.default_rng(RANDOM_SEED)
adata = adata[rng.choice(adata.n_obs, size=N_CELLS, replace=False)].copy()

# 标记全部基因为 robust，随后进行标准化等预处理
adata.var["robust"] = True
adata = ov.pp.preprocess(adata, target_sum=1e4)

print(adata)


# ============================================================
# 2. 准备通路/基因集 signatures
# ============================================================
# 作用：读取 GO Biological Process 2021 基因集，并抽取一条通路作为演示。
# 正式分析时，可以把 signatures 设置为多个通路组成的字典。

pathway_dict = ov.utils.geneset_prepare(
    GENESET_PATH,
    organism=ORGANISM,
)

# 单通路演示：字典格式为 {通路名: 基因列表或带权重 signature}
signatures = {GENESET_NAME: pathway_dict[GENESET_NAME]}

print(f"Selected geneset: {GENESET_NAME}")
print(f"Number of genes in geneset: {len(signatures[GENESET_NAME])}")


# ============================================================
# 3. 定义通用运行函数
# ============================================================
# 作用：减少重复代码；自动运行指定 enrichment score 方法，提取分数、记录耗时。

def run_es_method(adata, method, engine, signatures, geneset_name, tmin=3, **kwargs):
    """
    运行 OmicVerse ov.es 中的一种通路活性打分方法。
    返回当前 geneset 在每个细胞中的活性分数，以及运行耗时。
    """
    func = getattr(ov.es, method)

    start = time.perf_counter()
    func(
        adata,
        signatures=signatures,
        tmin=tmin,
        engine=engine,
        **kwargs,
    )
    elapsed = time.perf_counter() - start

    score_key = f"score_{method}"
    score = adata.obsm[score_key][geneset_name].astype(float).values

    return score, elapsed


# ============================================================
# 4. 批量运行 11 种 enrichment score 方法
# ============================================================
# 作用：统一比较 CPU/GPU 的结果与耗时。
# 注意：MDT 和 UDT 在原文中只运行 GPU，因为 CPU 逐细胞拟合树模型成本较高。

scores_cpu = {}
scores_gpu = {}
runtime_records = []

# 每个方法的额外参数；没有额外参数的留空字典
method_kwargs_cpu = {
    "aucell": {},
    "gsea": {"times": 50},      # CPU times>1 时常返回 NES
    "gsva": {},
    "ora": {},
    "ulm": {},
    "mlm": {},
    "waggr": {"times": 1},
    "zscore": {},
    "viper": {},
}

method_kwargs_gpu = {
    "aucell": {},
    "gsea": {"times": 1},       # GPU times=1 时通常返回 raw ES，和 CPU NES 尺度可能不同
    "gsva": {},
    "ora": {},
    "ulm": {},
    "mlm": {},
    "waggr": {"times": 1},
    "zscore": {},
    "viper": {},
    "mdt": {},                  # 仅 GPU
    "udt": {},                  # 仅 GPU
}

# 4.1 先运行有 CPU 基线的方法
for method, kwargs in method_kwargs_cpu.items():
    score, elapsed = run_es_method(
        adata,
        method=method,
        engine=ENGINE_CPU,
        signatures=signatures,
        geneset_name=GENESET_NAME,
        tmin=TMIN,
        **kwargs,
    )
    scores_cpu[method] = score
    adata.obs[f"{method} (CPU)"] = score
    runtime_records.append({"method": method, "engine": "CPU", "seconds": elapsed})

# 4.2 再运行 GPU 方法，包括 MDT 和 UDT
for method, kwargs in method_kwargs_gpu.items():
    score, elapsed = run_es_method(
        adata,
        method=method,
        engine=ENGINE_GPU,
        signatures=signatures,
        geneset_name=GENESET_NAME,
        tmin=TMIN,
        **kwargs,
    )
    scores_gpu[method] = score
    adata.obs[f"{method} (GPU)"] = score
    runtime_records.append({"method": method, "engine": "GPU", "seconds": elapsed})

runtime_df = pd.DataFrame(runtime_records)
print(runtime_df.pivot(index="method", columns="engine", values="seconds"))


# ============================================================
# 5. UMAP 可视化：比较每个方法的 CPU/GPU 空间分布
# ============================================================
# 作用：把通路活性分数像基因表达一样投影到 UMAP 上。

# 示例：绘制 AUCell CPU 与 GPU 结果
ov.pl.umap(
    adata,
    color=["aucell (CPU)", "aucell (GPU)"],
    cmap="Reds",
    frameon=False,
)

# 示例：绘制所有 GPU 方法；正式项目中可按需拆分，避免图太拥挤
all_gpu_colors = [f"{m} (GPU)" for m in scores_gpu.keys()]
ov.pl.umap(
    adata,
    color=all_gpu_colors,
    cmap="Reds",
    frameon=False,
)


# ============================================================
# 6. CPU 与 GPU 一致性评估
# ============================================================
# 作用：对同时具有 CPU/GPU 结果的方法计算 Pearson、Spearman 与 RMSE/range。
# 注意：GSEA 可能因为 NES/raw ES 尺度差异导致数值不完全重合。

def rmse_over_range(x, y):
    """计算 RMSE 除以数值范围，避免不同方法尺度不可比。"""
    rmse = np.sqrt(np.mean((x - y) ** 2))
    value_range = np.max(x) - np.min(x)
    return rmse / (value_range + 1e-12)

comparison_records = []
for method in scores_gpu:
    if method not in scores_cpu:
        continue

    x = np.asarray(scores_cpu[method], dtype=float)
    y = np.asarray(scores_gpu[method], dtype=float)

    # 如果分数没有变化，相关系数可能为 NaN，这是 ORA 等方法中可能出现的正常现象
    pearson = pd.Series(x).corr(pd.Series(y), method="pearson")
    spearman = pd.Series(x).corr(pd.Series(y), method="spearman")

    comparison_records.append({
        "method": method,
        "pearson_r": pearson,
        "spearman_r": spearman,
        "rmse_over_range": rmse_over_range(x, y),
    })

comparison_df = pd.DataFrame(comparison_records)
print(comparison_df)


# ============================================================
# 7. 绘制 CPU/GPU 分数散点一致性图
# ============================================================
# 作用：将不同方法的分数标准化后放到一张图中，直观看是否接近 y=x 对角线。

methods_with_both = [m for m in scores_gpu if m in scores_cpu]
cmap = plt.get_cmap("tab20")

fig, ax = plt.subplots(figsize=(5.5, 5.5))
for i, method in enumerate(methods_with_both):
    x = scores_cpu[method]
    y = scores_gpu[method]

    # Z-score 标准化，便于不同方法放在同一尺度比较
    xz = (x - x.mean()) / (x.std() + 1e-12)
    yz = (y - y.mean()) / (y.std() + 1e-12)

    ax.scatter(
        xz,
        yz,
        s=4,
        alpha=0.4,
        label=method,
        color=cmap(i % cmap.N),
    )

ax.plot([-6, 6], [-6, 6], "--", color="grey", lw=0.8)
ax.set_xlabel("CPU score, z-scaled")
ax.set_ylabel("GPU score, z-scaled")
ax.set_title("CPU vs GPU enrichment score consistency")
ax.legend(markerscale=3, fontsize=8, frameon=False)
plt.tight_layout()
plt.show()


# ============================================================
# 8. 保存结果，便于后续项目复用
# ============================================================
# 作用：将打分后的 AnnData 和统计表保存到本地。
# 路径可按自己的项目目录修改。

adata.write_h5ad("pancreas_ov_es_scores_demo.h5ad")
runtime_df.to_csv("ov_es_runtime_records.csv", index=False)
comparison_df.to_csv("ov_es_cpu_gpu_consistency.csv", index=False)

代码使用建议#

快速看通路空间分布： 优先跑 aucell、gsva、zscore。
做 regulon/转录因子活性： 优先考虑 viper、ulm、mlm，并确认 signature 是否带权重。
做正式结果： 不要只比较一条通路，建议批量导入 Hallmark、GO BP、KEGG、Reactome 或自定义口腔疾病相关基因集。
做口腔生信项目： 可把 GENESET_NAME 替换为炎症、TNF/NF-κB、IL-17、EMT、hypoxia、angiogenesis、T cell exhaustion、macrophage polarization 等通路。
复现注意： 原文中的 CPU/GPU 耗时只对应作者环境和 300 细胞示例，不应直接作为正式 benchmark 结论。