10分+ siVAE：单细胞转录组的可解释深度生成模型

导语

变分自动编码器（VAE）等神经网络对基因组数据的可视化和分析执行降维，但其可解释性有限：每个嵌入维度表示哪些数据特征是未知的。研究提出了siVAE，一种可通过设计解释的VAE，从而增强了下游分析任务。通过解释，siVAE还可以在没有明确基因网络推断的情况下识别基因模块和hub基因。本研究使用siVAE来鉴定其连通性与多种表型（如iPSC神经元分化效率和痴呆）相关的基因模块，展示了可解释的生成模型在基因组数据分析中的广泛适用性。

背景介绍

今天小编为大家带来一篇单细胞转录组可解释深度生成算法发表在10分+ Genome Biology的思路。题目为   siVAE: interpretable deep generative models for single-cell transcriptomes。  

数据介绍

本研究从ArrayExpress获得了胎儿肝脏图谱，从 10x 基因组学网站 获得了被称为“LargeBrainAtlas”的 130 万个脑细胞数据集，BrainCortex数据集（GSE126074），NeurDiff”数据集的 iPSC 神经元分化数据集，SEA-AD：西雅图阿尔茨海默病脑细胞图谱数据集。

研究设计

本研究提出了一种可扩展的、可解释的变分自动编码器（siVAE），它将变分自动编码器的非线性DR框架与线性PCA的可解释性相结合。siVAE是VAE的一种变体，它额外推断出用于解释细胞嵌入空间的基因组特征（基因或基因组区域）的特征嵌入空间。重要的是，通过使用非线性网络将单元和特征嵌入空间组合在一起，siVAE在不引入线性约束的情况下实现了可解释性，使其比LDVAE、scETM和VEGA更具表现力。与其他实现可解释的非线性DR的方法相比，siVAE要么更快，生成更准确的细胞低维表示，要么更准确地解释非线性DR，而不会引入线性限制或对先验知识的依赖。

siVAE推断单细胞基因组数据的可解释表示。a 对siVAE的输入是按特征分列的矩阵;这里显示的是八个基因的合成基因表达基质，其中四个受到严格调控（基因1-4），另外四个独立变化（基因5-8）。siVAE是一个神经网络，由一对编码器-解码器组成，它们共同学习细胞嵌入空间和特征嵌入空间。单元级编码器-解码器的作用类似于规范VAE，其中编码器的输入是单个单元c在所有输入特征（Xc,: ）上的测量值。细胞编码器使用输入细胞测量值来计算细胞嵌入空间中细胞位置的近似后验分布。特征编码器-解码器在所有输入训练单元 （X:,f）中获取单个特征f的输入测量值，并计算特征在特征嵌入空间中位置的近似后验分布。单元级和特征级编码器-解码器的解码器组合在一起，输出单元格c（Xc, f ）中特征f的表达水平。图b，d 从基因表达矩阵中学习的细胞和特征嵌入空间的可视化在图a。注意在d中，基因1、2、3和4的嵌入沿维度1都具有较大的幅度，但沿维度2没有，这表明基因 1、2、3和4解释了沿维度 1 的细胞嵌入空间的变化。基因5、6、7和8位于特征嵌入空间的原点，这表明它们不与其他特征共同变化。c 基因1的表达模式覆盖在细胞包埋空间中的细胞上。当从左到右检查细胞时，基因1的表达明显增加，这与显示基因1在维度1上具有大负载的特征嵌入空间一致。e 经过训练的siVAE模型可用于识别基因共表达网络中的枢纽和基因邻居，而无需明确推断共表达网络。

图1

结果解析

01siVAE准确生成细胞的低维嵌入

