【AI4Sci】预测高通量筛选中对复杂干扰的细胞反应

摘要

近期多重单细胞转录组学实验的进展，促进了药物和基因扰动的高通量研究。然而，穷尽式地探索组合扰动空间在实验上是不可行的。因此，需要计算方法来预测、解释和优先排序扰动。在此，我们提出组合扰动自编码器（compositional perturbation autoencoder，CPA），它结合了线性模型的可解释性与深度学习方法的灵活性，用于单细胞响应建模。CPA 能在计算机中（in silico）预测单细胞水平的转录扰动响应，适用于未见过的剂量、细胞类型、时间点和物种。利用新生成的单细胞药物组合数据，我们验证了 CPA 能预测未见过的药物组合，且优于基线模型。此外，该架构的模块化设计允许引入药物的化学表示，从而能够预测对完全未见过的药物的细胞响应。此外，CPA 也适用于基因组合筛选。我们通过在单细胞 Perturb-seq 实验中计算机填补 5,329 个缺失组合（占所有可能性的 97.6%）来证明这一点，其中包含多样的基因互作。我们设想 CPA 将通过在单细胞水平进行计算机响应预测，促进高效的实验设计和假设生成，从而加速单细胞技术的治疗应用。

• multiplexed single-cell transcriptomics：多重单细胞转录组测序，比如 sci-Plex、Perturb-seq，可以一次性测量成千上万个细胞在不同扰动下的基因表达。

• drug and genetic perturbations：药物扰动（不同化合物、剂量）和基因扰动（CRISPR 敲除、过表达）。

• high-throughput：高通量，即一次实验能测很多条件。

• 原始 CPA 只能预测训练集中出现过的药物。chemCPA 用分子指纹/化学结构编码器替代药物嵌入字典，可以预测全新药物。

• 不只是药物，基因扰动（CRISPR 组合）也能用 CPA 建模。

• Perturb-seq：CRISPR 筛选 + 单细胞 RNA 测序。

引言

单细胞 RNA 测序（scRNA-seq）能够刻画跨组织、跨物种的百万级细胞的基因表达。近期，新技术将这些测量扩展到高通量筛选（HTS），可以测量细胞对数千种独立扰动的响应。这些进展有望促进和加速药物开发。在单细胞水平应用 HTS，既能提供全面的分子表型，又能捕捉异质性响应——这是传统 HTS 无法识别的。

• scRNA-seq：单细胞转录组测序，测量每个细胞中所有基因的表达量。

• 引用的是 Tabula Muris（小鼠细胞图谱）、Domcke（人类胚胎细胞图谱）、Han（人类细胞图谱）等大型图谱项目。

• HTS（high-throughput screens）：高通量筛选，一次实验测很多条件。

• comprehensive molecular phenotyping：全面分子表型，即测所有基因的表达，而不是只看几个指标。

• heterogeneous responses：异质性响应，即同一扰动下不同细胞的反应可能不同。

• traditional HTSs：传统高通量筛选通常只测细胞存活率或少数蛋白，看不到异质性。

虽然“细胞哈希（cellular hashing）”等高通量方法能以低成本实现多样本 scRNA-seq，但这些策略需要昂贵的文库制备，且难以扩展到大量扰动。这些缺点在探索组合疗法或基因扰动时更加明显，因为对所有可能组合进行实验筛选是不可行的。

• cellular hashing：用脂质标签或寡核苷酸标签给不同样本的细胞“贴条形码”，混在一起测序，降低成本。

虽然人类细胞图谱（Human Cell Atlas）等项目旨在以可重复的方式全面绘制跨组织的细胞状态，但为扰动对基因表达的影响构建类似图谱是不可能的，因为可能性太多。既然暴力搜索组合空间不可行，就有必要开发计算工具来指导组合扰动空间的探索，从而提名有前景的候选组合疗法。一个成功的组合空间导航计算方法，必须能够预测细胞在新型扰动组合下的行为——这些扰动在原始实验中只被单独测量过。这些数据被称为分布外（OOD）数据。OOD 预测将使得研究以下情况成为可能：不同治疗剂量、组合疗法、多基因敲除、跨时间变化。

• Human Cell Atlas：人类细胞图谱，目标是绘制所有人类细胞类型。

• perturbation atlas：扰动图谱，即“所有扰动 × 所有细胞 × 所有剂量”的响应地图。

• OOD（Out-Of-Distribution）：训练集分布之外的数据。

  • 列举 OOD 预测的四类应用场景：

    1. 不同剂量：训练时只有低剂量，预测高剂量。

    2. 组合疗法：训练时只有单药，预测组合。

    3. 多基因敲除：训练时只有单敲除，预测双敲除。

    4. 跨时间：训练时只有早期时间点，预测后期。

近期，已经开发出若干计算方法来预测细胞对扰动的响应。

• 第一类方法利用机制建模，用来预测细胞存活率或少数选定蛋白的丰度。虽然机制模型在解释相互作用方面很强大，但它们通常需要纵向数据（实际中往往没有），且大多数无法扩展到全基因组测量，无法预测高维 scRNA-seq 数据。
• 线性模型没有这些可扩展性问题，但预测能力有限，且无法捕捉非线性的、细胞类型特异的响应。
• 相比之下，深度学习（DL）模型没有这些局限。近期，DL 方法已被用于建模 scRNA-seq 数据的基因表达潜在空间，以及描述和预测单细胞响应。然而，当前基于 DL 的方法也有局限：它们只能建模少量扰动；无法处理组合处理；无法纳入连续协变量（如剂量、时间）或离散协变量（如细胞类型、物种、患者）。因此，虽然当前 DL 方法已经能建模单个扰动，但还没有一个是为高通量筛选（HTS）设计的。

• mechanistic modeling：机制建模，基于生化反应网络、信号通路等已知机制。

  • cell viability：细胞存活率，如药物处理后细胞死没死。

  • abundance of a few selected proteins：少数选定蛋白的丰度，如磷酸化蛋白水平。

  • 这类模型只预测少数几个指标，不是全基因组表达。

  • interpret interactions：解释扰动之间的相互作用（如协同、拮抗）。

  • longitudinal data：纵向数据，即同一细胞/样本在不同时间点的多次测量。实际实验中很难获得。

• Linear models：线性模型

  • do not suffer from scalability issues：线性模型计算简单，能扩展到全基因组。

  • limited predictive power：预测能力有限，只能拟合线性关系。

• 深度学习（DL）模型

  • a handful of perturbations：少量扰动，通常只有几种。

  • combinatorial treatments：组合处理，即药物/基因组合。

  • continuous covariates：连续协变量，如剂量、时间。

  • discrete covariates：离散协变量，如细胞类型、物种、患者。

方法类型	代表工作	优点	缺点
机制模型	Fröhlich 2018, Yuan 2021	可解释相互作用	需纵向数据、无法扩展、无法预测高维
线性模型	Dixit 2016, Kamimoto 2023	可扩展	预测能力有限、无法捕捉非线性
深度学习	scGen, trVAE, DrVAE	可扩展、非线性	只能少量扰动、无法组合、无法纳入协变量

在此，我们提出组合扰动自编码器（CPA），一种用于预测 scRNA-seq 扰动响应的方法，涵盖剂量、时间、药物、基因敲除等多种条件的组合。CPA 借鉴了可解释线性模型的思想，并将其应用于灵活的深度学习模型中，以学习扰动和协变量的因子化潜在表示。给定一个 scRNA-seq 数据集、所施加的扰动、以及描述实验设置的协变量，CPA 将数据分解为一组嵌入（表示），分别对应细胞类型、扰动和其他外部协变量。借助对抗损失，这些嵌入彼此独立，因此可以在预测时重新组合，以预测新的扰动-协变量组合的效果。因此，通过探索新组合，CPA 可以指导实验设计，将假设引向实验学家感兴趣的表达模式。我们在六个公开数据集和一个新的非小细胞肺癌（A549）数据集上展示了 CPA 的实用性，该数据集包含 32 种单药和组合药物扰动，涵盖多项任务，包括对化合物、剂量、时间序列信息和基因扰动的响应预测与分析。

