用biolord解开单细胞数据

潜在优化作为解开中的归纳偏差

潜在优化是我们方法的关键组成部分。 典型的表示解缠方法使用编码器将原始数据样本映射为潜在代码。 这通常称为摊销推理。 虽然使用编码器网络将样本映射到代码很方便，但 Gabbay 和 Hoshen¹² 表明这种方法可能会取得低于标准的结果。 原因是在训练开始时，编码器（随机初始化）将所有样本属性映射到每个潜在代码（已知和未知）。 虽然损失函数鼓励解开（去除已知属性），但编码器的随机初始化导致优化从完美的纠缠状态开始。 后来的训练迭代很难完全消除这种纠缠。

相反，随机初始化的潜在代码通常不包含有关已知或未知属性的任何信息。 在训练时，与每个样本相对应的潜在代码在未知属性上变得更加丰富，而解开目标则确保它不会获得有关已知属性的信息。 直观地说，防止获取不需要的信息比丢失现有信息更容易。 总而言之，潜在优化通过为学习过程提供更好的初始化来帮助实现更多解缠结的潜在代码。

值得注意的是，参考文献中获得的结果。 ¹² 直接应用于此处介绍的生物环境。 如上所述，通过潜在优化解决的挑战是关于标记属性的。 由于这些是随样本（图像或单细胞测量）一起提供的标签，因此它们在两种设置中是相同的。 因此，潜在优化使我们能够获得关于目标属性的有意义的潜在代码。

生物领主模型

Biolord是一个深度学习生成框架，由多个共同优化的模块组成。 Biolord 的输入是一个数据集 (D={左[left({x}_{c},{y}_{c}right)right]}_{c=1}^{n})， 在哪里 n 是细胞的数量。 对于每个单元格 C, ({x}_{c}in {{mathbb{R}}}^{M}) 代表 中号 测量的特征（例如，基因表达计数的向量或来自 中号 基因），以及 y_C 是一组大小 K 代表已知的细胞属性，例如细胞类型标签、组织来源或年龄。 正如我们在下面详细阐述的那样， y_C，我们在构造其相应的潜在空间时区分分类属性和有序属性。 根据规定，各 K 集合中的元素 y_C 可能具有不同的维度。 给定输入数据集 D，biolord 管道由两个主要组件组成，同时定义和训练（构建细节在以下小节和补充说明中提供） 1）：

1. 1.

   分解的潜在空间——对于每个已知属性，构建一个专用的子网络。 每个子网络的架构是根据属性的类型（分类或有序）选择的，并且用户可以修改其他超参数。 我们表示 z_y 作为每个子网络的输出，它是对应于属性（分类或有序）的潜在空间 y_C， 和 z_你 作为未知属性的潜在空间（图 1）。 1b）。

2. 2.

   生成模块——生成器 G 将连接的分解潜在空间作为输入并输出对测量特征的预测。

值得注意的是，上述优化、分解的潜在空间和生成预测是联合完成的，使得分解的潜在空间中的嵌入相对于生成器的重建误差进行了优化。

已知属性潜在空间

给定已知的属性集， y_C，为每个子网构建一个专用子网 K 属性来表示其相应的潜在空间。 在这里，我们区分了分类属性（其中相似的单元共享类标签）和有序属性（其中属性特征之间的距离编码相似性）。 在有序属性的定义中，我们考虑连续变量和分类序数变量，因为重要的方面是属性的特征包含结构信息。 此外，测量的分类序数变量（例如年龄）通常代表连续变量的样本。 这样，我们构建不同的子网络如下：

• 分类属性子网 – 这些是使用嵌入模块定义的，以便潜在代码， z_y，在属于同一标签的所有单元格之间共享。 嵌入是直接优化的，即通过完整模型的目标函数应用潜在优化。

• 有序属性子网络——为了使用每个有序属性的结构，我们使用编码器； 多层感知器 (MLP)，默认值为深度 = 2，宽度 = 256。MLP 将输入特征映射到相应的潜在空间， z_y，使用完整模型的目标函数进行优化。

未知属性潜在嵌入

我们通过直接优化每个样本的嵌入来学习未知属性的表示。 我们使用正则化嵌入子网络，即高斯噪声的嵌入模块， 这，一个随机变量 (eta sim {mathscr{N}}(0,{sigma }^{2}I,))，具有固定的方差值 p，添加（补充说明 1）。 输出是一个独特的潜在代码， z_你，对于每个细胞，独立于基因表达或已知属性，在训练期间使用潜在优化进行优化。

