一、写在前面

Perturb-seq将单细胞 RNA 测序 （scRNA-seq） 与基于 CRISPR 的扰动相结合，可实现大规模单细胞检测，研究单个细胞的基因功能和细胞反应。例如我们此前分享的风靡一时的scTenifoldKnk单细胞虚拟基因敲除,准不准?，便是其中的一种。本次推送分享发表在预印刊上题为"Benchmarking AI Models for In Silico Gene Perturbation of Cells"的文章，全方位得测评了主流AI虚拟细胞干扰工具。

当下面临的局限性：

基因和细胞类型数量大：人类有大约 20,000 个蛋白质编码基因和大约 400 种主要细胞类型，几乎不可能探索所有细胞类型的所有基因扰动，更不用说探索大量的基因组合（基因调控网络）了

某些细胞无法在体外长时间存活：这些实验通常需要多天的细胞孵育，但许多细胞类型无法在培养物中长时间存活

这些限制刺激了许多用于计算模拟基因扰动的人工智能模型的发展，这些模型预测基因扰动后的细胞状态（通常表示为基因表达谱）

当下没有对模拟扰动模型的统一评估方法：

模拟基因扰动的算法通常是针对于特定场景的，在有限的数据集上进行验证，并使用不一致的指标进行评估。例如，GEARS 旨在预测已知细胞类型中未知的扰动（unseen perturbation），而 scGen 和 CPA 则侧重于预测对未知细胞（unseen cell）进行已知扰动之后细胞的反应

这些算法通常适合在特定的数据集和指标上进行评估，例如，Dynamo 和 CellOracle 利用代表生物发育轨迹的数据集（例如造血数据）进行验证，而 GEARS 和 scGPT 使用单细胞扰动数据集进行验证。

目前出现的评估模拟基因扰动算法的方法大部分具有局限性，例如，Wu et al.20为各种机器学习方法提出了一个Perturation Benchmark，但它不评估算法预测未知扰动的能力。

目前缺乏一个全面的基准测试框架，提供标准化数据集和评估指标来公平地比较各个模拟基因扰动的算法。

 如果需要单细胞数据分析教学、生信热点全文复现、自测数据个性化分析辅导、实验科研服务和常态化实验学习，欢迎联系[Biomamba_zhushou]。

二、文章工作

文章中整理并过滤了 scPerturb 的基因扰动数据集，保留了 17 个单细胞数据集，包括 984,000 个细胞（cellcount）和 3,190 个扰动（Number of perturbation）

本研究设计了四种关键评测场景，全面评估模型性能：

1. 未见过的扰动预测（Unseen perturbation transfer）： 评估模型对已知细胞类型中未见过的基因扰动效应的预测能力。研究团队从scPerturb数据库筛选整理了17个单细胞数据集，包含984,000个细胞和3,190个基因扰动。包括 CRISPR-Cas9 敲除、CRISPR干扰 （CRISPRi） 和 CRISPR 激活（CRISPRa）。

2. 未见过的细胞类型预测（Unseen cell type transfer）： 测试模型将已知扰动效应迁移到新的细胞类型的能力。研究使用了来自K562和RPE1两种细胞系的490个共同扰动数据。来自 K562 细胞系的数据集用于训练，而 RPE1 数据集用作测试数据。

3. 零样本迁移预测（Zero-shot transfer）： 评估模型在没有针对性训练数据的情况下，预测CMAP框架下的基因扰动效应的能力。研究收集了6种细胞系的21,361个扰动数据。

4. 细胞状态转换预测（Cell state transition prediction）： 考察模型预测关键基因驱动细胞状态转换的能力。团队选取了4个代表不同生物过程的数据集进行评测。

Fig.1

考虑到只有一部分基因受到特定扰动的影响，本研究中将每个扰动的基因表达谱中的基因分为三组：

（1） 差异表达基因 （DEG），它们在扰动前和扰动后细胞之间表现出显著的表达差异

（2） 稳定高表达基因（SHG），不差异表达，在扰动前后细胞中均表现出高表达

（3） 稳定低表达基因 （SLG），不差异表达，在扰动前后细胞中均表现出低表达

对于差异表达基因，我们使用 PCC （DEG_PCC）和 ACC （DEG_ACC） 值评估性能，而对于稳定高表达和低表达基因集，我们使用 MAE 值评估性能，分别称为 SHG_MAE 和 SLG_MAE

研究团队设计了一套全面的评测指标：