研究首先在细胞嵌入空间推理的背景下评估了siVAE，目标是生成细胞的低维表示，其中相同细胞类型的细胞聚集在一起。我们使用由177,376个细胞组成的涵盖40种细胞类型的胎儿肝脏图谱。将siVAE与其他可解释和不可解释的降维方法进行了基准测试。我们在5倍分层交叉验证框架中测量了每种方法的准确性，首先使用训练折叠来学习细胞嵌入空间，然后使用已知的细胞类型标签和嵌入空间内的细胞坐标训练k-NN（k=80）分类器。然后，我们将验证折叠中的单元格通过经过训练的单元格编码器以生成验证单元格嵌入，这些嵌入用于对验证折叠单元格进行分类。我们将更高的k-NN精度与更准确的细胞嵌入空间相关联，其中相同类型的细胞聚集在一起。 本研究将siVAE与VAE以及LDVAE和scETM进行了比较，其中所有四个VAE框架都使用相同大小的单元编码器解码器，潜在维度为2，VAE和siVAE使用相同的激活函数。总体而言，我们发现siVAE的细胞嵌入空间在精度上与VAE相当，这表明siVAE特征编码器-解码器的引入不会影响siVAE在其细胞嵌入空间方面的性能。siVAE细胞包埋空间的2D可视化显示，同类型细胞聚集在一起，与VAE的细胞包埋空间惊人相似（图2a）。此外，siVAE在整个胎儿肝细胞图谱（图2b）以及每种细胞类型上的VAE平衡分类准确性方面具有竞争力。因此，siVAE在生成细胞嵌入空间方面与VAE具有竞争力，但具有可解释性的额外优势，我们将在下面探讨。相比之下，LDVAE和scETM方法可像siVAE一样解释，但执行线性DR，其分类准确性显着降低（图2b），并且生成的可视化结果与VAE和siVAE相比，不同的细胞类型混合在一起更突出（图2a）。因此，LDVAE和scETM获得可解释性，但对细胞嵌入空间的准确性有很大的代价。 接下来，我们构建了一组siVAE的模型变体，以确定siVAE的哪些方面导致其优于LDVAE的性能。LDVAE与传统VAE大致相似，但有两个关键区别。首先，LDVAE解码器仅限于使用线性激活函数以实现可解释性;因此，LDVAE执行线性降维。其次，LDVAE损失函数在输入特征（基因）上使用负二项式或零膨胀负二项分布，而不是规范VAE中使用的高斯分布。原则上，与通常用于对数变换数据的高斯分布相比，NB或ZINB观察模型更适合单细胞转录组学数据。因此，我们构建了两种siVAE变体，称为siVAE-NB和siVAE-linear。siVAE-NB与siVAE相同，只是它对观察层使用负二项分布，同时在其解码器中保持非线性激活函数以实现非线性DR.siVAE-linear与siVAE相同，只是它限制特征和细胞解码器使用线性激活函数，如LDVAE，并且不实现可解释性项。图2b表明siVAE-NB的性能比具有高斯分布（siVAE）的相应模型差，这表明使用NB输出层不会导致更准确的细胞嵌入空间。siVAE线性比LDVAE更准确（图2b），表明siVAE的特征编码器-解码器总体上有利于降维。然而，siVAE-线性性能比siVAE更差，验证非线性激活函数有利于降维。 最后，siVAE本身允许在模型中进行批量校正，类似于其他方法。我们在iPSC神经元分化（NeurDiff）数据集上测试了我们的模型，其中253,381个iPSC衍生的细胞在开始分化为神经元之前和之后使 用10x Chromium测 序。我们特别关注分化前第11天的样品，因为我们观察到pool_id有很强的批次效应（图2c）。当我们训练siVAE并在训练期间提供批次信息时，消除了按批次的聚类，而按细胞类型进行的聚类仍然保留（图2d）。这些结果表明，siVAE是现有降维方法的可行替代方案，可用于执行降维、可视化和批量校正等常见任务。

图2

02siVAE 比现有的特征归因方法更准确、更快地解释细胞嵌入空间

在证明siVAE产生与规范VAE竞争的细胞嵌入空间之后，我们接下来验证了嵌入维度的siVAE解释是准确的。同样，我们将细胞嵌入空间的解释定义为大小为F×K的特征加载（或更一般地说，归因）V的矩阵，其中F是特征（例如基因）的数量，K是细胞嵌入维度的数量，V f, k 的大小表示原始数据空间中细胞嵌入维度k和特征f之间的关联强度。