因此，为每个单元优化唯一的代码可能会阻碍我们解开纠缠的努力； 该模型可以使用未知属性的潜在代码对整个表达信息进行编码，并忽略特定于属性的编码。 关注 Gabbay 和 Hoshen¹²，为了确保已知属性信息不会泄漏到未知属性的表示中，我们将其正则化为两种方式。 首先，我们向嵌入引入加性高斯噪声，其次，我们向损失添加激活惩罚项，限制嵌入的幅度，从而引入最小性损失项，

在哪里 我 是此项的超参数权重。 这些共同强化了未知属性和已知属性的表示之间共享信息的最小化。 也就是说，未知属性的表示被优化以最小化其包含的关于已知属性的信息。

发电机模块

发电机 G 被构造为解码器网络，参数化为 我，它将级联分解的潜在空间作为输入，并输出测量特征的表达分布的参数化（由均值和方差给出），

根据作为模型输入提供的数据、预处理的对数归一化数据、原始计数或峰值、分布、 磷，可以分别遵循高斯分布、零膨胀负二项式或泊松分布^14,26。 为了定义重建和完整性损失项，我们对每个分布使用各自的负对数似然损失， ({rm{NLL}}left(x|{G}_{theta }right))。 遵循 Gabbay 和 Hoshen 提出的原始模型¹²，我们包括一个关于预测平均值的均方误差项， 米_我，由提供， G_我, ({rm{MSE}}left(x,{mu }_{theta }right))调整为 t (‘reconstruction_loss’) 超参数。 这使我们能够直接优化参数分布建模的所有选择的均值预测（补充说明 1）。 因此我们可以将完整性项写为：

模型优化

结合以上内容，我们可以将完整的模型目标写为两项的组合。 第一项通过优化生成器的准确性来引入完整性，第二项强制已知和未知属性的表示之间共享信息的最少性，

由于上面定义的不同组件是联合优化的，因此分解的潜在空间内的嵌入以及生成器的预测输出受到输入测量以及已知属性标签的影响。

Biolord-classify：用部分标签解开生物表征

为了执行半监督解缠，在这种情况下，我们缺少细胞子集的标签，我们采用参考文献中提出的推导。 ¹³。 除了上述生物领主模型组件（分解的潜在空间和生成模块）之外，我们还包括一个分类器， (Cin {mathcal{C}})，对于每个分类属性和一个回归量， (Rin {mathcal{R}})，对于每个有序属性，它们与之前的组件一起训练。

分类器（回归器）将基因表达作为输入并输出类标签/特征。 对于缺少标签的单元格，分类器（回归器）的输出用于完成分解的潜在表示（扩展数据图 1）。 1）。 为了训练分类器（回归器），我们在现有损失函数中添加了一项，以鼓励对可用标签的样本进行正确预测。 对于分类器，我们使用分类交叉熵损失， (H左(y,C左(x右)右))。 对于回归量，我们使用输出和提供的特征之间的均方误差损失， ({rm{MSE}}left(,y,Rleft(xright)right))。 在所有情况下，仅在提供标签的单元格上评估损失（由集合表示） ({X}^{,S},{Y}^{,S})）。 分类损失由以下公式提供：

在哪里 ({Y}_{C}^{,S}左({Y}_{R}^{,S}右)) 表示与相应分类器关联的标签（特征）集， C （回归量， 右）。 ({{mathcal{L}}}_{{rm{cls}}}) 被添加到生物领主目标中，以便所有组件现在都得到联合训练，

通过将分类模块（分类器和回归器）作为生物领主训练过程的一部分，我们允许在半监督环境中训练生物领主模型，因为分类器和回归器用于估算用作分解潜伏输入的缺失标签细胞。 此外，估算的标签可用于数据的下游分析（图 1）。 1c）。

Biolord 支持的下游分析

Biolord 使用解耦的潜在表示、生成管道和 Biolord 分类模块来实现多种下游分析任务（图 1）。 1c）。 在 Biolord 框架内，我们提供实用函数来实现这种分析。 下游任务如下：

1. 1.

   潜在空间表示——潜在空间嵌入提供了对特定属性内以及不同属性之间的结构的洞察。 潜在表示是一组向量，将已知属性的状态映射到 n_潜 维度状态。 现在可以应用任何下游分析工具来揭示状态之间的属性和关系，例如相关性分析、聚类或低维表示。 潜在表示可用于探索给定属性的不同标签之间的结构，例如使用相关矩阵，或通过考虑级联表示来研究不同属性之间的交互。