• autoencoder：自编码器，一种神经网络结构（编码器 + 解码器）。

• across combinations of conditions：跨条件组合，比如“药物 A + 剂量 1μM + 时间 24h + 细胞类型 A549”。

• borrows ideas from interpretable linear models：借鉴线性模型的可分解性（加性结构）。

• flexible DL model：灵活的深度学习模型，能捕捉非线性。

• factorized latent representations：因子化潜在表示，即把数据分解成多个独立的“因子”（扰动因子、协变量因子等）。

• Given...：输入包括三部分：

  1. scRNA-seq 数据（基因表达矩阵）

  2. 扰动标签（用了什么药、敲了哪个基因）

  3. 协变量（细胞类型、物种、时间、剂量等）

• adversarial loss：对抗损失，即判别器与编码器对抗，迫使基础状态嵌入不包含扰动/协变量信息。

• guide experimental design：CPA 预测“药物 A + 药物 B”会产生强烈的免疫激活表达模式，实验学家就可以优先做这个组合。

CPA 的核心思想图解：

输入：
  scRNA-seq 数据 + 扰动标签 + 协变量（细胞类型、剂量、时间等）
        ↓
编码器 + 对抗损失
        ↓
分解为独立嵌入：
  [基础状态嵌入] + [扰动嵌入] + [协变量嵌入]
        ↓
预测时自由重组：
  [细胞类型 A 的基础状态] + [药物 B 的嵌入] + [剂量 1μM 的嵌入]
        ↓
解码器
        ↓
输出：预测的基因表达

CPA 完整工作流程图解

输入：
  x_i（基因表达）+ d_i（扰动/剂量）+ c_i（协变量）

编码器 f_enc：
  x_i → z_i^basal（基础状态）

对抗训练：
  判别器试图从 z_i^basal 预测 d_i, c_i
  编码器努力让判别器猜不出
  → z_i^basal 不含扰动/协变量信息

组合：
  z_i = z_i^basal + V^perturbation · f(d_i) + V^covariate · c_i
        ↑              ↑                    ↑
     基础状态      扰动嵌入（非线性缩放）   协变量嵌入

解码器 f_dec（非线性）：
  z_i → x_i'（预测的基因表达）