除了通过设计构建可解释细胞嵌入空间的siVAE和LDVAE等方法外，文献中还有两种类型的特征归因方法可以帮助解释推理后的细胞嵌入空间。首先，通用神经网络特征归因方法可以量化神经网络的每个输出节点对网络每个输入节点（特征）的依赖性，包括DeepLIFT，显著性图，grad × input，积分梯度，Shapley值等方法。这些方法的优势之一是它们原则上可以应用于任何经过训练的神经网络，使它们具有高度的可推广性。其次，诸如基因相关性之类的方法已被专门开发用于解释任何DR方法的细胞潜在空间，包括那些不基于神经网络的方法，并且可以在学习细胞嵌入空间后应用。

我们首先使用神经网络特征归因方法作为金标准，将siVAE与基因相关性进行比较，因为它们已经在其他应用中得到了广泛的验证。图3a显示了 siVAE、基因相关性和三种神经网络特征归因方法（显著性图、grad × input和 DeepLIFT ）归因之间的平均成对相关性。siVAE负载与神经网络特征归因方法高度相关（中位数Spearman ρ=0.73，P=1.1×10−15），siVAE与DeepLIFT非常一致（中位数Spearman ρ=0.98，P=2.2×10−16）。相比之下，虽然基因相关性产生的特征归因在其参数选择中是一致的（中位数Spearman ρ=0.84，P= 3.10×10−22），但它们与神经网络特征归因方法（中位数Spearman ρ=0.11，P=2.1×10-6）和siVAE（中位数Spearman ρ=0.14，P=3.9×10−6）的相关性较差。这些结果表明，与siVAE相比，基因相关性在解释VAE架构的细胞嵌入空间方面不太准确。

图3

针对上述结果，我们将神经网络归因方法应用于siVAE解码器，生成真实特征归因。以前的工作建议将归因方法应用于编码器以提高执行速度。在这里，我们发现在siVAE编码器上运行归因方法会产生截然不同的解释，这些解释与解码器的解释存在强烈分歧，这表明对编码器的解释是不合适的。考虑到编码器的主要作用是计算每个单元的潜在嵌入的近似后验分布，而不是解码器，解码器负责将点从单元格嵌入空间直接映射到原始数据特征空间，这些结果是有意义的。因此，我们的结果表明，特征归因应该应用于VAE的解码器，而不是编码器。在我们解释细胞嵌入空间的实验中，很明显，许多神经网络特征归因方法的执行在计算上是昂贵的。由于这些特征归因方法分别对网络的每个嵌入维度或每个输出节点执行计算，因此在VAE解码器上运行时，它们的运行时间与嵌入维度或特征的数量成线性关系。我们通过检查每个批次的执行时间是否随时间不变，在较小的规模上证实了这一点。我们认为，执行时间在未来将变得更加重要，因为随着数据集中单元格数量的增加，嵌入维度的数量预计将更大，以适应数据集中更多的异质性。此外，对于诸如scATAC-seq之类的分析，特征的数量预计将很大，这些分析可以分析数十万个或更多的基因组区域。

为了确定每种方法的特征归因计算中最耗时的步骤，我们选择了LargeBrainAtlas数据集的一个子集，用于将嵌入维度的数量从20到512变化，并选择了BrainCortex数据集的一个子集，用于将特征数量从28k变化到240k。siVAE平均每个嵌入维度0.0073天（图3c），每10k个特征0.0027天（图3d），表明siVAE执行时间对细胞数量和特征都是稳健的。另一方面，我们发现当嵌入维度的数量或输入特征的数量很大时，神经网络归因方法的扩展性很好，但当它们都很大时，则不然。例如，DeepLIFT、grad ×input和显著性图分别在每个嵌入维度执行 0.014、0.0053 和 0.0012 天（图3c），但在输入特征数量方面缩放不良，每 10k 个特征分别执行 2.9、0.95和0.94天（图3d）。将渐变×输入和显著性图切换到正向模式导致输入特征数量的执行时间较快（每 10k 个特征分别为 5.3×10 天和 6.5×10 天）（图 3d），但导致嵌入维度数的缩放性较差（每个嵌入维度分别为0.18天和0.17天）（图 3c）。较慢的归因方法（如积分梯度和沙普利值）因其执行时间不可行而被排除在外。总之，神经网络归因方法的缩放性很差，无论是嵌入维度的数量还是输入特征的数量，这取决于使用的是正向模式还是反向模式。因此，如果特征的数量和嵌入维度都很大，这使得它们的执行时间相对于siVAE变慢，我们预计这种情况在不断增长的单细胞数据集中会越来越频繁地发生。