2. 2.

   不确定性评估——不确定性度量提供了评估模型普遍性的代理。 我们使用参考文献中提出的评估指标。 ⁵ 当已知与属性相关的其他协变量（例如药物属性的通路关联）时，它允许量化属性相对于其潜在表示的不确定性，并在 Biolord 包中提供其实现。 不确定性的定义是无法从变量中预测协变量（药物的途径）。 k-属性的潜在空间表示的NN图。 正式地，我们定义，

   $${u}_{i}=sum _{jin {{mathscr{N}}}_{i}}frac{1}{log dleft(i,jright)} times Hleft({C}_{{{mathscr{N}}}_{i}}right),$$

   在哪里 ({{mathscr{N}}}_{i}) 是值邻居的集合 我, d 是一个距离度量并且 H 是香农熵，并且 ({C}_{{{mathscr{N}}}_{i}}) 是与邻居相关的协变量向量 我 基于潜在表征。
3. 3.

   反事实预测——生物领主模块可以采用特定的细胞实例并修改其已知属性。 Biolord 获得的未知属性嵌入捕获了细胞特定的嵌入。 因此，当将细胞的测量特征以及感兴趣的已知属性的不同标签作为输入传递时，细胞特定的属性表示将保持不变（未知属性嵌入），并且只有修改后的已知属性的嵌入才会改变。 由于嵌入是提供预测的生成模块的输入，因此所有观察到的变化都是由已知属性的修改引起的。 为了在实践中获得反事实预测，我们采用一组参考单元，例如感染数据集中的控制单元，并使用它们的测量特征以及已知属性的任意组合（例如，考虑感染情况修改状态标签）：输入到经过训练的生物领主模型。 这使我们能够首先对看不见的生物状态进行采样，更重要的是获得一组受控样本，在这些样本中我们可以保证测量特征中所有观察到的变化都是修改属性的结果（图 1）。 1天）。

4. 4.

   特征与状态的关联——将反事实预测与统计测试配对，使我们能够恢复一组编码给定观察状态的特征（例如基因）。 在这里，我们显式地将修改后的属性与底层单元状态解耦； 因此，观察到的预测变化是由修改的属性引起的。

5. 5.

   分类——biolord-classify 模块可以扩展部分标记属性的标记。 这提供了数据的完整标记，原则上可以进一步检查并用作各种下游分析管道的输入。

数据集、培训和评估

科学丛3

sci-Plex 3 数据集¹⁶ 包含来自三种人类癌细胞系（A549、MCF7 和 K562）7,561 个基因的 649,340 个细胞的测量结果，以及四种不同剂量（10 nM、100 nM、1 μM 和 10 μM）的 188 种药物的扰动。 我们使用参考文献中提供的预处理 anndata 文件。 ⁵，下载自 https://f003.backblazeb2.com/file/chemCPA-datasets/sciplex_complete_middle_subset.h5ad。 在下载的 anndata 文件中，我们添加 RDKit 功能^{17 号} 使用 chemprop 包²⁷ 以及分配外拆分，保留九种未见过的药物进行验证——Dacinostat、Givinostat、Belinostat、Hesperadin、Quisinostat、Alvespimycin、Tanespimycin、TAK-901 和 Flavopiridol。

训练参数

我们在处理后的基因表达上训练生物领主模型。 我们使用 RDKit 化学信息特征嵌入药物^{17 号}，以及剂量作为有序属性。 细胞系作为分类属性传递。 我们使用权重和偏差²⁸ 用于实验跟踪和超参数调整。 补充说明中提供了超参数详细信息 3。

评估和基准

遵循参考文献提供的设置。 ⁵，我们使用决定系数来评估预测精度 r² （r² 分数），在模型的反事实预测和所有基因的真实测量之间计算。

包含的基准如下：

1. 1.

   朴素基线—— r² 在对照、未受干扰的细胞（每个细胞系）和相应的药物处理的细胞之间评估分数。

2. 2.

   化学注册会计师⁵– 使用报告的最佳超参数直接在单细胞数据上训练药物编码网络的独立设置⁵。

3. 3.

   chemCPA-预⁵– 预训练模型，药物编码网络通过批量 RNA 高通量筛选 (L​​1000) 进行训练²⁹。 chemCPA 的作者与我们分享了预训练模型⁵。 对所有对手参数进行了超参数调整。

4. 4.

   扰动网络⁸—该模型由三个网络组成：扰动表示网络、细胞表示网络和映射网络。 对于扰动表示网络，我们使用 ref 提供的预训练模型。 ⁸ 在 ZINC 数据集上进行训练³⁰。 其余网络按照 PerturbNet 在线 Github 存储库中提供的示例进行训练。 细胞表示网络是通过上述 anndata 文件进行训练的。 映射网络是根据两个训练模型提供的潜在表示进行训练的。