• 差异表达基因的皮尔逊相关系数（DEG_PCC）：衡量预测值与真实值的相关性

ΔX：真实表达量；Δ~X：预测值；μ：真实表达量的平均值；μ~：预测值的平均值；σ：真实表达量的标准差；σ~：预测值的标准差

（这个公式可能出现了错误，我认为其中E应该替换为∑求和符号）

• 差异表达基因的方向准确率（DEG_ACC）：评估预测变化方向的准确性

在公式(5)中，分子的含义是将II(x)的计算结果求和，而II(x)的值在X=0时为1，其余情况II(x)=0

在数学和计算机运算中，sign()的功能是取某个数的符号（正或负）：

当x>0，sign(x)=1;

当x=0，sign(x)=0;

当x<0， sign(x)=-1；

当方向准确时II(x)之内的数字应该等于0（1-1或者(-1)-(-1)），其余情况均为方向不准确

所以ACC就是一个介于0和1之间的分数，ACC越大，预测的准确性越高

• 稳定高表达/低表达基因的平均绝对误差（SHG_MAE/SLG_MAE）：考察模型对非差异表达基因的预测准确度

这个公式的意义就是取预测值绝对误差的平均值，以此衡量预测的准确性

模型简述：

Linear（存疑）：是一种线性模型，但文章method部分对其计算方法的描述不够清晰

Linear-scGPT：将训练scGPT的1000万个细胞的单细胞数据作为矩阵嵌入到Linear模型中，可以视为经过了大量数据的预训练

CellOracle：基于染色质可及性数据（ATAC-seq）的正则化线性模型，利用scATAC-seq数据构建基因调控网络，利用线性模型模拟扰动效应

GEARS：

基因嵌入与扰动嵌入：对于每个基因，GEARS初始化一个基因嵌入向量和一个基因扰动嵌入向量。

图神经网络（GNN）：利用基因共表达图和基因本体（Gene Ontology, GO）图，通过图神经网络整合基因之间的关系信息。

扰动预测：将基因嵌入与扰动嵌入结合，通过跨基因层和基因特异性输出层，预测扰动后的基因表达。

scGen：scGen通过变分自编码器（VAE）学习细胞基因表达数据的低维潜在表示。它利用潜在空间中的向量差来模拟扰动效应。

CPA：自动编码器

CellOT：CellOT（Cell Optimal Transport）是一种基于神经最优传输理论的方法，用于学习和预测单细

胞扰动响应。其核心在于通过最优传输理论建立未处理细胞状态与扰动后细胞状态之间的映射关系，从而预测细胞在特定扰动下的响应

scGPT：scGPT（single-cell Generative Pre-trained Transformer）是一种基于生成式预训练Transformer架构的单细胞多组学基础模型，旨在通过大规模单细胞数据的预训练，捕捉基因和细胞之间的复杂关系，为多种下游任务提供强大的基础

scFoundation：transformer架构，与scGPT类似，有以下不同

• scGPT：

强调Transformer架构的灵活性（scGPT可以将不同组学层的特征视为独立的token，从而实现多模态数据的整合。例如，它能够同时处理转录组和表观组学数据，捕捉其他方法未能识别的细胞亚群）。

提供了丰富的下游任务支持，包括跨物种建模。

• scFoundation：

更注重“开箱即用”的零样本学习能力，无需微调即可在多种任务中表现良好。

模型规模较大（如1亿参数），在处理高维数据时表现出色。

Correlation 方法：首先构建一个基因-基因相关矩阵，其中每个元素代表两个基因之间的 PCC。当一个基因受到干扰时，其他基因的表达会根据相应 PCC 确定的比例因子而变化。如果基因被敲低，则其表达直接设置为 0。如果基因过表达，则在其原始标准化表达值上增加 1。

在相同细胞中进行不同的扰动

经过严格的数据过滤和标准化处理（Methods），我们纳入了17 个 scRNA-seq 数据集进行分析。

本文计算了各个纳入数据集的E距离（E-dis），这是一种点云之间距离的统计度量。

E-distance是一个统计学指标，用于度量两个点云分布之间的距离。在基因扰动场景下，它用来衡量扰动前后细胞状态的变化程度。

Fig.2

本次测试中，Linear、Linear-GPT、CellOracle、GEARS、scGPT 和 scFoundation共六种模型进行测试。

并且加入了一种基础方法——KnownAverage，该方法在所有已知扰动中对所有细胞的基因表达量取平均值，作为对未知扰动的预测。

值得注意的是，KnownAverage 方法在所有四种类型的指标中始终表现出一些最佳的整体性能。