03共表达基因在特征嵌入空间中聚类

文献中广泛使用了PCA等线性DR方法的负载，以深入了解基因共表达网络（GCN）的结构。在这里，我们探讨了siVAE负载矩阵在多大程度上可以用来深入了解GCN结构。GCN是其中节点表示基因，边缘表示一对基因共表达的图形。GCN捕获细胞群中成对（或更多）基因之间基因表达测量的协变异。GCN之所以令人感兴趣，是因为它们可用于鉴定（1）细胞群特异性基因模块，代表高度共表达的基因组，因此可能以细胞类型特异性方式共同发挥作用，以及（2）基因枢纽，它们是与异常大量基因相连的基因（高度中心性）， 并且通常代表细胞中的关键功能基因。虽然GCN推理对于询问细胞中的基因调控模式很有价值，但GCN推理是一项众所周知的困难且容易出错的任务。 由于siVAE输入数据在所有特征中均匀居中和缩放，因此siVAE被迫学习输入特征之间的协变模式，从而可以从低维表示中准确重建输入数据。因此，很自然地会问，经过训练的siVAE模型是否可以深入了解训练数据的基因共表达网络结构，而无需明确的基因网络推断。 先前的工作表明，特征基因（通过PCA负载捕获的基因）代表基因共表达网络中的网络模块。我们假设通过siVAE特征负载捕获的siVAE基因也可能代表网络模块，此外，训练数据中的共表达基因在siVAE特征嵌入空间中近端。为了探索共表达基因组如何在特征嵌入空间中组织，我们构建了一个合成GCN，由五个社区组成，每个社区有50个紧密相关的基因，以及另外一组50个不相连的基因（图4a）。每个社区都遵循一个中心辐射模型，其中中心基因与社区中的每个其他基因相连，而社区中的每个基因又仅连接到中心。没有边缘连接来自不同群落的基因。基于该基因网络，我们对由5000个细胞和300个基因组成的单细胞基因表达数据集进行了采样。采样表达基质用于训练siVAE将基因嵌入其特征嵌入空间中。

图4

我们发现属于同一群落的基因在特征嵌入空间中共定位，但有趣的是，枢纽节点嵌入在其相应群落的不同位置（图4b）。我们推断，鉴于细胞嵌入空间的容量有限，siVAE倾向于专门维护有关每个枢纽的信息，因为它们的高度中心性。另一方面，同一群落中的非枢纽基因在特征嵌入空间中共定位，因为siVAE的有限能力迫使非枢纽基因在保留有关枢纽的信息下进行类似的预测。有趣的是，网络中的50个断开连接的基因在特征嵌入空间中紧密聚集但靠近起源，而作为群落一部分的基因聚类但距离起源更远。这可能是因为两个原因。首先，siVAE的特征编码器-解码器的KL发散项将倾向于将基因吸引到起源。其次，由于根据定义，断开的基因不会与其他基因共同变化，因此有关其表达模式的信息在压缩过程中往往会丢失，导致VAE倾向于预测解码器中该基因的平均表达水平（由于数据集中，所有断开连接的基因的平均表达水平将为0）。这反过来又鼓励特征嵌入位于原点，因为可解释性项鼓励特征嵌入的线性乘积与特征载荷来预测基因的表达模式，因此如果特征位于特征嵌入空间中的原点，则会导致预测的表达为 0。