补充说明中提供了有关所有框架的更多详细信息 3。

遗传扰动

扰动序列（单基因）

Perturb-seq 数据集¹⁹ 包含 5,060 个基因的 65,899 个细胞的测量结果，其中包括 81 个单基因扰动和对照细胞。 我们使用 GEARS 提供的预处理后的 anndata⁶。 为了获得有意义的特征（代表遗传扰动），我们使用 GEARS 基因本体 (GO) 图中的扰动边缘。 GO 图最初是通过在共享大量 GO 项的基因之间添加加权边来生成的⁶。 最后，对于训练，我们仅考虑每个扰动和控制细胞的平均表达。

扰动序列（双基因）

Perturb-seq 数据集²⁰ 包含跨 5,045 个基因的 89,357 个细胞的测量结果，包括 131 个双基因扰动、105 个单基因扰动和对照细胞。 我们使用 GEARS 提供的预处理后的 anndata 对象⁶。 如上所述，我们利用 GEARS 的 GO 图来获取代表遗传扰动的有意义的特征。 对于训练，我们仅考虑单基因扰动和对照细胞的平均表达。 为了获得对两个基因扰动的预测，我们将表达差异近似为每个单基因扰动的预测差异之和。

训练参数

我们使用训练集中扰动的平均表达式来训练生物领主模型。 我们遵循 GEARS 中定义的设置，该设置考虑了在未见扰动集中不同的五种不同的训练-测试-验证分割。 对于双基因扰动设置，我们对扰动之间的五个分裂之一进行区分，其中两个、一个或零个双基因扰动在训练期间是看不见的。 我们使用 GO 术语特征作为模型的有序属性。 我们使用权重和偏差²⁸ 用于实验跟踪和超参数调整（补充说明 3）。

评估和基准

遵循参考文献中建议的程序。 ⁶，我们评估未见扰动预测中的归一化均方误差。 标准化是针对“无扰动”设置中的预测进行的，即预测执行扰动没有影响； 因此，未受扰动的细胞状态与受扰动后的细胞状态相同。

为了进行基准测试，我们将我们的性能与 GEARS 进行比较⁶，使用其再现性存储库中提供的设置运行评估（https://github.com/yhr91/gears_misc/blob/main/paper/fig2_train.py）。

时空单细胞图谱 疟原虫 肝脏阶段

研究疟疾寄生虫的肝脏阶段 疟原虫，阿夫里亚特等人。²¹ 在感染后的五个时间点（2、12、24、30 和 36）对数千个感染和未感染的肝细胞进行了分子特征分析。 我们从 Zenodo 下载了作者提供的预处理注释数据³¹。 数据注释包括以下内容：

• rough_time：表示收集细胞（或对照）时的 hpi 数。

• eta_normalized：基于分区标记基因的空间分区评分，用于将细胞分类为门静脉周围/中央周围。

• 伪时间：使用 Monocle 对受感染肝细胞 PBA 基因子集的标准化数据进行计算。

• 状态：通过肝细胞的FACS分选推断的感染状态。

• 流产：根据宿主转录组的聚类，将 36 hpi 时的细胞分类为流产/生产细胞。

训练参数

我们定义了两种生物领主设置，如下所述。 补充说明中提供了报告结果的超参数 4。

对完整数据集的感染状态分析

生物领主模型是针对注射小鼠（感染和未感染）以及对照小鼠（数据集不包括模拟和蚊虫叮咬样本）的肝细胞定义的。 作为输入，我们使用宿主转录组（限于原始出版物中使用的 8,355 个基因）以及状态分类（感染/未感染/对照）、空间区域（门静脉周围/中心周围）和时间（2、12、24、30、36） hpi 或控制）。

流产状态分类

在 24、30 和 36 hpi 时对受感染的肝细胞训练生物领主分类模型。 宿主转录组（仅限于高度可变的基因）以及空间区域（门静脉周围/中央周围）、时间（24、30 和 36 hpi）、stress_score（使用 scanpy 计算）³² 具有应激基因的函数“scanpy.tl.score_genes()”²¹）以及 36 hpi 的部分流产状态分类（流产/生产）。 我们引入了stress_score来解开原始出版物中报告的压力信号²¹，来自失败的签名。

报告摘要

有关研究设计的更多信息，请参阅 自然投资组合报告摘要 链接到这篇文章。