预测新组合：
  取 z_i^basal + 新的 V^perturbation · f(d') + c_i
  → 解码 → 预测新组合的基因表达

结果

多重扰动作为基因表达潜在空间中的组合过程

以往的工作把扰动对基因表达的影响建模为相互独立的过程。差异表达分析将每个条件与对照分别比较，而用条件变分自编码器（conditional variational autoencoder，cVAE）建模联合潜在空间则高度不可解释，且不适合预测条件组合的效应。我们的目标是将神经网络的潜在空间因子化，使其成为可解释的、组合性的模型。如果潜在空间是线性的，我们就可以把观察到的基因表达描述为一个因子模型，其中每个成分是一个单独的扰动。

• 差异表达：逐个比较，无法组合。如“药物 A vs 对照”，逐个条件单独比较

• cVAE：联合建模（把所有条件混在一起建模），但不可解释（因为潜在空间里的维度没有明确含义）、不能组合（因为 cVAE 没有把扰动分解成独立因子）

• CPA 目标：

  • factorize：因子化，把潜在空间分解成多个独立的“因子”。

  • interpretable：可解释，每个因子对应一个明确含义（如“药物 A 的效应”）。

  • compositional：组合性，因子可以自由组合，预测新组合。

• 如果潜在空间是线性的，那么“药物 A + 药物 B”的效应 = “药物 A 的效应” + “药物 B 的效应”。这样就能简单相加，预测组合。

然而，基因表达潜在空间，尤其在复杂组织中，是非线性的，最好用图或非线性嵌入来近似描述。在 scRNA-seq 数据集中，细胞群体的基因表达谱通常在多种扰动（如药物、基因敲除、疾病状态）下观察，并带有标注的协变量（如细胞系、患者、物种）。每个细胞都标注了其实验条件和扰动，其中实验协变量用分类标签表示，扰动用连续值表示（例如，药物以不同剂量施加）。这假设每个条件有足够多的细胞，以便用大型神经网络估计对照和扰动状态下的潜在空间。

我们不在基因表达空间中假设因子模型，而是在非线性潜在空间中建模扰动效应的非线性叠加，其中我们约束叠加是加性的（见材料与方法）。我们解耦扰动和协变量的效应，并通过将其编码为非线性变换的标量权重（即学习到的药物-响应曲线）来允许连续效应（如药物剂量）。线性潜在空间因子模型使得我们能够解释这个空间，通过将协变量驱动的潜在空间方差与每个扰动引起的方差分开。在评估时，我们不仅能插值并解释已观察到的扰动组合，还能预测其他组合，甚至可能是在不同的协变量设置下。

• categorical labels：分类标签，如细胞类型 = A549 / MCF7 / K562。

• continuous value：连续值，如剂量 = 10nM / 100nM / 1μM / 10μM。

• decouple：解耦，把扰动和协变量的效应分开。

• nonlinearly transformed scalar weight：非线性变换的标量权重，即 f(d)。

• learned drug-response curve：学习到的药物-响应曲线，如 log-sigmoid。

我们介绍了CPA（见材料与方法），该方法结合了自然语言处理和计算机视觉中的思想，来预测扰动组合对单细胞基因表达的影响。给定一个包含多种扰动和协变量的单细胞数据集，CPA 首先用编码器神经网络将细胞的基因表达分解为一系列可学习的、加性的嵌入，分别对应其基础状态、观察到的扰动和观察到的协变量。关键的是，CPA 编码器学习到的细胞基础状态的潜在表示，与扰动和协变量对应的嵌入是解耦的（不包含关于它们的信息）。这种解耦是通过训练一个判别器分类器，与 CPA 的编码器网络竞争来实现的。CPA 中编码器网络的目标是学习一个表示细胞基础状态的嵌入，使得判别器网络无法从中预测出扰动或协变量的值。为了考虑连续的时间或剂量效应，每个扰动的学习嵌入通过一个神经网络进行非线性缩放，该网络接收每个细胞的连续协变量值（如时间或剂量）。在将细胞基础状态、扰动和协变量的学习嵌入线性整合为一个统一嵌入后，CPA 使用非线性神经网络解码器恢复细胞的基因表达向量。解码器的非线性使得模型能够捕捉复杂的、细胞类型特异的、非加性的组合处理效应。

里不能有扰动或协变量的信息；如果混入了扰动信息，那么预测新组合时，基础状态就被“污染”了，无法自由组合。

考虑一个简单的例子：一个细胞受到两个基因敲除扰动。CPA 通过线性组合每个扰动的单扰动嵌入和基础状态，来重建组合处理的整体基因表达效应。这种约束使模型能够学习单个处理在与其他处理组合时的行为模式，从而预测训练中未见过的组合。当模型观察到某个单扰动与其他扰动组合的多样化训练数据时，该扰动组合的预测性能可能会提高。相反，当模型从未见过某个单扰动的任何组合训练数据时，它可能会对包含该处理的组合生成虚假预测。

• diverse training data：多样化的训练数据。

  • 例如：如果训练集中有“药物 A + 药物 B”“药物 A + 药物 C”“药物 A + 药物 D”，那么预测“药物 A + 药物 E”会更准。因为模型学到了“药物 A 在组合中的行为模式”。

• encoder neural network：编码器，把高维基因表达压缩到低维潜在空间。

• learnable, additive embeddings：可学习的加性嵌入。

• basal state：基础状态，细胞未受扰动时的状态。

• observed perturbation：观察到的扰动（如药物、基因敲除）。

• observed covariates：观察到的协变量（如细胞类型、物种）。

• 对抗训练的目标

  • 编码器的目标：让  不含扰动/协变量信息。

  • 判别器的目标：尽力从 中猜出扰动/协变量。

  • 两者对抗，最终编码器胜出， 变得“干净”。

与许多神经网络模型类似，CPA 使用反向传播，基于重建误差和判别器误差进行训练，以调整编码器网络、解码器网络、每个扰动和协变量的嵌入、以及剂量/时间非线性缩放器的参数。学习到的嵌入允许我们衡量不同扰动和协变量之间的相似性，即它们在对基因表达的影响上的相似程度。CPA 的主要特点是它在评估（预测）时的灵活性。在获得与某些观察到的基因表达、扰动和协变量值对应的解耦嵌入后，我们可以进行干预，将扰动嵌入替换为我们选择的任何其他扰动嵌入。这种操作实际上是在估计一个反事实问题的答案：如果这个细胞接受了不同的处理，它的基因表达会是什么样？这种方法在预测未见过的扰动组合及其对基因表达的影响方面特别有用。CPA 还可以可视化与扰动和协变量相关的转录相似性和不确定性，如后文所示。

• backpropagation：反向传播，神经网络训练的标准方法，通过梯度下降更新参数。

• reconstruction error：重建误差，即“输入 x_i”和“解码器输出 x_i'”之间的差距，越小越好。

• discriminator error：判别器误差，判别器试图从 z_basal 中预测扰动/协变量，它的误差越小，说明 z_basal 里的信息越“泄露”。

• transcriptional similarity：转录相似性，即扰动之间在基因表达影响上的相似程度。

CPA允许单细胞扰动实验的预测性和探索性分析

我们首先在三个小型单细胞数据集上展示了 CPA 的性能和功能：

1. PBMC + IFN-β 刺激（Kang et al, 2018）

2. 人肺癌细胞 + 4 种药物扰动（Srivatsan et al, 2020）

3. 跨物种、纵向的 LPS 处理吞噬细胞（Hagai et al, 2018）

（详见材料与方法）

这些数据集代表了该模型的不同潜在应用场景：

1. 不同细胞类型中的二元扰动，即“有刺激 / 无刺激”（IFN-β 加或不加）。

2. 多种剂量，如 sci-Plex 的 7 种剂量。

3. 多个物种 + 时间变化（而非剂量变化）

我们将每个数据集划分为三组：

1. train（训练集）：用于模型训练

2. test（测试集）：用于调整模型参数

3. OOD（分布外集）：训练和调参时从未见过，用于衡量模型的泛化能力

• 如果只用 train/test，无法判断模型是真的学会了“组合规律”，还是只是“背下了训练数据”。

• OOD 集强制模型预测完全没见过的条件（如最高剂量、某组合），这样才能证明模型的真正泛化能力。

在这里，我们考虑了来自狼疮患者样本的 PBMC，这些细胞被 IFN-β 处理。在这种情况下，刺激是二元的，没有任何连续协变量（如剂量或时间）与之关联。

为了评估 CPA 在提供协变量信息时能否解耦协变量和扰动信息，我们训练了两个模型：

1. 模型 (i)：训练时只提供扰动标签

2. 模型 (ii)：训练时同时提供扰动和细胞类型标签

然后我们检查了这两个模型得到的基础潜在表示，即在通过对抗训练将协变量和/或扰动信息转移到各自嵌入后剩余下来的潜在信息。正如预期：

• 模型 (i) 的潜在表示显示扰动混合良好，但细胞类型信息被保留（附录图 S2B）。

• 模型 (ii) 的潜在表示显示细胞类型和扰动都混合良好，因为在这种情况下，两者的标签都提供给了模型，模型成功地将这些信息嵌入到了相应的潜在因子中（附录图 S2C）。

为了进一步展示两个模型之间的差异，我们解码了 B 细胞的基础潜在表示，不加入细胞类型和扰动嵌入，并观察了 CD74 和 CD37 的基因表达，这两个基因是该数据集中 B 细胞的标记基因（附录图 S2D）。从模型 (i) 得到的细胞显示这些基因高表达，这是因为细胞类型信息没有被解耦，被保留在了基础潜在空间中；而模型 (ii) 则不是这样，因为细胞类型信息已成功从基础潜在空间中移除（附录图 S2E）。

• binary stimulation：二元刺激，即“有 IFN-β / 无 IFN-β”。

• without continuous covariate：没有连续协变量（如剂量、时间）。

• 这意味着：这个数据集只能测试离散协变量（细胞类型）的解耦，不能测试剂量/时间的处理。

• 模型 (i)：只告诉模型“有没有 IFN-β 刺激”，不告诉它“这是什么细胞类型”。

• 模型 (ii)：同时告诉模型“有没有 IFN-β 刺激”和“这是什么细胞类型”。

• 这是一个对照实验：看提供细胞类型标签后，CPA 是否能把它从基础状态中解耦出去。

• without factoring the cell type and perturbations embeddings：不加扰动/协变量嵌入，只看基础状态本身能解码出什么。如果 里保留了细胞类型信息，那么解码后应该能看到 B 细胞标记基因高表达。

此外，我们利用这个数据集来评估模型是否能够学习细胞类型特异的扰动响应。我们训练了多个 CPA 模型，每次留出某一种细胞类型的扰动细胞，然后预测这些 OOD 条件下扰动细胞的基因表达。CPA 成功建模了不同细胞类型中各种基因对刺激的响应，即使在细胞类型特异的响应情况下也是如此（附录图 S3）。Srivatsan 等人（2020）的 sci-Plex2 数据集包含人肺腺癌细胞系（A549）在四种药物扰动和递增剂量下的测量结果（附录图 S4A）。在这种情况下，模型学习泛化到未见过的药物剂量。为了展示 OOD 特性，我们留出了所有药物中第二高剂量的细胞。之所以这样选择，是因为在最高剂量下，大多数药物导致绝大多数细胞死亡，我们没有足够的细胞来测试 CPA 模型的泛化能力。由于 CPA 学习到的潜在空间表示仍然是高维的，我们可以使用各种降维方法来可视化它。我们选择了使用余弦相似性核计算的 Kernel PCA（附录图 S4C）。为了展示 CPA 对单个基因的剂量-响应动态的捕捉能力，我们观察了所有扰动下排名前 2 的差异表达基因（附录图 S5）。剂量-响应曲线与观察数据吻合良好。我们通过计算预测基因表达与真实基因表达的均值和方差之间的 R² 分数来评估预测的优劣。我们在整个基因表达向量上计算 R²，也在该条件下排名前 50 的差异表达基因（DEGs）上单独计算，以确保我们捕捉到了感兴趣的基因的响应（附录图 S4B）。为了有一个参照点，我们构建了一个基线，即计算 OOD 条件与训练数据集中随机子集之间的 R² 分数（见材料与方法中的基准部分）。优于这个基线表明模型确实学到了扰动和协变量信息，而不仅仅是建模了训练数据的平均表示。

• cell-type-specific responses：细胞类型特异的响应，即同一种扰动（IFN-β）在不同细胞类型（B 细胞、T 细胞、单核细胞）中引起的基因表达变化不同。

• 具体做法：

  • 模型 1：留出 B 细胞，用其他细胞类型训练，预测 B 细胞的扰动响应。

  • 模型 2：留出 T 细胞，用其他细胞类型训练，预测 T 细胞的扰动响应。

  • 以此类推。

• generalize to unseen dosages：泛化到未见过的剂量。例如：训练时只有 10nM、1μM、10μM，预测 100nM。

• 高维空间无法直接可视化，需要降到 2D。选择 Kernel PCA 是因为潜在空间是非线性的，普通 PCA 不够用。

• dose-response dynamics：剂量-响应动态，即基因表达随剂量变化的曲线。

• R² score：决定系数，衡量预测值与真实值的吻合程度，越接近 1 越好。

• top 50 DEGs：排名前 50 的差异表达基因。为什么要单独看 DEGs？因为全部基因中大部分没变化，R² 可能被“稀释”；看 DEGs 才能确认关键基因预测准不准。

我们额外使用未见条件的嵌入与观察流形中最近嵌入之间的距离，作为不确定性的代理指标（见材料与方法）。对于训练时观察到的条件，这个距离等于 0；对于未呈现给模型的扰动和协变量空间中的点，这个距离增大，例如训练数据集中采样剂量之间的剂量值（附录图 S4D）。对于药物组合，这个距离增大（附录图 S4E），这一发现与以下事实一致：模型在训练期间从未见过某些药物组合，且此类预测是在与训练时观察到的条件更远的条件下进行的。

作为第三个例子，我们研究了 Hagai 等人（2018）的跨物种数据集。在这里我们展示，协变量的动态可以是非单调函数，例如时间而非剂量-响应。在这个例子中，来自小鼠、大鼠、兔子、猪的骨髓来源单核吞噬细胞被 LPS 扰动（附录图 S6A）。CPA 能够建模感兴趣的基因随时间的响应（如原始研究中所指出的）（附录图 S6B）。

• distance equals 0：距离为 0，说明查询条件在训练集中出现过。

• increases for points not presented：未出现过的点，距离增大。

• dosage values between those sampled：训练时采样剂量之间的剂量值（如训练有 10nM 和 1μM，查询 100nM）。

• time instead of dose response：时间是协变量，而不是剂量。

• 剂量-响应通常是单调的（剂量越高，效应越强）。

• 时间响应通常是非单调的（如炎症因子先升后降）。

CPA在大规模单细胞高通量筛选中发现可解释的潜在空间

最近提出的 sci-Plex 检测（Srivatsan et al, 2020）通过核哈希（nuclear hashing）在单个实验中刻画数千种独立扰动。利用这个高通量筛选，188 种化合物在三种癌细胞系中进行了测试。这个化合物组合被选中，以靶向多样化的靶点和分子通路，涵盖转录和表观遗传调控因子以及多种作用机制。被筛选的细胞系 A549（肺腺癌）、K562（慢性髓性白血病）、MCF7（乳腺腺癌） 分别暴露于这 188 种化合物的四种剂量（10 nM、100 nM、1 μM、10 μM），共生成了 29 万个细胞的 scRNA-seq 图谱（图 2A）。与之前一样，我们将数据集划分为三个子集：train、test、OOD。对于 OOD 情况，我们留出了在三种细胞系中效应最强的 36 种药物的最高剂量（10 μM）。OOD 情况下的药物、剂量和细胞系组合已从训练集和测试集中移除。

• nuclear hashing：核哈希，用寡核苷酸标签给不同样本的细胞核“贴条形码”，混在一起测序，降低成本。

• OOD 集 = 36 药物 × 3 细胞系 × 最高剂量。

• sci-Plex3 数据集包含 188 种化合物 × 3 种癌细胞系 × 4 种剂量，共约 29 万细胞；OOD 集留出了 36 种效应最强药物的最高剂量（10μM），这些药物-剂量-细胞系组合完全从训练和测试集中移除，用于严格测试 CPA 对未见条件的泛化能力。

CPA 能以出人意料的准确度外推到未见过的 OOD 条件，因为它捕捉到了某个化合物在对照和处理条件之间的差异——即使它从未见过该化合物的最高剂量示例。举一个例子，Momelotinib 与对照相比有强烈的差异响应，可以从排名前 5 的差异表达基因的分布中看出（图 2B）。尽管在三种细胞系中都没有见过 Momelotinib 在最高剂量下的效应，CPA 仍能合理地推断这些基因的均值和分布（图 2B）。CPA 在建模未见过的扰动方面表现良好，因为在 OOD 条件下，真实值与预测值的相关性总体上优于不同化合物之间目标细胞与已有细胞的相关性（图 2D）。观察单个条件时（图 2C），CPA 能很好地重现 OOD 条件下低表达和高表达基因。此外，我们在这个数据集上将 CPA 的性能与 scGen 进行了比较。scGen 是另一种用于单细胞数据集扰动预测的深度学习方法。由于该模型无法处理与扰动相关的连续协变量，我们在基准测试中只保留对照和第二高剂量的细胞。我们比较了预测均值和方差的 R²，以及逐基因计算的 Wasserstein 距离，以衡量整体分布。所有指标都在整个基因向量和DEGs上计算。将 CPA 的预测与 scGen 比较时，我们观察到在更直接的均值预测上有 1.54% 的提升，而高阶矩上的提升更大。同时，CPA 在方差预测上显著优于 scGen，提升了 35.85%，在整体分布预测上也同样优于 scGen。这个基准测试是一个非常简化的场景，现有模型可以在此基准上测试，它们在均值预测上表现良好，但无法捕捉整体分布的变化——而分布变化在单细胞数据中至关重要，因为它捕捉了细胞异质性（见附录图 S7）。

• 因为 scGen 不能处理连续剂量，只能把剂量当离散类别，所以只选两个剂量。

• higher moments：高阶矩，即方差、偏度、峰度等。

当预测实验中包含更多未见协变量时，CPA 的表现更差。为了评估模型泛化到未见条件的能力，我们在 28 种不同的划分上训练 CPA，每种划分留出不同条件，其中某个剂量在 1–3 个细胞系中被留出（图 2E）。我们在这里看到，K562 是最难泛化的细胞系——当考虑用两个细胞系训练、泛化到另一个细胞系时。我们还发现，外推到最高剂量比插值中间剂量更难，这与以下事实一致：预测更高剂量的实验效应比用 S 形曲线拟合中间剂量更困难。

训练后，CPA 学习到 188 种化合物的压缩表示，每种药物由一个 256 维向量表示（图 2G）。为了测试学习到的药物嵌入是否有意义，我们考察作用机制相似的化合物在潜在空间中是否相似。这在大量主要作用机制上都成立：我们发现表观遗传、酪氨酸激酶信号、蛋白质合成类化合物被模型聚在一起，这表明这些机制的药物在这三种癌细胞系上有效，与原始发表（Srivatsan et al, 2020）的发现一致。

我们进一步证明，模型学习到了化合物之间的普适关系，这些关系在不同数据集和模态上都成立。使用 sci-Plex 数据集中测试的同一组化合物加上 853 种其他化合物（共 1,000 种化合物），我们在 L1000 bulk 扰动测量数据上训练 CPA，涵盖 82 种细胞系（Musa et al, 2019）。我们观察到 CPA 在 bulk RNA-seq 数据上同样表现良好，而且 sci-Plex 中出现的表观遗传和酪氨酸激酶信号化合物在潜在表示中彼此接近，这表明学习到的模型相似性跨数据集适用（图 2H）。这对其他学习到的嵌入也成立：对协变量嵌入（这里是 82 种细胞系）应用相同的相似性度量，我们观察到模型学习到的细胞系嵌入也反映了细胞系在扰动响应上的相似性——来自血液组织的细胞系聚在一起（图 2F）。

• extrapolating to the highest dose：外推到最高剂量（如训练有 10nM、100nM、1μM，预测 10μM）。

• interpolating intermediate doses：插值中间剂量（如训练有 10nM、1μM，预测 100nM）。

• sigmoidal behavior：S 形曲线，剂量-响应通常呈 S 形。

图2：从188种药物在癌细胞系中的大规模单细胞筛选数据中学习药物和细胞系的潜伏表示

图2：

• A：sci-Plex 样本（n=290,889）的 UMAP 表示，包含 A549、K562、MCF7 三种细胞系，按化合物靶向的通路着色。
• B：A549 细胞在 Momelotinib（一种 JAK 抑制剂）最高剂量处理下，排名前 5 的差异表达基因在真实、对照和 CPA 预测细胞中的分布。
• C：5,000 个基因和排名前 50 的 DEGs在 CPA 预测细胞和真实细胞之间的平均基因表达，其中排名前 5 的 DEGs 用红色高亮，针对四种模型未见过最高剂量的化合物。
• D：R² 分数箱线图，针对 36 种化合物和 108 个独特留出条件（跨不同细胞系）的预测和真实细胞。基线表示将每个条件与从随机采样细胞导出的均值进行比较。
• E：所有基因和排名前 50 DEGs 的 R² 分数箱线图。每一列代表一个场景：特定剂量的所有化合物在某个细胞系或细胞系组合上被留出训练，之后被预测。
• F：CPA 从 L1000 数据集学到的 82 种细胞系潜在表示，部分癌细胞系按组织来源着色。
• G：CPA 学到的潜在药物嵌入的二维表示。与表观遗传调控、酪氨酸激酶信号、蛋白质折叠/降解通路相关的化合物按其已知通路着色。右上角小图显示数据中三种细胞系的潜在协变量嵌入，表明没有特定的相似性偏好。
• H：在 L1000 数据集中测试最多的前 1,000 种化合物的 bulk RNA 细胞系谱上训练的 CPA 模型的潜在药物嵌入。与 (A) 中 sci-Plex 实验重叠的化合物按与 (G) 相同的通路标签着色。

(E 和 F) 中的箱线图表示中位数（中心线）和四分位距（铰链），须线表示最小值和最大值。

我们通过在 A549 细胞中使用原始 sci-Plex 中的 13 种化合物进行新的组合实验，进一步验证了在 sci-Plex 数据上训练的模型。我们利用 CPA 在 sci-Plex3 数据（图 2）上训练时预测的扰动响应，挑选了高响应的扰动（这也反映在模型的扰动嵌入中——它们形成独立于对照的聚类），用来进行额外的验证实验。我们选择的组合覆盖多种通路和响应强度，使用原始实验中的第二高剂量，以捕捉最大的细胞变异性。我们看到这些组合分为两个行为聚类（图 3A），其中较小的聚类主要由 Alvespimycin 的转录响应主导。然后我们评估了 CPA 模型预测留出的扰动组合的能力（图 3B）。

CPA 成功预测了以下情况的转录响应（图 3B）：

1. 与对照相似的化合物

2. 由一种化合物主导的组合

3. 包含两种化合物转录响应的组合

在图 3B 中可以看到，Panobinostat + Alvespimycin 组合（图 3A 左上聚类）被预测的 R² 为 0.81，尽管模型没有见过任何与之相似的细胞——因为该聚类只包含 Panobinostat + Alvespimycin 这一个条件，因此完全被留出。CPA 表现优于之前描述的对照随机基线模型和一个线性模型（见材料与方法中的基准部分），并准确预测了高变基因的表达水平（图 3C–E）。然后我们可以通过查看从 CPA 扰动潜在空间组合单个扰动向量得到的潜在空间，来重建组合之间的表示（图 3F）。我们看到三个“效应聚类”，分别由 Givinostat、Panobinostat 和 Alvespimycin 代表（图 3F）。Givinostat 和 Panobinostat 都是组蛋白去乙酰化酶（HDAC）抑制剂，但通过不同的作用机制起作用。

图3：在新的大规模药物组合数据集上对预测结果进行验证

图3：

• A：combosciplex 数据集的 UMAP 表示，包含 n=63,378 个细胞和 32 种扰动及扰动组合，在 A549 细胞中。左侧 UMAP 突出显示了后续结果所用的划分。右侧 UMAP 显示了选出的五个分布外（OOD）条件及其在聚类中的表达模式。
• B：(A) 中显示的五个分布外条件及每个条件的模型性能。条件名称左侧的圆圈表示对预测的定性难度评估。绿色标签的条件由两种化合物中的一种主导，模型应该相对容易预测；而包含 Alvespimycin 的组合在转录上与训练中见过的条件更不相似。
• C：CPA vs. 线性模型 vs. 前述随机基线的基准测试，用 R² 衡量，分别在高变基因和排名前 50 的差异表达基因上计算。箱线图表示中位数（中心线）和四分位距（铰链），须线表示最小值和最大值。
• D：预测 vs. 真实的处理后表达值，前 5 个 DEGs 用红色标出。
• E：每个分布外条件中排名前 2 的 DEGs的小提琴图，以及刺激前、刺激后和 CPA 预测的表达值。
• F：CPA 学到的组合潜在向量的 UMAP 表示。四个单独条件及它们出现的组合被高亮显示。

CPA 是一种用于预测单细胞扰动的可扩展框架

CPA 采用的模块化架构的好处之一是它的灵活性和可扩展性。虽然扰动字典对潜在空间中的组合结构很有效，但它仅限于训练集中出现的化合物。因此，预测实验中未筛选过的化合物的扰动响应是不可行的。为了实现对未见药物的预测，Hetzel 等人（2022）提出了一个最近的扩展，称为 chemCPA（见材料与方法）。chemCPA 引入了一个扰动网络，使用已知的化学描述符来编码小分子（图 4A）。这个扰动网络取代了 CPA 中的扰动嵌入字典。它包含一个预训练的分子编码器 G、一个扰动编码器 M，以及一个摊销的剂量缩放器 S。M 学习从通用化学嵌入映射到化合物 d 的潜在扰动效应 z_di，而剂量缩放器 S 学习在给定化学嵌入 h_di 和施加剂量 s_i 的情况下映射到有效剂量 ŝ_i。

• perturbation dictionary：扰动字典，即存储每个扰动嵌入向量的“字典”。

• compositional structure：组合结构，即潜在空间中的加性组合。

• M：化学嵌入 → 潜在扰动效应。

• S：化学嵌入 + 施加剂量 → 有效剂量。

• 这样，即使用户没见过的药物，只要它有化学结构，就能生成嵌入。

我们将 chemCPA 应用于 sci-Plex3 数据集（图 2）和新的组合数据集（图 3）。在单药预测实验中，我们留出了 9 种化合物作为 OOD（与 Hetzel et al, 2022 类似）——Dacinostat、Givinostat、Belinostat、Hesperadin、Quisinostat、Alvespimycin、Tanespimycin、TAK-901 和 Flavopiridol。由于 CPA 的扰动字典仅限于训练集中观察到的化合物，当这些药物完全从训练中排除时，无法比较 CPA 和 chemCPA。我们在训练和验证集中保留了最低两个剂量的观测，以便在具有挑战性的设置中进行比较。因此，排除了 1μM 和 10μM 的剂量值。这种设置下的结果说明，化学先验在整个基因集和DEGs上、跨剂量都改善了扰动预测（图 4B）。此外，我们观察到 chemCPA 对属于组蛋白去乙酰化通路的化合物泛化特别好（图 4C），这与原始 sci-Plex 发表（Srivatsan et al, 2020）以及该扰动数据集中的普遍支持一致。我们报告了 Hetzel 等人（2022）的结果，其中 OOD 药物完全从训练和验证中排除（图 4D），并观察到包含低剂量值反而显著降低了 chemCPA 的性能，这表明高剂量观测对 CPA 和 chemCPA 的训练至关重要。

• CPA：如果药物完全没见过，字典里没有它的嵌入，无法预测。

• chemCPA：即使药物完全没见过，只要有化学结构，就能生成嵌入。

• 妥协方案：让 CPA 至少见过这些药物的低剂量数据，这样才能对比；但高剂量（1μM、10μM）仍然留出。

• chemical prior：化学先验，即分子结构信息。

我们对药物组合数据集进行了类似的分析。在这里，OOD 集由所有包含 Panobinostat 的 9 种药物组合组成，分别与 SRT3025、PCI-34051、Sorafenib、Dasatinib、SRT2104、SRT1720、Crizotinib、Alvespimycin 和 Curcumin 之一组合。也就是说，Panobinostat 的单药观测保留在训练和验证集中，这使得 CPA 与 chemCPA 可以对比。虽然两个模型在这个场景下都表现良好，但 chemCPA 整体上优于 CPA。我们应该强调，Alvespimycin 和 Panobinostat 组合引起的明显分布偏移（图 3A）无法被任一模型完全识别，两个模型在 DEGs 上的 R² 都只有约 0.3。相比之下，其他组合在 chemCPA 模型下 DEGs 的 R² 高于 0.6（图 4E）。对 CPA 提升最显著的是 Panobinostat 和 SRT2104 的组合（图 4F）。通过 chemCPA 的扰动网络，可以将单药扰动从训练中排除。这个 OOD 设置的结果（图 4G）表明，在组合场景中，获取单药影响是至关重要的。虽然 DEGs 的中位 R² 从 0.85 下降到 0.38，但 Panobinostat 的单药扰动仍然可以被高精度预测（图 4H）。

• chemCPA 的优势来自化学先验。

图四：CPA的可扩展性使得能够预测对未见过药物的反应

图四：

• A：为 CPA 引入化学先验知识的架构修改方案。分子编码器 G 将化合物的化学信息映射到潜在药物嵌入 h_di。这个模块可以基于预训练的图编码器或分子指纹（如这里所用的）。训练时只优化扰动模块 M 和摊销剂量缩放器 S。
• B：CPA 和 chemCPA 的对比，包括一个忽略所有药物诱导扰动效应的基线。分数在 sci-Plex3（n=290,889） 数据上计算，测试集由九种化合物组成。整个基因集（左）和 DEGs（右）都显示了。
• C：(B) 中模型在两种组蛋白去乙酰化药物（Givinostat 和 Dacinostat）上的性能对比，跨所有剂量和细胞系（阴影）。
• D：(B) 中模型在最高剂量 10μM下的中位分数，包括 chemCPA 在九种测试化合物完全未见、且从训练和验证中排除的设置下的结果。
• E：CPA 和 chemCPA 在新的 combosciplex（n=63,378）上的对比，包括与 (B) 相同类型的基线。测试集由所有包含 Panobinostat 的九种药物组合组成。
• F：(E) 中模型在 DEGs 上 R² 差异最大的两个条件上的详细性能对比。
• G：(E) 中模型的中位分数，包括 chemCPA 在单药 Panobinostat 观测也被留出的设置下的结果。DEGs 上的中位分数从 0.85 降到 0.38。
• H：展示 chemCPA 在 Panobinostat 被留出时预测其单药效应的能力。这与 Panobinostat 包含在训练集时（E）的分数进行对比。

箱线图表示中位数（中心线）和四分位距（铰链），须线表示最小值和最大值。

CPA 允许建模组合基因扰动

组合药物疗法被认为可以解决单一疗法的疗效有限问题（Menden et al, 2019），并防止癌症治疗中的耐药性（Jia et al, 2009; Menden et al, 2019; Adam et al, 2020）。然而，少数基因的组合表达往往驱动细胞层面的复杂性，导致新特性、新行为和多样化细胞类型的出现（Norman et al, 2019）。为了研究这种基因相互作用（GIs），最近的扰动 scRNA-seq 检测允许我们测量细胞对单个或组合基因扰动的基因表达响应（Dixit et al, 2016; Datlinger et al, 2017）。虽然实验方法对于评估组合疗法的效果是必要的，但在实践中，没有计算预测，实验上探索所有可能的组合是不可行的。

为了实现这一目标，我们将 CPA 模型应用于从 Perturb-seq（单细胞 RNA 测序混合 CRISPR 筛选）收集的 scRNA-seq 数据，以评估 105 个基因的单扰动或组合相互作用（即单基因 x、单基因 y、以及 x+y 对）如何影响 K562 细胞的生长（Norman et al, 2019）。总的来说，这个数据集包含 284 个条件，测量了约 108,000 个单细胞，其中 131 个是独特的组合对（即 x+y），其余是单基因扰动或对照细胞。我们观察到，潜在基因相互作用流形将诱导已知且相似基因程序的 GIs 放置在彼此接近的位置（图 5A）。我们进一步将我们的潜在空间聚类与平均基因表达嵌入进行了比较，类似于原始发表（附录图 S8）。总的来说，CPA 潜在空间与平均基因表达嵌入一样能分组相似扰动，但实现了更细粒度的聚类，从而获得更好的聚类指标。接下来，我们试图评估预测特定基因相互作用的能力。我们研究了 CBL 和 CNN1 之间的协同相互作用，该作用驱动红细胞分化，此前已被验证（Norman et al, 2019）。我们训练了一个 CPA 模型，将 CBL+CNN1 从训练数据中留出。任一基因的过表达都会导致细胞从对照向单扰动和双扰动细胞（附录图 S9A）向红细胞基因程序进展。CBL 和 CNN1 双过表达上调了已知基因标记（Norman et al, 2019），如血红蛋白（见图 5B 中的 HBA1/2 和 HBG1/2）。我们观察到，我们的模型成功预测了这种协同相互作用，重现了与真实数据相似的模式，并与原始发现一致（图 5C）。我们进一步评估 CPA 预测此前报道的（Norman et al, 2019）DUSP9 和 ETS1 之间的遗传上位相互作用，导致双扰动细胞中 DUSP9 表型占主导（图 5C）。

• Perturb-seq：单细胞 RNA 测序 + 混合 CRISPR 筛选。

• overexpression：过表达（不是敲除）。

为了系统评估 CPA 的泛化行为，我们训练了 13 个不同的模型，每次留出约 10 个独特组合的所有细胞，覆盖数据集中所有 131 个双扰动条件，这些条件在训练后被预测。我们将 CPA 的预测与上一节描述的线性基线和随机基线进行了比较（图 5D）。令人惊讶的是，基线方法取得了与 CPA 和线性模型相当的准确预测。这一观察表明，不同扰动引起的转录组效应非常相似，且该数据集在全转录组水平上包含有限的非线性效应。因此，线性模型或简单基线可以在均值上提供准确预测，这与之前展示的组合药物效应不同。报告的 R² 值显示，大多数扰动的预测是稳健的；较低分数出现在以下情况：由于样本稀缺（n<100）导致评估有噪声，或当其中一个扰动在数据中只有单扰动细胞时，导致模型无法预测未见过的组合（图 5E，见附录图 S9 和 S10）。为了进一步理解 CPA 性能何时下降，我们首先在训练时没有见过任何组合的子集上训练它，然后逐渐增加训练中见过的组合数量。我们发现，当模型用更多组合训练时，整体预测准确度提高；而当用较少组合训练时，它可能无法预测 DEGs（见图 5F 中 n < 71 组合）。

因此，一旦用足够大且多样化的训练数据训练，CPA 就能稳健地预测未见过的扰动。我们接下来问，我们的模型能否泛化到测量过的组合之外，并在计算机中生成所有 5,329 个组合——这些组合在真实数据集中未被测量，但占所有可能性的约 98%。为了研究这些预测的质量，我们训练了一个所有组合都在训练中出现过的模型，以实现最大的训练数据和样本多样性。然后我们为所有缺失的组合预测了 50 个单细胞。我们对所有测量和生成的数据的平均基因表达向量进行了嵌入（图 5G），同时为每个条件报告了不确定性值（图 5H）。我们假设，生成的嵌入越接近测量数据，它就越可能探索测量数据周围相似的遗传流形空间。有了这些信息，我们基于基因流行度对 MIS 聚类进行了注释，发现单基因（即基因 x）与其他基因（即 y）配对形成的组合（即 x+y）是聚类分离的主要驱动因素（图 5I）。没有测量过双扰动的基因，用新预测的转录组表达不太可能分离为独立聚类（附录图 S11A），这表明它们的互作特异性效应比训练数据中至少有一个双扰动的基因更少变化。同时，远离测量数据的点可能指示新的基因互作行为，例如 KLF1 与共扰动基因的情况（图 5J）。然而，这需要额外的考虑和验证步骤。

• gene prevalence：基因流行度，即某个基因出现在多少组合中。

图五：CPA 预测组合基因扰动

图五：

• A：使用 CPA 对来自 Perturb-seq（n=108,497） 的 281 个单基因和双基因扰动推断的 UMAP 潜在空间。每个点代表一个基因扰动。着色表示来自原始发表的与扰动基因相关的已知基因程序。
• B：与协同基因对（CBL + CNN1）相关的细胞的测量和 CPA 预测基因表达。基因名称取自原始发表。
• C：与 (B) 相同，但针对上位基因对（DUSP9 + ETS）。显示了 DUSP9 + ETS 共扰动细胞 vs 对照细胞的前 10 个 DEGs。
• D：测量和预测细胞的平均基因表达 R² 值，针对所有基因（黑色）或前 100 个 DEGs，预测所有 131 个组合，对比 CPA、线性模型和基线（13 个训练模型，每次约 10 个测试组合）。
• E：预测和真实平均基因表达的 R² 值 vs 真实数据中的细胞数量。
• F：预测和真实细胞的 R² 值 vs 训练中见过的组合数量。
• G：每个测量（n=284，红点）和 CPA 预测（n=5,329，灰点）扰动组合的平均基因表达 UMAP。
• H：与 (G) 相同，但显示测量不确定性（余弦相似度）。
• I：与 (G) 相同，但显示 Leiden 聚类中的主导基因（25 个或更多观测）。KLF1 聚类用红色虚线圆圈高亮。
• J：以 KLF1 作为常见共扰动基因的细胞平均表达 UMAP。
• K：KLF1 与共扰动基因之间线性回归相关指标的层次聚类，在测量和预测细胞中。每一行表示用单扰动作为预测变量拟合双扰动线性回归时得到的汇总参数，以及系数 c1 和 c2 之间的关系。定义见材料与方法。
• L：RF 选择的基因（x 轴）在测量（紫色）和预测（黄色）扰动（y 轴）中的缩放基因表达变化（vs 对照）。标题表示逐基因回归系数，以及 Norman 等人（2019）建议的互作模式。

箱线图表示中位数（中心线）和四分位距（铰链），须线表示最小值和最大值。

为了研究新预测条件之间的相互作用类型，我们比较了双扰动和单扰动相对于对照细胞的差异，并据此注释了它们的相互作用模式，改编自原始发表（Norman et al, 2019）。在每个基因特异的聚类中，我们观察到这些值存在变异性，表明我们的预测包含了超出单基因扰动效应的细粒度信息，且无法用二维嵌入完全解析。在使用这些指标和可用实验数据水平（附录图 S11B）对基因扰动进行整理后，我们决定当至少一个基因有双扰动测量时，基于这些值预测和注释互作模式。例如，我们观察到 KLF1 及其伙伴基因扰动对仅从这些指标上就在大多数测量数据点中形成聚类，表明存在多种互作模式，这些模式无法用二维嵌入完全描述（图 5K）。基于这些值及其变异性提取基因对，使我们能够提取互作类型，并通过阈值化将其标记为可解释的类别，如 Norman 等人之前提出的。当我们进一步检查每个共扰动中的差异表达基因时，这些指标验证了此前报道的上位相互作用（CEBPA），并提出了新的互作：KLF1 主导行为（NCL）、基因协同（FOXA3）、上位（PTPN13）等（图 5L）。对所有测量和预测的双扰动重复此分析，我们发现了具有潜在互作流行度的基因（附录图 S11C）。在反复响应多种扰动的基因中，我们发现了方向和幅度上的共同基因表达趋势（附录图 S11D），表明变异受保守的基因调控原则调节，这些原则可能被我们学习到的模型捕捉到了。总之，我们的分析表明，CPA 可以通过利用基因扰动数据生成双扰动测量，当与不确定性指标结合时，使我们能够在预测的基因-基因扰动中生成和解释基因调控规则。

讨论

在计算机中预测细胞对扰动的响应行为，对于最优实验设计和识别有效药物和治疗至关重要。通过 CPA，我们引入了一种通用且可解释的方法，用于在单细胞分辨率下建模细胞行为。CPA 实现为一个使用随机梯度下降训练的神经网络，可扩展到数百万个细胞和数千个基因。

我们将 CPA 应用于各种数据集和任务，从预测单细胞响应到学习嵌入，以及重建化合物的表达响应，涵盖可变的药物-剂量组合。具体来说，我们展示了跨剂量水平和时间序列的扰动建模，并展示了在药物扰动研究以及多基因敲除的基因扰动检测中的应用，揭示了由我们模型预测值推断的潜在基因-基因相互作用模式。CPA 结合了线性分解模型的可解释性与非线性嵌入模型的灵活性。

虽然 CPA 在我们的实验中表现良好，但众所周知，机器学习中没有一种方法适用于所有情况，与任何其他机器学习模型一样，如果测试数据与训练数据差异很大，CPA 就会失败。为了提醒 CPA 用户注意这些情况，量化模型不确定性至关重要。为此，我们实现了一个基于距离的不确定性分数来评估我们的预测。当前的不确定性估计启发式源自 CPA 的组合公式。这种公式导致与训练数据中观察到的协变量/扰动向量组合差异很大的组合会有更高的距离，因此比训练中观察到的协变量/扰动具有更高的不确定性（见附录图 S4 和 S12）。然而，与其他计算机预测模型一样，CPA 的预测不能取代实验验证；相反，它可以作为高效进行实验的潜在指南。此外，可扩展的贝叶斯不确定性模型是未来工作的有前景的替代方案（Gal & Ghahramani, 2016）。虽然我们选择实现确定性自编码器方案，但向变分模型（Lopez et al, 2018; Lotfollahi et al, 2020）以及均方误差以外的损失函数（Eraslan et al, 2019）的扩展是直接的。除了模型不确定性，数据不确定性和偏差在将模型偏向特定预测区域方面也起着重要作用。例如，与其他基因共同测量的基因往往比在组合中只有一个或两个基因都未被共同测量的基因显示出更复杂的预测；在后一种情况下，模型倾向于预测单基因效应（显性）（见附录图 S11D）。这些多变量偏差的本质，以及它是纯粹数据驱动的还是与生物学相关的，应该成为额外工作的焦点。我们希望这一见解有助于未来组合筛选的设计。最后，我们使用了一个固定架构，超参数非常相似。然而，为了达到最佳性能，需要系统性的超参数扫描；例如，在存在许多扰动和协变量的情况下，可能需要不同的解耦损失权重，以确保基础状态不含扰动和协变量信息。

• 如果训练数据有该基因的组合数据 → 预测更复杂、更准确。

• 如果训练数据没有该基因的组合数据 → 模型倾向于预测“单基因主导”。

除了 CPA，现有方法（Lotfollahi et al, 2020; Russkikh et al, 2020），如 scGen（Lotfollahi et al, 2019），也已被证明能够在数据集不包含组合处理或剂量依赖扰动时，预测单细胞扰动响应。因此，在扰动较少的简单场景下，将 CPA 与这些方法进行基准对比可能是有益的。然而，考虑到当前生成大规模扰动研究的趋势（Dixit et al, 2016; Norman et al, 2019; Srivatsan et al, 2020），这种方法可能不切实际。

目前，该模型仅基于基因表达，因此无法直接捕捉其他层次的相互作用或效应，例如转录后修饰、信号传导或细胞间通讯。然而，由于基于神经网络方法的灵活性，CPA 可以扩展到包括其他模态，例如通过多模态单细胞 CRISPR（Frangieh et al, 2021; Papalexi et al, 2021）结合 scRNA-seq 和 ATAC-seq（Clark et al, 2018; Chen et al, 2019）以及 CUT&Tag（Kaya-Okur et al, 2019; Wu et al, 2021）。特别地，我们期望空间转录组学（Rodriques et al, 2019; van den Brink et al, 2020）成为 CPA 扩展的宝贵来源，因为它样本数量高，且深度学习模型在计算机视觉中占主导地位。

• 逻辑：

  • 空间转录组学数据量大（像图像）。

  • 深度学习在图像领域非常成熟。

  • 所以空间转录组学是 CPA 扩展的天然方向。

CPA 模型不限于单细胞扰动。虽然我们选择单细胞设置是因为可用样本数量高，但 CPA 可以轻松应用于大规模 bulk 队列，其中的协变量可能是患者 ID 或转录因子扰动。这些以及任何其他可用的属性都可以独立控制（Lample et al, 2017），以实现组合性、可解释的预测。任何 bulk 组合模型都可以与较小规模的单细胞模型结合，以构成真正的多尺度观测方差模型。DL 设置的灵活性也将允许在扰动或协变量潜在空间上添加约束。例如，这些约束可以是化合物的化学相似性（Mater & Coote, 2019），或临床协变量引起的患者 ID 差异。CPA 与普通自编码器的关键区别在于其潜在空间解耦，以及由此产生的在细胞状态和协变量背景下扰动的可解释性。最终，生物学的任何目标都不仅是盲目预测，而是机制理解。这一目标体现在 DL 模型在序列基因组学中的方向——不仅预测新相互作用，还解释学习到的基因调控密码。因此，我们相信 CPA 不仅可以作为计算机筛选的假设生成工具，还可以作为整体数据近似模型。偏离我们假设的数据生成过程（见材料与方法）将告诉我们给定数据集中缺失的信息和/或因子模型中缺失的方面。通过包含多层调控，所得模型可以在预测和机制理解（例如协同基因调控或其他相互作用）方面增长灵活性。

bulk RNA-seq 是把一群细胞（通常是成千上万甚至上百万个细胞）混合在一起测序，得到的是平均基因表达。

最后，我们期望 CPA 能促进基于表达的扰动筛选的新机会，不仅学习最优药物组合，还能根据个体细胞响应定制实验和治疗，实现个性化。

材料与方法

数据生成过程

我们考虑一个数据集 ，其中每个  描述细胞 i 的 G 个基因的基因表达。扰动向量  包含元素 ，描述施加到细胞 i 上的药物 j 的剂量。如果 d_{i,j} = 0，这意味着扰动 j 没有施加到细胞 i。除非另有说明，后续假设为列向量。类似地，向量之向量  包含额外的离散协变量，如细胞类型或物种，其中每个协变量本身是一个向量。更具体地说，c_{i,j} 是一个 K_j 维的 one-hot 向量。

• D：数据集。

• N：细胞数量。

• x_i：第 i 个细胞的基因表达向量。

• R^G：G 维实数空间。

• G genes：G 个基因（如 5000 个高变基因）。

作用：定义数据集的基本结构——每个细胞一个表达向量。

• perturbation vector：扰动向量。

• M：扰动数量（如 188 种药物）。

• d_{i,j} ≥ 0：剂量，非负。

• dose of drug j：药物 j 的剂量。

我们假设一个未知的生成模型产生了我们的数据集 D。这个生成过程的三个初始组件是：细胞 i 的潜在（未观察）基础潜在状态 ，以及它的扰动向量 d_i 和协变量向量 c_i。我们假设基础潜在状态独立于扰动向量 d_i 和协变量向量 c_i。接下来，我们形成潜在（也未观察）的扰动潜在状态 z_i：

 （1）

在这个方程中，矩阵  的每一列表示 d_i 中 M 个可能扰动之一的 d 维嵌入。类似地，矩阵  的每一列表示第 j 个离散协变量的 d 维嵌入，表示为 K_j 维 one-hot 向量 c_{i,j}。

函数 f_j : R → R 对扰动向量 d_i 中的每个 d_{i,j} 进行非线性缩放，因此实现了 M 条独立的剂量-响应（或时间-响应）曲线。最后，我们假设生成过程通过一个未知的解码分布 返回观察到的基因表达 x_i。这个过程构建了观察 (x_i, d_i, c_i)，然后被包含在我们的数据集 D 中。

• z_i：扰动后的潜在状态。

• z_i^basal：基础状态。

• V^perturbation：扰动嵌入矩阵。

• f_j(d_{i,j})：剂量-响应曲线。

• V^{cov_j}：协变量嵌入矩阵。

• c_{i,j}：协变量 one-hot 向量。

• 这是 CPA 的核心公式：基础状态 + 扰动效应 + 协变量效应。

【生成过程】
  1. 采样基础潜在状态 z_i^basal（与扰动/协变量独立）
  2. 采样扰动向量 d_i（剂量）
  3. 采样协变量向量 c_i（细胞类型、物种等）
  4. 形成扰动潜在状态：
     z_i = z_i^basal + V^perturbation · f(d_i) + V^cov · c_i
  5. 通过解码分布 p(x_i | z_i) 生成基因表达 x_i
  6. 得到观察 (x_i, d_i, c_i)，加入数据集 D

【CPA 的目标】
  反向推断：
  - 从 x_i 推断 z_i^basal（编码器）
  - 学习 V^perturbation、V^cov、f_j
  - 学习解码分布 p(x_i | z_i)（解码器）

CPA

假设上述生成过程，我们的目标是训练一个机器学习模型 ，使得给定数据集三元组 (x_i, d_i, c_i) 以及目标扰动 d'，估计基因表达 x'_i。x'_i 表示反事实分布——即具有协变量 c_i 的基因表达 x_i，如果被 d' 而非 d_i 扰动，会是什么样子。

• M((x_i, d_i, c_i), d')：模型 M，输入是三元组 + 目标扰动。

• dataset triplet (x_i, d_i, c_i)：数据三元组（表达、扰动、协变量）。

• target perturbation d'：目标扰动，即想要预测的新扰动。

• estimates the gene expression x'_i：估计基因表达。

• 反事实问题：如果处理不同会怎样？例如：细胞 X 在药物 A 下的表达已知，那它在药物 B 下的表达会是什么样？

给定数据集和学习目标，我们现在准备描述我们提出的模型——CPA。在下文中，我们将分别描述如何训练和测试 CPA 模型。

Training（训练）

CPA 模型的训练本质上是对数据三元组 进行自编码。也就是说，训练时 CPA 并不试图回答反事实问题。相反，训练过程包括：

1. 将基因表达 x_i 编码为估计的基础状态 ，且不含任何关于 (d_i, c_i) 的信息；

2. 将 与关于 (d_i, c_i) 的可学习嵌入结合，形成估计的扰动状态 ẑ_i；

3. 将 ẑ_i 解码回观察到的基因表达 x_i。

更具体地说，CPA 模型的前向传播首先将观察到的基因表达 x_i 编码为估计的基础状态：

估计的基础状态被用作辅助分类器（见下一段）的输入，也用于计算估计的扰动状态 ẑ_i：
(2)

与公式 (1) 相反，这个表达式引入了三个额外的可学习组件：

1. 扰动嵌入 V̂^perturbation

2. 协变量嵌入 V̂^cov

3. 可学习的剂量-响应曲线 f̂_1, ..., f̂_M，这里实现为小神经网络，并约束满足 f̂_j(0) = 0。

作为最后一步，解码器 f̂^dec 接受估计的扰动状态 ẑ_i，返回 f̂^dec_μ(ẑ_i) 和 f̂^dec_σ²(ẑ_i)，即反事实基因表达 x'_i 的估计均值和方差。

为了训练 CPA 模型，我们使用三个损失函数。第一个是重建损失函数，采用高斯负对数似然：
(3)
其中 s(σ²) = log(exp(σ² + 10^{-3}) + 1) 对方差施加正性约束并增加数值稳定性。这个损失函数奖励端到端自编码过程，如果它产生了观察到的基因表达 x_i。

第二个，根据我们对数据生成过程的假设，我们感兴趣的是从 ẑ_i^basal 中移除关于 (d_i, c_i) 的信息。为了实现这种信息移除，我们遵循对抗方法（Lample et al, 2017）。特别地，我们设置以下辅助损失函数：

函数 f̂^adv_d, f̂^adv_{c_{i,j}} 是一组神经网络分类器，试图在给定估计的基础状态 ẑ_i^basal 的情况下预测 (d_i, c_i)。

给定这组损失函数，训练过程是一个优化问题，重复以下两个步骤：

1. 采样 ，通过更新 f̂^adv_d 和 f̂^adv_{c_{i,j}} 的参数来最小化 l_i^d + Σ_j l_i^c_{i,j}，对所有 j = 1, ..., K；

2. 采样 (x_i, d_i, c_i) ~ D，通过更新编码器 f̂^enc、解码器 f̂^dec、扰动嵌入 V̂^perturbation、协变量嵌入 V̂^cov_j（对所有 j = 1, ..., K）以及剂量-响应曲线估计器 f̂_1, ..., f̂_M 的参数来最小化 l_i - λ · (l_i^d + Σ_j l_i^c_{i,j})。

Testing（测试/预测）

给定一个观察 (x_i, d_i, c_i) 和一个反事实处理 d'，我们可以使用训练好的 CPA 模型来回答：具有协变量 c_i 的基因表达 x_i，如果被 d' 而非 d_i 扰动，它的反事实分布会是什么样？为此，我们遵循以下流程：

1. 计算估计的基础状态 ；

---

懒得写了

---

Lotfollahi, M., Klimovskaia Susmelj, A., De Donno, C. et al. Predicting cellular responses to complex perturbations in high‐throughput screens. Mol Syst Biol 19, MSB202211517 (2023). https://doi.org/10.15252/msb.202211517