我们还使用胎儿肝细胞图谱数据确认了共表达基因在特征嵌入空间中的聚类。与上述模拟不同，对于胎儿肝脏图谱，没有可用于鉴定真正共表达的基因，这些基因是同一潜在基因群落的一部分。因此，我们在具有40种细胞类型的整个胎儿肝脏图谱上训练siVAE，并认为相同细胞类型的标记基因是共表达基因的基本事实集。我们从MSigDB中选择了四种细胞类型进行可视化（图4c），基于观察结果，它们建立的标记基因集重叠最少（表明它们是不同的细胞类型），并且标记基因集与CellTypist数据库一致。对于每种选定的细胞类型，我们通过组合MSigDB中与靶细胞类型相对应的所有可用基因集来创建MSigDB元标记集。在siVAE学习的结果特征嵌入空间中，我们看到相同细胞类型的元标记倾向于按预期聚集在特征嵌入空间中（图4d）。我们还观察到，可视化来自密切相关细胞类型的标记会产生来自不同细胞类型的元标记的更显着混合，因为密切相关的细胞类型在元标记基因集中往往会有很大的重叠，我们也希望它们的元标记基因总体上也更紧密地共表达。我们的结果总体上表明，共表达基因倾向于在siVAE特征嵌入空间中共定位。

04siVAE隐式识别底层共表达网络中的基因hubs

我们观察到社区中的hub基因被siVAE区别对待，这促使我们假设我们可以从经过训练的siVAE模型中识别hub基因，而无需推断GCN。Hub基因通常在GCN推断后被识别出来，因为它们在网络结构和基因组本身方面都起着至关重要的作用，并且通常是与疾病相关的遗传变异的靶标。我们推断，由于hub基因与许多其他基因相连，siVAE更有可能将枢纽基因的表达模式存储在细胞嵌入中，以用于重建自身和其他基因表达模式。因此，我们在siVAE模型中使用基因特异性重建精度作为程度中心性的无GCN度量。作为程度中心性的基本事实测量，我们计算了每个基因预测基因组中所有其他基因表达水平的个体能力，推理hub基因应该预测网络中许多其他基因。我们观察到度中心性和siVAE重建精度之间存在很强的关系，正如预期的那样。注意对于GCN相关分析，我们将潜在嵌入的数量设置为64，而不是用于可视化嵌入空间的2，因为这里我们不需要直接可视化嵌入空间。 图5a比较了siVAE对度中心性的估计与使用许多现有GCN推理算法推断的GCN估计的度中心性的精度。总体而言，siVAE的预测度中心性和地面真实中心性之间的相关性最高（Spearman ρ=0.90，P=2.2×10−16），明显大于其他方法（中位数Spearman ρ=0.36，P=9.0×10−11）。对于siVAE，也观察到很强的相关性，可解释性项（γ）降低到0。作为单独的性能比较，我们查看了前20个中心的平均度中心性，其中较大的值表示中心更强大。siVAE前20个中心的平均度中心性为0.092，大于GCN推理方法，GCN推理方法的前20个预测中心的平均度中心性为0.074（图5b）。在比较2000个变化最大的基因时，结果也是一致的。这些结果表明，使用siVAE，我们可以比首先推断GCN然后识别枢纽基因更准确地识别高度中心性基因。

图5

05在降维和网络推理方法之间观察到基因邻居的系统差异