对于 SHG_MAE 和 SLG_MAE，KnownAverage 方法在几乎所有数据集中都再次取得了最佳性能，并且Linear-scGPT 在大多数情况下始终排名前三。尽管 GEARS、scGPT 和 scFoundation 等深度学习方法并没

有始终获得最高排名，但它们的指标并没有明显差于表现最好的方法。

在 17 个数据集中模型的 DEG_PCC 排名热图上进行了分层聚类（不包括 CellOracle，因为它与所有数据集不兼容）

结果揭示了不同的性能模式，将数据集分为三个集群：G1、G2 和 G3

G3 组表现出最高的平均扰动数 （527.8） 和平均 E 距离 （24.2），表明这些数据集包括广泛的扰动，在扰动前具有显着的表型效应。单细胞基础模型，如 scGPT 和 scFoundation，在 G3 数据集上表现最佳，这可能是由于训练数据中存在丰富的异构扰动，这支持了大型深度学习模型的学习能力。

G2 组显示出中等的平均扰动数 （87.2） 但较低的平均 E 距离 （2.2），这表明这些数据集中的扰动仅引起表达谱的细微变化。有趣的是，scFoundation 在 G2 数据集上取得了最佳性能，证明了它能够比其他方法（包括 scGPT）更有效地捕获轻微的扰动效应。

G1 组的平均扰动数最低 （25.2）。在这些数据集中，KnownAverage、Linear 和 Linear-scGPT 等简单方法的性能优于深度学习方法。这种性能差距可能归因于扰动数量有限，这限制了深度学习模型有效学习扰动效应的能力。

该测试的结果

scGPT 和 scFoundation成为预测未知扰动结果任务中表现最好的方法。有了足够的扰动数据进行训练，他们能够准确预测 DEGs 的表达变化，同时保留非 DEGs 的表达水平。

有趣的是，最简单的方法 AverageKnown 在各种数据集中也表现得非常出色。通过使用聚类热图可视化该数据集中扰动的平均表达变化，我们观察到超过一半的扰动基因被分组到一个簇中（图 S3）（文章中无fig.S3，无Extended data figures and tables，无Supplementary information，下载PDF后依然无fig.S3），因此，AverageKnown 方法可以通过简单地平均这些受扰动的单元来获得良好的性能。

将相同扰动应用于不同细胞中

本研究使用了同一研究中的 Replogle2022_K562 和 Replogle2022_RPE1 数据集，其中K562数据集用于训练模型，RPE1数据集用于测试模型。

K562，一种源自慢性粒细胞白血病患者的淋巴母细胞系，以及 RPE1，一种视网膜色素上皮细胞系。

经过标准化数据处理，我们在两个数据集之间确定了 490 个共享的基因扰动和 3,144 个重叠的基因。

Fig.3

本文在这个任务中选择了 scGen 、 CPA 、 CellOT 、 scGPT 和 scFoundation这几个模型进行测试。

此外，研究者还加入了一个简单模型：DirectTransfer，它直接将扰动的基因表达变化从Replogle2022_K562 转移到 Replogle2022_RPE1（Xa-Xa perturb=Xb-Xb perturb）。

在评估模型性能之前，该研究首先评估了两个数据集之间共享扰动的相似性。对于每个数据集中的每次扰动，我们计算了受扰动细胞和未受扰动细胞之间平均基因表达谱的基因表达改变，计算了Replogle2022_K562 和 Replogle2022_RPE1 中每次扰动的基因表达改变（Alter_PCCs）的PCC（皮尔逊相关系数）（fig3b）。简单的来说，就是两个数据集中Xa-Xa perturb与Xb-Xb perturb的相似性（a与b代表两个细胞系）。

在图中Alter_PCCs 的分布呈现双峰模式：大部分扰动聚集在 0.05 附近，表明扰动对两个细胞系之间基因表达的影响不同，而小部分扰动聚集在 0.45 附近，表明扰动对两个细胞系可能具有相似的影响（PCC值越低，表明两个对象之间越不呈线性相关）。结论是Xa-Xa perturb与Xb-Xb perturb大部分不相似。

fig.3c log3(E-dis)，用于度量两个点云分布之间的距离。在基因扰动场景下，它用来衡量扰动前后细胞状态的变化程度。研究者根据log3(E-dis)将 490 个共享扰动分为四组：P1（在两个数据集中均高表达）、P2（在K562 中高表达但在 RPE1 中低表达）、P3（在 K562 中低表达但在 RPE1 中高表达）和 P4（在两个数据集中都低表达）。这些分类反映了扰动的不同影响，这些影响与不同的Alter_PCC值相关联。

测试结果

DirectTransfer 在各种指标（包括 DEG_PCC、DEG_ACC、SHG_MAE 和 SLG_MAE）上都优于所有其他方法。（fig.3e）

P1：DirectTransfer 在 P1 组中取得了最高性能。CellOT 和 CPA 在 DEG_PCC 和 DEG_ACC 方面表现出中等性能，而 CellOT 在 SHG_MAE 和 SLG_MAE 方面表现出下降。相比之下，scGPT 、 scFoundation 和 scGen的整体性能较差。

基因本体论 （GO） 富集分析显示，P1 扰动与 mRNA 加工和 RNA 剪接等基本生物过程相关，属于种属特异性比较低的基因，解释了它们在 K562 和 RPE1 细胞系中显着且一致的影响。（fig.3g）

P2：对于 P2 中的扰动，scFoundation 在 DEG_PCC 中实现了接近 DirectTransfer 的性能，而 scGen、CPA、CellOT 和 scGPT 产生了极低的 DEG_PCC 值。（fig.3f）

P2 中的扰动仅与 5 个基因相关，例如 GAB2，它参与白细胞和肥大细胞脱颗粒的正调节。这表明 P2 中的基因可能在 K562 中发挥重要作用，但在 RPE1 中作用较弱。（fig.3g）

P3：在 P3 扰动上，scGPT 在所有指标上的表现与 DirectTransfer 相当，而其他算法的性能明显较差。

P3 中的扰动在一些与蛋白质生产和分泌相关的生物过程中富集，这可能反映了 RPE1 细胞的上皮性质，RPE1 细胞参与膜和分泌蛋白的更多加工和运输。

P4：DirectTransfer 在 P4 扰动上仍然保持了最好的整体性能，其次是 scGPT 和 scFoundation。

P4 中的扰动对应于对 K562 和 RPE1 细胞影响较小基因。这些扰动在一些在细胞中不占主导地位的生物过程中富集，例如线粒体翻译和线粒体基因表达。

上述所有结果表明，将相同扰动应用于不同细胞的任务中，高级机器学习或统计方法未能胜过最简单的模型DirectTransfer。

基于 transformer 的模型（如 scGPT 和 scFoundation）的性能通常优于基于自动编码器的方法（如scGen、CPA 和 CellOT），但仍明显低于 DirectTransfer 的性能。但 DirectTransfer 完全忽略了不同细胞类型的特异性，不适合实际使用。

零样本迁移预测

我们使用 CMAP 数据作为基准测试的测试数据。我们收集了单细胞数据集与CMAP数据重叠的6个细胞系：A375、PC3、MCF7、HT29、HepG2 和 A549。验证了单细胞数据中基因表达谱与相应CMAP数据之间的一致性。（fig.4b）

对于 CMAP 数据集中的每个扰动，5 级数据提供了一个聚合的 z 分数向量，表示基因表达变化。正 z 值表示扰动后表达增加（上调），而负 z 值表示表达降低（下调）。为了改进分析，我们选择了 z 分数大于 1.5 或小于 -1.5 的基因作为显著影响，并排除了单细胞数据集中表达水平小于 0.1 的基因。这些基因被确定为最高显著基因，我们使用这些基因计算方向准确性（ACC）以评估模型性能。

Fig.4

这些单细胞数据集被处理为模拟扰动算法的输入，以直接预测由基因扰动引起的基因表达变化，而无需任何额外的训练。使用方向精度 （ACC） 评估模型性能，通过将预测结果与 CMAP 数据的真实数据进行比较来计算。该分析总共包括了 6 个细胞系的 21,361 次扰动。

Correlation 方法：使用基因之间的 PCC 作为模拟扰动的相关变化因子;当基因受到干扰时，相对于它具有较高 PCC 的基因受到更显着的影响。

Method中对Correlation的描述：首先构建了一个基因-基因相关矩阵，其中每个元素代表两个基因之间的PCC。 当一个基因受到干扰时，其他基因的表达量会根据相应的 PCC 确定的比例因子发生变化。 如果一个基因被敲除，其表达量会直接设为 0。

DirectTransfer：利用了两个来源：K562 细胞 （DirectTransfer_K562） 的基因组规模 Perturb-seq 数据，以及来自 L1000 平台的数据

Random：基因随机上调或下调