本研究还探索了可以在多大程度上识别在GCN中共享边缘的相邻基因，而不必明确推断GCN。基因邻居往往具有相似的功能，相互作用，和/或属于同一基因群落。因此，基因邻居的鉴定有助于识别细胞中的协同功能基因。 在这里，我们假设我们可以使用经过训练的siVAE模型直接识别基因邻居，而不必首先推断明确的GCN。GCN 推理方法通常在网络中的节点对之间输出边权重，其中较大的权重对应于两个节点在底层 GCN 中共享一条边的可能性更大。对于siVAE，我们计算了两组不同的边缘权重：（1）siVAE-Euc，其中两个基因之间的边缘权重设置为它们在特征嵌入空间中的欧氏距离，较小的距离对应于更近的距离;（2）siVAE-GCN，我们首先从训练的siVAE模型中采样一个新的scRNA-seq数据集，该数据集与训练数据的大小相匹配，然后在采样的scRNA-seq数据集上运行GCN推理方法（ARACNE，MRNET，CLR和GRNBOOST2）以计算基因之间的边缘权重。为了定量评估使用每种方法进行邻域识别的准确性，我们测量了给定查询基因在通过最接近查询基因边缘权重排名的20个基因的表达水平预测时解释的方差百分比。直观地说，与其他基因相比，查询基因的真正邻居应该更能预测查询基因的表达水平。在我们的评估中，我们只考虑了152个查询基因，这些基因被预测在所有测试方法中具有高度中心性。 总体而言，大多数方法识别出的邻居同样可以预测152个查询基因的表达水平，除了scETM（图5c）。不包括LDVAE和ARACNE，每种方法解释的方差中位数百分比为79.9%，±0.84 s.d。排除LDVAE和ARACNE，方法之间解释的方差百分比的成对差异平均仅为0.013%。值得注意的是，我们观察到LDVAE和ARACNE解释的方差百分比较低（平均而言，解释的方差分别为77.2%和78.3%）。LDVAE较差的结果与上述较差的分类性能结果一致。 当考虑通过不同方法选择的邻居的重叠时，令人惊讶的是，降维方法（scVI，siVAE，LDVAE）作为一个组强烈聚类，而基于GCN推理的方法作为一个组强烈聚类，但这两个组之间的重叠明显较少（图5d）。这部分令人惊讶，因为邻域集都具有大致相同的预测性能（图5c），这表明DR方法正在系统地识别不同的邻居，这些邻居与GCN方法识别的邻居集一样共同表达。特别是，考虑到siVAE-GCN使用GCN推理方法识别基因邻居，但应用于siVAE生成的数据集（而不是原始训练数据集）。图5d表明，与GCN推理方法相比，由siVAE-GCN识别的邻域基因仍然与siVAE识别的邻域基因更相似，这表明DR方法识别的独特邻域是DR方法学习的共表达模式的属性，而不是由于识别邻域基因的方法。如果我们考虑相邻集之间表达的平均成对相关性，而不是测量基因的直接重叠，DR和GCN方法之间的不良重叠也成立（图5e）。更具体地说，与基于神经网络的邻域集之间的平均皮尔逊相关系数（平均皮尔逊ρ=0.46）相比，GCN定义的邻域集之间的平均皮尔逊相关性（平均皮尔逊ρ=0.67）具有更高的平均皮尔逊相关性。DR和GCN邻居集之间的平均相关性也很低（平均Pearsonρ = 0.39） 。因此，我们的结果表明，由于GCN和降维识别的邻域集是系统不同的，但对相邻基因的预测大致相同，那么这两种方法都应该用于在网络中寻找共表达的基因。

06iPSCs中线粒体基因的共表达与神经元分化效率有关

接下来，我们假设可以使用siVAE隐式学习多个细胞群中的GCN，从而找到网络结构与细胞群水平表型之间的关联。由于稳健的单GCN推理已经具有挑战性，因此尚未对比较多个GCN的方法进行广泛的研究。如前所述，NeurDiff 数据集包括在分化前（11 天，主要由两种祖细胞类型，中脑底板祖细胞 （FPP） 和增殖 FPP （P-FPP） 细胞）以及开始分化为神经元后（第 30 天和第 52 天）分析的 215 个 iPSC 细胞系中收集的 scRNA-seq 数据。通过计算第52天被鉴定为成熟神经元的测序细胞的比例，每个细胞系都可以估计它们分化为神经元的效率。发现效率是高度可重复的，最初的研究发现第11天多能细胞的单基因表达水平与效率之间存在有希望的关联。在这里，我们假设第11天iPSCs中基因亚群的共表达模式也与分化效率有关（图6）。我们分别对P-FPP和FPP细胞进行了GCN分析。

图6

讨论

降维是许多单细胞基因组分析任务中关键的第一步。siVAE代表了使非线性降维方法可解释的重大进步，从而提高了对驱动输入数据变化的因素的洞察力。我们展示了siVAE的可解释性特征如何在不需要明确的网络推断的情况下识别基因共表达网络的新特征。这些范围从识别基因中心和基因邻居到识别其网络连接本身与分子和临床表型相关的基因模块。siVAE具有可扩展性，并广泛适用于单细胞转录组学以外的其他单细胞数据模式。