我们在单细胞数据中观察到 CMAP 的重要变化基因有一种独特的表达模式：CMAP 中的上调基因在单细胞数据中往往表现出低表达水平，而 CMAP 中的下调基因在单细胞数据中往往表现出高表达水平。 如图 4d 和图S4 所示，在所有分析的细胞系中，这一趋势是一致的。

基于这一观察结果，研究者引入了平均预测法（AverageSC）（存疑），该方法只需对单细胞数据中所有关键基因的表达值取平均值即可生成预测结果。

原文：Building on this observation, we introduced an AverageSC method for this scenario, which simply averages the expression values of all landmark genes in single cells to generate predictions.

结果

CellOracle 无法预测所有细胞系中基因表达变化的方向。同样，CPA 表现不佳，在某些细胞系（如HepG2）中，结果略差于 Random。Correlation 方法显示比 Random 略有改善。在 DirectTransfer 方法中，DirectTransfer_K562 在 Random 级别执行，而 DirectTransfer_L1000 方法实现了最高的总体准确度，平均值为 0.692。这表明从同一平台（例如 L1000 平台）转移基因表达变化可以提供可靠的参考。

总体而言，只有 DirectTransfer_L1000 和 AverageSC 方法实现了明显优于 Random 的性能。但是，DirectTransfer_L1000 依赖于来自同一平台的数据，而这些数据在实际应用中通常不可用，从而限制了其实用性。同时，仅平均基因表达值的 AverageSC 方法不考虑基因扰动机制，因此对于理解或建模扰动效应不切实际。大多数其他方法无法预测重要基因的变化方向，这凸显了当前方法在零镜头情景下预测扰动效应的局限性。虽然单细胞基础模型显示出一些超过随机性能的潜力，但它们在单个细胞系上的最佳准确度约为0.55，仅略高于随机，仍远未实际适用性。

细胞状态转换预测

许多生物过程，例如器官发育和细胞重编程，是由单个基因或一小部分关键基因驱动的。评估计算方法是否能够准确预测受这些基因影响的细胞状态转换，对于评估它们对研究人员自身数据的适用性和指导实验设计至关重要。

为了建立基准框架，我们收集了四个代表不同生物过程的单细胞数据集：

Sun202229验证 PDCD1 的敲除可以促进耗竭的 T 细胞的恢复，恢复其效应器和记忆功能；

Ainciburu202330验证 SPI1 和 GATA1 作为驱动造血细胞分化的经典双基因切换基序，在这个研究中SPI1被证实可促进干细胞向粒细胞-单核细胞祖细胞 （GMP）谱系的分化，GATA1可促进干细胞向巨核细胞-红细胞祖细胞 （MEP） 谱系的分化；

Nair202331用于验证成纤维细胞可以通过四种转录因子（TF）重编程为干细胞：POU5F1、SOX2、KLF4 和MYC；

Steele202032用于验证 PTF1A 的敲除可诱导腺泡细胞转化为导管样细胞（腺泡管化生，ADM）。

并非所有方法都可以应用于所有四个数据集。CellOracle 仅限于具有已知结合基序的 TF 的计算机扰动，而Sun2022 数据集中的 PDCD1 不包括在其 TF 列表中。CPA 和 DirectTransfer 依赖于具有匹配扰动的参考数据，这些数据对于 Nair2023 数据集和 Steele2020 数据集中的 PTF1A 不可用

Fig.5

Sun2022 和 Steele2020 数据集中的细胞状态表现出离散的细胞状态，而 Ainciburu2023 和 Nair2023 中的细胞形成连续的轨迹。想要准确预测具有离散细胞状态的数据集中的转换，需要模型在基因扰动后捕获细胞状态的更实质性变化。在流线图中，大多数方法无法捕捉到正确的过渡趋势

在 Ainciburu2023 数据集中，Correlation 方法和 CellOracle 都成功预测了 GATA1 敲除后从 MEP 到 LMPP 的转变，以及 SPI1 敲除后从 GMP 到 LMPP 的转变（fig.5d-e）。然而，其他方法未能捕获这些特定的谱系转换。在 Nair2023 数据集中，CellOracle 在流线图中提供了从成纤维细胞样细胞到干细胞样细胞转变的最准确可视化（fig.5f）。scGPT 在绘图中也表现出准确的趋势，尽管它不太明显。相比之下，Correlation 方法出现了错误，认为成纤维细胞样细胞转变为角质形成细胞样细胞（正确结果应该为干细胞样细胞）。