论文总结1、提出VEGA基因生物信息增强的变分自编码器。线性解码器稀疏由生物网络提供信息。摘要变分自编码器等深度学习架构彻底改变了转录组学数据的分析。然而这些变分自编码器的潜空间几乎没有解释性。为了提供更多生物学见解我们引入了一种新型稀疏的变分自编码器架构VEGAVae Enhanced by Gene Annotations其解码器接线灵感来自先验的生物抽象能够直接解释潜在变量。我们通过整合多种生物抽象来源如通路、基因调控网络和细胞类型身份展示VEGA在多样生物学环境中的可解释性和灵活性。我们展示了该模型能够重现细胞特异性反应机制、主调控因子的状态并共同研究发育细胞中的细胞类型和细胞状态。我们设想该方法可作为开发和药物治疗实验等背景下的解释性生物学模型。引言单细胞RNA测序scRNA-Seq技术的最新进展使得细胞状态的表征达到了前所未有的规模和分辨率。在众多广泛使用的单细胞复杂转录组模式分析框架中人工神经网络ANN如自编码器AEs2已成为强大的工具。AE是旨在将输入数据集转换为解码表示的神经网络同时最小化信息丢失3。其架构设计的多样性使AE适合应对scRNA-Seq分析中诸多重要挑战如降维4、聚类5和数据补值6。近年来深度生成模型如变分自编码器7VAEs已被证明在单细胞转录组的概率建模中极为有用如scVI或scGen 8–10。虽然深度生成模型在其专门建模任务中表现出色但它们通常缺乏可解释性因此无法提供生物学意义上的潜在转录组表示。例如用 scGen 提取的潜在扰动载体不能直接关联基因模块变异10。这里我们提出了VEGA基因注释增强的VAE这是一种变分自编码器其线性解码器稀疏由生物网络提供信息。该模型旨在提供一个可解释的潜在空间以表示各种生物学信息例如生物通路的状态或转录调控因子的活性。具体来说我们模型的范围有两个方面1在可解释的潜在空间上编码数据2推断样本外数据的生物信息。结果VEGA的网络结构设计为了创建易于解释的VAE我们提出了一种受生物先验启发的新架构如图1A所示。图1。设计一种具有可解释潜在空间的新型VAE架构。A VEGA模型概述。VEGA由深度非线性编码器和掩蔽线性解码器组成能够在先验指定的可解释潜空间中编码单细胞转录组数据。B VEGA潜空间的UMAP嵌入保留了Kang等人的生物信号。PBMC数据集14。C推断干扰素-α/β信号通路活性将受刺激细胞与对照群体分离。D双变量因子图显示模型恢复色氨酸分解代谢活性的能力这是一种先天免疫细胞对扰动的特异性反应。E火山图显示受刺激和控制先天免疫细胞之间活性生物因素的差异。红点表示具有|loge贝叶斯因子|3且潜空间至少为5的均绝对差MD。F VAE贝叶斯因子与GSEA的比较 -log10FDR。点的大小表示基因组的大小。红色、蓝色和紫色象限分别对应我们模型独有的显著命中点、GSEA独有且两者共有的命中点。在许多标准VAE实现中编码-解码器架构的信息瓶颈通常将潜在变量建模为多元正态分布。尽管VAE潜在变量对输入数据提供了高度有用的表示但通常难以解释。为解释VAE潜在变量Svensson等人提出研究对原始变量的相应权重类似于PCA11等标准因子模型所提供的可解释性。虽然这种方法在解释单细胞生物学方面提供了宝贵见解但需要对解码器权重进行进一步的统计分析。我们模型的生成部分译码器直接关联潜在变量生物抽象如基因调控网络与原始基因特征当且仅当基因k作为生物抽象l的一部分先验注释时基因k与潜在变量zl之间存在依赖关系。这依赖于一个基本假设一组已知的生物模块支撑并协调给定scRNA-Seq数据集的基因表达模式。注释可以来自经过策划的数据库如MSigDB12、推断基因调控网络13或任何其他预定义的基因集。这种设计有两个主要优势1潜在变量先验初始化因此可直接解释为生物模块的活性;2生物模块规范的灵活性使其可推广到不同的生物抽象如通路、基因调控网络甚至细胞类型、方法。在初始化此类模型架构时可以建模生物模块与推理网络编码器、生成网络译码器中原始特征之间的依赖关系或两者兼有。我们选择使用深度编码器、稀疏线性解码器结构可以基于提升模型推理能力的目标当近似生物模块的真实后验分布Pzi|习时可以捕捉可解释潜在空间中的更复杂模式而潜在变量与原始基因特征之间的联系则由神经网络中的稀疏线性解码器部分保证补充图1。在可解释的潜在空间上重现生物信息为验证我们提出的模型重现生物通路状态的能力我们将该模型应用于已发表的外周血单核细胞PBMCs数据集该数据集由趋化因子干扰素β14刺激METHODS。我们首先证明模型能够利用反应组通路作为潜在变量捕捉细胞类型和刺激状态图1B。我们进一步研究了单个通路的激活状态具体证明干扰素-α/β信号活性会分离受刺激和初生细胞证实模型能够在其潜在空间捕捉通路活性图1C图1D。我们还进一步研究了参与干扰素诱导免疫细胞激活的其他已知生物通路发现某些细胞过程具有特定细胞类型的激活。例如色氨酸对干扰素的分解反应将先天免疫细胞树突状细胞、FCGR3A单核细胞、CD14单核细胞与适应性免疫细胞NK细胞、T细胞CD8、T细胞CD4、B细胞分离开来如先前研究的那样图1D。这些结果共同展示了模型能够将细胞投射到可解释的空间中从而在细胞过程层面研究特定细胞类型的模式。通常情况下对刺激的预期反应往往并非事先已知。差异基因表达已成为研究样本组间基因表达模式差异的标准工具。该程序也适用于基因表达建模的深度生成网络9。我们提出类似方法基于模型潜伏空间研究生物模块活性差异我们可以提出备选假设并利用模型的变分后验分布对应贝叶斯因子BF方法计算其后验概率的比值。当应用于刺激组与对照组的先天免疫细胞时我们表明本方法发现预期在被刺激组中被激活的通路干扰素信号、色氨酸分解代谢具有显著性|logeBF| 3图1E。我们将差异因子活性方法的结果与标准GSEA工具包的假发现率FDR值进行了比较方法图1F。虽然两种方法均在受刺激组中激活了干扰素-α/β信号通路但GSEA未发现先天免疫细胞中色氨酸分解代谢激活见图1F。总体而言VEGA似乎较少受到GSEA固有的基因集大小偏差问题影响图1F补充图2鉴于数据集的生物学特性VEGA显示出更具上下文特异性且显著的命中。细胞系药物治疗生物反应的大规模研究接着我们研究了模型是否能捕捉并重现癌症细胞系大规模实验中的药物反应模式如近期实验方案如MIX-Seq 17。为此我们收集了97个癌细胞系在5种不同条件下的单细胞数据24小时DMSO治疗对照、24小时曲美替尼治疗MEK抑制剂、24小时达曲替尼治疗突变BRAF抑制剂、24小时纳维托克拉斯治疗Bcl-2抑制剂和24小时BRD3379治疗未知作用机制的工具化合物MoAMETHODS。我们通过结合药物治疗数据集和对照组DMSO数据集初始化潜伏空间训练了每种不同药物治疗的模型共4个模型。总体而言每个模型均能重现细胞系对各自的生物反应治疗图2A补充图3。对于曲美替尼而言G2M检查点活性的重要变化治疗条件的下降与MEK抑制剂19,20的预期MoA相符图2B。接下来我们试图研究是否能重现对照组与处理条件下的生物反应模式。对于每一对我们计算了模型中每个潜在变量通路的贝叶斯因子。所得热图可用于理解和解释所有实验条件下的反应模式图2C。在研究每个数据集的低维嵌入时发现曲美替尼在所有研究药物中产生了最强的转录反应。值得注意的是曲美替尼特异性干扰素α和干扰素γ反应被我们的模型正确重现正如之前的实验21和原始MIX-Seq作者17所示。此外我们发现达布洛芬尼治疗的BRAF突变黑色素瘤细胞系表现出强烈的转录反应与曲美替尼细胞系聚集如MIX-Seq研究中报道的那样图2.C补充图3。总体而言这里呈现的结果与该数据集上之前的基因集分析结果一致展示了我们模型在大规模实验中重现药物反应模式的能力。胶质母细胞瘤的基因调控分析揭示了肿瘤细胞的分层现象如前所述我们模型的一个主要优势是潜在空间规范的灵活性任何生物模块都可以与基因集相关从而编码在我们的模型中。生物学中最重要的方面之一是转录因子对基因表达的严格调控22。分析这些主调控因子的活性对于理解细胞类型或疾病等生物状态非常重要因为它们的活性失调会对基因表达程序和表型产生显著影响23,24。为此我们研究了使用主调控因子指定模型潜伏空间是否有助于理解单细胞胶质母细胞瘤GBM数据集中潜在的基因调控网络GRN25。我们使用了24年报告的GBM ARACNe 13网络来指导模型的结构设计。具体来说转录因子用于初始化模型的潜伏空间转录因子-靶点连接初始化解码器接线。训练模型后我们表明预注释的细胞类型可以在潜伏空间中被重现图2D。我们研究了STAT3和OLIG2的活性这两种分别是间充质MES和前神经PNGBM亚型的主调控因子。我们确认它们在潜伏空间中推断的活性在肿瘤细胞中存在反相关性图2E。此外作为寡突胶质细胞分化的已知主调控因子OLIG2被推断为在寡突胶质细胞前体细胞OPCs中被激活。这些结果表明我们方法在潜在空间和解码器布线的规范上具有灵活性和可扩展性。此外这些结果也证明了我们模型能够忠实地还原生物信息。结合细胞类型和细胞状态表征精炼皮层类器官发育分析现代细胞生物学的一个重大挑战是识别和定义细胞类型和细胞状态以便在共同词汇下系统性地研究稳态和疾病发展。为此我们研究了是否能将细胞类型以标记基因列表定义的标记集形式和细胞状态反应组通路的信息结合在模型的潜在空间中以实现细胞类型和细胞状态的解离对应表示。我们将模型应用于Field等人27页中皮层类器官早期发育期间的细胞数据集包括潜伏空间中研究中定义的主要细胞类型信息图3A。训练后每个细胞类型标记集的活动能够正确区分其对应细胞类型正如原作者所注图3B-D。此外在针对每种细胞类型群体的“一与静止”差分因子分析中对应标记集的活性显著富集|logeBF| 3这可能为研究人员利用细胞类型标记集的潜在活性来注释未知簇提供了机会图3E。为了研究我们的模型是否能将细胞类型身份与细胞状态如分裂细胞群与静止细胞群体区分开来我们将数据集预测为两个组成部分1神经上皮标记集活动一种早期脑前驱和2细胞周期有丝分裂通路活动图3F。如前所述神经上皮标记集活动将神经上皮细胞与数据集其他部分区分开来而细胞周期有丝分裂通路的活性则区分了两个前驱群体桡神经胶质细胞和神经上皮中静止细胞与活跃分裂细胞。为验证被识别为分裂细胞的增殖我们研究了细胞周期有丝分裂通路活动与MKI67基因表达之间的相关性MKI67基因是典型的增殖标志细胞周期有丝分裂通路集中不存在外部验证器见图3G。总体而言MKI67的表达与细胞周期有丝分裂通路的推断活动高度相关R2 0.67。综合来看这些结果展示了我们模型在不同细胞群体中联合推断细胞类型和细胞状态身份的潜力因为结合模型潜在空间中不同的信息来源通路、主调控因子、细胞类型标记可以揭示单细胞身份的不同方面。推理过程推广到样本外数据我们假设模型可用于推断训练时未见的数据的可解释潜在表征我们称之为样本外数据。为确认其可推广性我们在两种情况下评估模型1推断的“生物学泛化”即在Kang等人PBMC数据集训练时遗漏了某一对细胞类型/状态类似于scGen10数据集内设置以及2推理的“技术泛化”即在Kang等人14 PBMC数据集上训练的模型我们称之为研究A在另一个模型上进行评估仅包含Zheng等人对照细胞的PBMC数据集28研究B数据集间设置。对于“生物学推广”我们首先检查了样本外刺激的CD4 T细胞中干扰素-α/β信号通路活性分布与样本内CD4 T细胞推断的活性相符图4A。为了更系统地比较样本外与样本内细胞之间推断的潜伏空间我们采用了贝叶斯微分因子活性过程METHODS分别对1个受刺激的样本内细胞和针对特定细胞类型的对照单元用整个数据集训练的模型和2刺激同一细胞类型的样本外细胞和对照单元训练模型去除一对细胞类型/条件并检查前50个差异激活因子的重叠程度图4B。结果显示样本内与样本外差异激活因子之间存在一致性平均重叠率为72%。为进一步评估数据重建能力我们在样本内和样本外设置中测量了原始数据与解码数据之间的R2见图4C。我们表明在样本外环境中R2仅有微小下降证实了模型能够推广到类似实验环境中产生的未见数据。对于“技术推广”我们再次检查了研究B编码的对照CD4 T细胞的干扰素-α/β通路活性分布与研究A对照的CD4 T细胞一致见图4D。我们还调查了研究A对照细胞在“一对一”差异条件下每种细胞类型的前50个差异因子是否与研究B对照细胞的类似程序重叠图4E。我们显示研究A中平均67%的前50个差异因素与研究B重叠表明模型可以在训练时未见的研究中推广。随后我们研究模型是否能利用推断出的潜空间准确重建两项研究的原始表达特征。我们表明研究B的原始细胞与重建细胞之间的R2值虽低于研究A但在大多数细胞类型中研究A与研究B表达谱的基线相关性有所改善见图4F。我们认为这些结果共同展示了模型能够推断训练时未见数据的可解释潜在潜在空间但部分局限性来自潜在的批处理效应。通过重新训练模型来轻松缓解这个问题因为这个过程相当快。讨论本研究引入了一种新型VAE架构VEGA其解码器接线灵感来自已知生物学用于推断单细胞层面各种生物模块的活性。通过将单细胞转录组数据编码到事先指定的可解释潜空间中我们的方法为分析不同背景下各种生物抽象的活性提供了一种快速高效的方法。我们展示了其在理解特定细胞类型群体对不同扰动反应中的应用提供了对生物模块活动的可解读见解。我们设想我们的模型可用于根据癌症中通路表达优先排序药物因为研究特定细胞群体的反应有助于理解药物的敏感性和耐药性。将药物反应预测模型与我们这样的解释模型整合可能对设计新颖的治疗策略非常有益。此外我们推测我们可解释的VAE将有助于理解一般生物学问题。 潜空间的灵活性为分析生物模块的活性铺平了道路如通路、转录调控因子甚至细胞类型特异性。本研究中展示了模块。我们还进一步证明我们的模型可用于共同探索不同细胞亚群的细胞类型和细胞状态无论是否涉及扰动实验。此外解码器连接的权重提供了潜在变量与原始特征之间关系的直接解释性。例如它可以用于对比推断基因调控网络中的相互作用置信度或通过数据驱动的方式对基因在特定生物模块中的重要性进行排名。 当前架构的明显局限在于模型中稀疏的单层解码器。事实上这种架构设计阻碍了普适性和鲁棒性的进一步提升。因此我们的VAE的生成能力受到了削弱。例如理论上我们的模型可以用潜在向量算术进行可解释的响应预测类似于 scGen10但我们方法有限的生成能力牺牲了预测性能以换取可解释性。因此我们的设计选择在潜在空间的生物学解释性上妥协了这些方面。我们相信网络生物学中的先进见解例如能够更全面描述调控机制的多层GRNs能够缓解这些局限。这将为对特定细胞群体进行靶向的、计算机内激活和抑制生物程序的可能性打开了可能以研究其对发育或疾病进展的影响。 另一方面线性解码器的硬编码连接在上下文需要时没有空间纠正生物因素的先验知识这在其他潜在变量模型如f-scLVM 29中存在。事实上从现有数据库如MSigDB获得的先前生物学知识可能不完整或非上下文特定因为额外的未注释基因在某些生物学因素中可能起重要作用。我们认为通过使用标准正则化方法对解码器的权重进行正则化而非掩蔽这一点可以更容易改善。我们把这个问题留给后续探索。方法VEGA的网络结构我们在此介绍用于生物模块分析的变分自编码器VAE架构。我们的VAE是一个深度生成模型旨在最大化单细胞数据集X在生成过程下出现的可能性过程7,10描述为其中θ是神经网络的可学习参数。这里我们选择让潜在变量Z明确表示生物学抽象如通路、基因调控网络或单元类型标记集。为实现这一点我们将VAE架构中的解码器部分修改为单层、掩码、线性解码器。具体来说该层连接潜在节点zj与输入基因特征之间的连接通过掩码M指定M是一个二元矩阵其中训练过程中我们会归零与掩码权重相关的梯度使反向传播仅适用于来自先验注释的权重其余权重保持为0。此外解码器被限制为正权重w ≥ 0以保持对生物模块活动方向性的可解释性。在模型生成部分明确指定了基因与潜在变量之间的联系后我们确保潜在空间代表生物模块的活动及其可解释性。我们选择将潜在变量Z的集合建模为多元正态分布由我们的推理网络参数化参数可学习φ如下这种变分分布的选择很常见并且在以往的单细胞研究中已被证明效果良好 9,10。训练模型时目标是下界证据ELBO类似于许多标准VAE实现 7,10其中变分分布的期望可以通过对一小批数据的蒙特卡洛积分近似且当我们设先验时库尔布莱克-莱布勒散度项的解是闭形式。在采样变分分布时使用了重参数化技巧7以便在训练模型时应用标准反向传播。 为了保留预先注释的生物网络中不存在的基因信息我们在模型的潜在空间中增加了额外的全连通节点。这有两个效果1它能够模拟未注释基因的表达这对于训练期间数据的良好重建至关重要;2它可以帮助捕捉现有生物因素无法解释的数据额外方差显著提升模型的训练效果。贝叶斯微分因子活动研究两组细胞间因子激活的差异及其意义具有重大意义。为此我们从贝叶斯差异基因中汲取灵感表达式程序在9中引入并提出了类似的差分因子分析方法。对于给定的因子k一对细胞xa xb及其对应的基群IDsa sb例如两种不同的处理条件我们有两个互斥的假设这直观上可以看作是测试一个细胞的生物因子激活平均值是否高于其他细胞期望值代表经验频率。我们通过研究定义为对数贝叶斯因子K来评估最可能的假设与9类似假设单元格独立我们可以计算从各组随机抽样的多个单元格对的平均贝叶斯因子。这有助于我们判断某一因子在某一群体中被激活频率更高。通过本文我们认为如果K的绝对值大于3相当于比值比为≈20因子会显著差异激活。数据集与预处理Kang等人数据集Kang等人14的数据集包含两组PBMCs一组对照组一组用干扰素-β刺激组。我们选择使用与 scGen 作者 10 所描述的相同的预处理步骤。简而言之细胞通过与八个原始细胞类型簇中最大相关性进行注释平均取前20个簇基因的平均值。随后对数据进行了筛选去除表达基因少于500个和表达在5个或更少细胞中的基因。然后对每个单元格的计数进行了归一化和对数变换我们选择了顶部6998个高度变异的基因最终数据集为18,868个细胞。原始数据可在GSE96583获取。Zheng 等人数据集Zheng 等人.28数据集包含来自健康供体的3000个PBMCs。过滤后每个细胞的计数被归一化并进行对数转化。然后我们对基因进行了子集使用Kang等人PBMC数据集中的相同6998个基因。最终数据集包含2623个细胞和6998个基因。原始数据可在 https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.1.0/pbmc3k 年获取。MIX-Seq 数据集MIX-Seq17数据集来自 https://figshare.com/s/139f64b495dea9d88c70我们利用实验3的数据收集了足够的单元实现模型的平滑训练。对于5个可用数据集97个分别用DMSO、曲美替尼、达布拉非尼、纳维托克莱克斯和BRD3379处理的细胞系我们切除了表达基因少于200个的细胞以及表达少于3个细胞的基因。然后我们对每个单元格的计数数进行了归一化并对数变换了数据。最后每个药物治疗实验数据集与对照数据集DMSO处理结合提取出前5000个高度变异基因。最终数据集规模为曲美替尼DMSO数据为16732细胞4999基因;达布拉非尼DMSO数据为16942细胞5000基因;NavitoclaxDMSO数据为14507细胞5000基因;BRD3379DMSO数据为15304细胞5000基因。Darmanis 等人数据集Darmanis等人提供的原始胶质母细胞瘤数据25来自 http://www.gbmseq.org/ 并预处理如下我们移除表达基因少于200个的细胞以及表达在3个或更少细胞中的基因。每个单元格的计数被归一化然后对数变换数据。最后我们将转录组限制为前6999个高度变异基因。最终数据集共计3566个细胞。原始数据可在GSE84465获取。Field 等人数据集Field等人27的皮层类器官数据处理方式与胶质母细胞瘤数据集类似。经过归一化和高度变异基因选择后数据集共有4378个细胞6999个基因。原始数据GSE106245可用。VEGA潜空间基因注释的选择在初始化模型的潜在空间时我们选择使用分子签名数据库MSigDB12中的预注释基因集。特别是我们选择使用标志基因集注释50个基因集或Reactome数据库674个基因集。该研究的MIX-Seq分析部分采用了反应组MIX-Seq分析部分采用了标志性H。对于胶质母细胞瘤细胞的基因调控网络分析我们从胶质母细胞瘤的整体RNA-Seq样本中推导出ARACNe13,30网络。具体来说该网络来自一篇先前发表的论文31并重新用于研究胶质母细胞瘤单细胞转录组学谱。对于皮层类器官分析中的细胞类型标记基因我们联系了作者以获取用于注释这些细胞类型的最相关基因。标记列表作为补充表提供。可视化降维方法在可视化数据集时我们使用了UMAP Python软件包中实现的UMAP算法32。我们使用默认参数除了min_dist参数设置为0.5。我们还使用了 sklearn python 包中实现的 t-SNE 33默认参数。与GSEA的比较我们用Python中的gseapy包中的预秩函数运行了基因集富集分析12GSEA。简而言之我们使用Wilcoxon秩和检验计算了对照组和处理组每个基因的差异表达得分该检验体现在Scanpy包34的 rank_genes_groups 功能中。我们根据基因的测试统计数据对基因进行了排名并以以下设置运行了GSEA预排序最小基因集大小为5最大基因集大小为1000以及1000个排列。我们根据FDR对基因组进行了排名并在FDR≤0.05时考虑了显著的命中。当GSEA返回的FDR为0时我们用1e-5替代以避免对数计算时出现数学错误。评价指标我们计算了轮廓分数以评估模型潜伏空间中的生物分离情况。我们利用潜空间中的欧几里得距离计算每个单元格i的轮廓系数定义为其中ai和bi分别是单元i的平均簇内距离和单元格的最近簇平均距离。我们使用了Kang等人14中的刺激标记或细胞类型标记来评估模型潜伏空间的生物学相关性见补充说明1。sklearn package35 silhouette_score实现用于计算。
受生物网络启发的可解释变分自编码器
发布时间：2026/5/19 5:00:32
论文总结1、提出VEGA基因生物信息增强的变分自编码器。线性解码器稀疏由生物网络提供信息。摘要变分自编码器等深度学习架构彻底改变了转录组学数据的分析。然而这些变分自编码器的潜空间几乎没有解释性。为了提供更多生物学见解我们引入了一种新型稀疏的变分自编码器架构VEGAVae Enhanced by Gene Annotations其解码器接线灵感来自先验的生物抽象能够直接解释潜在变量。我们通过整合多种生物抽象来源如通路、基因调控网络和细胞类型身份展示VEGA在多样生物学环境中的可解释性和灵活性。我们展示了该模型能够重现细胞特异性反应机制、主调控因子的状态并共同研究发育细胞中的细胞类型和细胞状态。我们设想该方法可作为开发和药物治疗实验等背景下的解释性生物学模型。引言单细胞RNA测序scRNA-Seq技术的最新进展使得细胞状态的表征达到了前所未有的规模和分辨率。在众多广泛使用的单细胞复杂转录组模式分析框架中人工神经网络ANN如自编码器AEs2已成为强大的工具。AE是旨在将输入数据集转换为解码表示的神经网络同时最小化信息丢失3。其架构设计的多样性使AE适合应对scRNA-Seq分析中诸多重要挑战如降维4、聚类5和数据补值6。近年来深度生成模型如变分自编码器7VAEs已被证明在单细胞转录组的概率建模中极为有用如scVI或scGen 8–10。虽然深度生成模型在其专门建模任务中表现出色但它们通常缺乏可解释性因此无法提供生物学意义上的潜在转录组表示。例如用 scGen 提取的潜在扰动载体不能直接关联基因模块变异10。这里我们提出了VEGA基因注释增强的VAE这是一种变分自编码器其线性解码器稀疏由生物网络提供信息。该模型旨在提供一个可解释的潜在空间以表示各种生物学信息例如生物通路的状态或转录调控因子的活性。具体来说我们模型的范围有两个方面1在可解释的潜在空间上编码数据2推断样本外数据的生物信息。结果VEGA的网络结构设计为了创建易于解释的VAE我们提出了一种受生物先验启发的新架构如图1A所示。图1。设计一种具有可解释潜在空间的新型VAE架构。A VEGA模型概述。VEGA由深度非线性编码器和掩蔽线性解码器组成能够在先验指定的可解释潜空间中编码单细胞转录组数据。B VEGA潜空间的UMAP嵌入保留了Kang等人的生物信号。PBMC数据集14。C推断干扰素-α/β信号通路活性将受刺激细胞与对照群体分离。D双变量因子图显示模型恢复色氨酸分解代谢活性的能力这是一种先天免疫细胞对扰动的特异性反应。E火山图显示受刺激和控制先天免疫细胞之间活性生物因素的差异。红点表示具有|loge贝叶斯因子|3且潜空间至少为5的均绝对差MD。F VAE贝叶斯因子与GSEA的比较 -log10FDR。点的大小表示基因组的大小。红色、蓝色和紫色象限分别对应我们模型独有的显著命中点、GSEA独有且两者共有的命中点。在许多标准VAE实现中编码-解码器架构的信息瓶颈通常将潜在变量建模为多元正态分布。尽管VAE潜在变量对输入数据提供了高度有用的表示但通常难以解释。为解释VAE潜在变量Svensson等人提出研究对原始变量的相应权重类似于PCA11等标准因子模型所提供的可解释性。虽然这种方法在解释单细胞生物学方面提供了宝贵见解但需要对解码器权重进行进一步的统计分析。我们模型的生成部分译码器直接关联潜在变量生物抽象如基因调控网络与原始基因特征当且仅当基因k作为生物抽象l的一部分先验注释时基因k与潜在变量zl之间存在依赖关系。这依赖于一个基本假设一组已知的生物模块支撑并协调给定scRNA-Seq数据集的基因表达模式。注释可以来自经过策划的数据库如MSigDB12、推断基因调控网络13或任何其他预定义的基因集。这种设计有两个主要优势1潜在变量先验初始化因此可直接解释为生物模块的活性;2生物模块规范的灵活性使其可推广到不同的生物抽象如通路、基因调控网络甚至细胞类型、方法。在初始化此类模型架构时可以建模生物模块与推理网络编码器、生成网络译码器中原始特征之间的依赖关系或两者兼有。我们选择使用深度编码器、稀疏线性解码器结构可以基于提升模型推理能力的目标当近似生物模块的真实后验分布Pzi|习时可以捕捉可解释潜在空间中的更复杂模式而潜在变量与原始基因特征之间的联系则由神经网络中的稀疏线性解码器部分保证补充图1。在可解释的潜在空间上重现生物信息为验证我们提出的模型重现生物通路状态的能力我们将该模型应用于已发表的外周血单核细胞PBMCs数据集该数据集由趋化因子干扰素β14刺激METHODS。我们首先证明模型能够利用反应组通路作为潜在变量捕捉细胞类型和刺激状态图1B。我们进一步研究了单个通路的激活状态具体证明干扰素-α/β信号活性会分离受刺激和初生细胞证实模型能够在其潜在空间捕捉通路活性图1C图1D。我们还进一步研究了参与干扰素诱导免疫细胞激活的其他已知生物通路发现某些细胞过程具有特定细胞类型的激活。例如色氨酸对干扰素的分解反应将先天免疫细胞树突状细胞、FCGR3A单核细胞、CD14单核细胞与适应性免疫细胞NK细胞、T细胞CD8、T细胞CD4、B细胞分离开来如先前研究的那样图1D。这些结果共同展示了模型能够将细胞投射到可解释的空间中从而在细胞过程层面研究特定细胞类型的模式。通常情况下对刺激的预期反应往往并非事先已知。差异基因表达已成为研究样本组间基因表达模式差异的标准工具。该程序也适用于基因表达建模的深度生成网络9。我们提出类似方法基于模型潜伏空间研究生物模块活性差异我们可以提出备选假设并利用模型的变分后验分布对应贝叶斯因子BF方法计算其后验概率的比值。当应用于刺激组与对照组的先天免疫细胞时我们表明本方法发现预期在被刺激组中被激活的通路干扰素信号、色氨酸分解代谢具有显著性|logeBF| 3图1E。我们将差异因子活性方法的结果与标准GSEA工具包的假发现率FDR值进行了比较方法图1F。虽然两种方法均在受刺激组中激活了干扰素-α/β信号通路但GSEA未发现先天免疫细胞中色氨酸分解代谢激活见图1F。总体而言VEGA似乎较少受到GSEA固有的基因集大小偏差问题影响图1F补充图2鉴于数据集的生物学特性VEGA显示出更具上下文特异性且显著的命中。细胞系药物治疗生物反应的大规模研究接着我们研究了模型是否能捕捉并重现癌症细胞系大规模实验中的药物反应模式如近期实验方案如MIX-Seq 17。为此我们收集了97个癌细胞系在5种不同条件下的单细胞数据24小时DMSO治疗对照、24小时曲美替尼治疗MEK抑制剂、24小时达曲替尼治疗突变BRAF抑制剂、24小时纳维托克拉斯治疗Bcl-2抑制剂和24小时BRD3379治疗未知作用机制的工具化合物MoAMETHODS。我们通过结合药物治疗数据集和对照组DMSO数据集初始化潜伏空间训练了每种不同药物治疗的模型共4个模型。总体而言每个模型均能重现细胞系对各自的生物反应治疗图2A补充图3。对于曲美替尼而言G2M检查点活性的重要变化治疗条件的下降与MEK抑制剂19,20的预期MoA相符图2B。接下来我们试图研究是否能重现对照组与处理条件下的生物反应模式。对于每一对我们计算了模型中每个潜在变量通路的贝叶斯因子。所得热图可用于理解和解释所有实验条件下的反应模式图2C。在研究每个数据集的低维嵌入时发现曲美替尼在所有研究药物中产生了最强的转录反应。值得注意的是曲美替尼特异性干扰素α和干扰素γ反应被我们的模型正确重现正如之前的实验21和原始MIX-Seq作者17所示。此外我们发现达布洛芬尼治疗的BRAF突变黑色素瘤细胞系表现出强烈的转录反应与曲美替尼细胞系聚集如MIX-Seq研究中报道的那样图2.C补充图3。总体而言这里呈现的结果与该数据集上之前的基因集分析结果一致展示了我们模型在大规模实验中重现药物反应模式的能力。胶质母细胞瘤的基因调控分析揭示了肿瘤细胞的分层现象如前所述我们模型的一个主要优势是潜在空间规范的灵活性任何生物模块都可以与基因集相关从而编码在我们的模型中。生物学中最重要的方面之一是转录因子对基因表达的严格调控22。分析这些主调控因子的活性对于理解细胞类型或疾病等生物状态非常重要因为它们的活性失调会对基因表达程序和表型产生显著影响23,24。为此我们研究了使用主调控因子指定模型潜伏空间是否有助于理解单细胞胶质母细胞瘤GBM数据集中潜在的基因调控网络GRN25。我们使用了24年报告的GBM ARACNe 13网络来指导模型的结构设计。具体来说转录因子用于初始化模型的潜伏空间转录因子-靶点连接初始化解码器接线。训练模型后我们表明预注释的细胞类型可以在潜伏空间中被重现图2D。我们研究了STAT3和OLIG2的活性这两种分别是间充质MES和前神经PNGBM亚型的主调控因子。我们确认它们在潜伏空间中推断的活性在肿瘤细胞中存在反相关性图2E。此外作为寡突胶质细胞分化的已知主调控因子OLIG2被推断为在寡突胶质细胞前体细胞OPCs中被激活。这些结果表明我们方法在潜在空间和解码器布线的规范上具有灵活性和可扩展性。此外这些结果也证明了我们模型能够忠实地还原生物信息。结合细胞类型和细胞状态表征精炼皮层类器官发育分析现代细胞生物学的一个重大挑战是识别和定义细胞类型和细胞状态以便在共同词汇下系统性地研究稳态和疾病发展。为此我们研究了是否能将细胞类型以标记基因列表定义的标记集形式和细胞状态反应组通路的信息结合在模型的潜在空间中以实现细胞类型和细胞状态的解离对应表示。我们将模型应用于Field等人27页中皮层类器官早期发育期间的细胞数据集包括潜伏空间中研究中定义的主要细胞类型信息图3A。训练后每个细胞类型标记集的活动能够正确区分其对应细胞类型正如原作者所注图3B-D。此外在针对每种细胞类型群体的“一与静止”差分因子分析中对应标记集的活性显著富集|logeBF| 3这可能为研究人员利用细胞类型标记集的潜在活性来注释未知簇提供了机会图3E。为了研究我们的模型是否能将细胞类型身份与细胞状态如分裂细胞群与静止细胞群体区分开来我们将数据集预测为两个组成部分1神经上皮标记集活动一种早期脑前驱和2细胞周期有丝分裂通路活动图3F。如前所述神经上皮标记集活动将神经上皮细胞与数据集其他部分区分开来而细胞周期有丝分裂通路的活性则区分了两个前驱群体桡神经胶质细胞和神经上皮中静止细胞与活跃分裂细胞。为验证被识别为分裂细胞的增殖我们研究了细胞周期有丝分裂通路活动与MKI67基因表达之间的相关性MKI67基因是典型的增殖标志细胞周期有丝分裂通路集中不存在外部验证器见图3G。总体而言MKI67的表达与细胞周期有丝分裂通路的推断活动高度相关R2 0.67。综合来看这些结果展示了我们模型在不同细胞群体中联合推断细胞类型和细胞状态身份的潜力因为结合模型潜在空间中不同的信息来源通路、主调控因子、细胞类型标记可以揭示单细胞身份的不同方面。推理过程推广到样本外数据我们假设模型可用于推断训练时未见的数据的可解释潜在表征我们称之为样本外数据。为确认其可推广性我们在两种情况下评估模型1推断的“生物学泛化”即在Kang等人PBMC数据集训练时遗漏了某一对细胞类型/状态类似于scGen10数据集内设置以及2推理的“技术泛化”即在Kang等人14 PBMC数据集上训练的模型我们称之为研究A在另一个模型上进行评估仅包含Zheng等人对照细胞的PBMC数据集28研究B数据集间设置。对于“生物学推广”我们首先检查了样本外刺激的CD4 T细胞中干扰素-α/β信号通路活性分布与样本内CD4 T细胞推断的活性相符图4A。为了更系统地比较样本外与样本内细胞之间推断的潜伏空间我们采用了贝叶斯微分因子活性过程METHODS分别对1个受刺激的样本内细胞和针对特定细胞类型的对照单元用整个数据集训练的模型和2刺激同一细胞类型的样本外细胞和对照单元训练模型去除一对细胞类型/条件并检查前50个差异激活因子的重叠程度图4B。结果显示样本内与样本外差异激活因子之间存在一致性平均重叠率为72%。为进一步评估数据重建能力我们在样本内和样本外设置中测量了原始数据与解码数据之间的R2见图4C。我们表明在样本外环境中R2仅有微小下降证实了模型能够推广到类似实验环境中产生的未见数据。对于“技术推广”我们再次检查了研究B编码的对照CD4 T细胞的干扰素-α/β通路活性分布与研究A对照的CD4 T细胞一致见图4D。我们还调查了研究A对照细胞在“一对一”差异条件下每种细胞类型的前50个差异因子是否与研究B对照细胞的类似程序重叠图4E。我们显示研究A中平均67%的前50个差异因素与研究B重叠表明模型可以在训练时未见的研究中推广。随后我们研究模型是否能利用推断出的潜空间准确重建两项研究的原始表达特征。我们表明研究B的原始细胞与重建细胞之间的R2值虽低于研究A但在大多数细胞类型中研究A与研究B表达谱的基线相关性有所改善见图4F。我们认为这些结果共同展示了模型能够推断训练时未见数据的可解释潜在潜在空间但部分局限性来自潜在的批处理效应。通过重新训练模型来轻松缓解这个问题因为这个过程相当快。讨论本研究引入了一种新型VAE架构VEGA其解码器接线灵感来自已知生物学用于推断单细胞层面各种生物模块的活性。通过将单细胞转录组数据编码到事先指定的可解释潜空间中我们的方法为分析不同背景下各种生物抽象的活性提供了一种快速高效的方法。我们展示了其在理解特定细胞类型群体对不同扰动反应中的应用提供了对生物模块活动的可解读见解。我们设想我们的模型可用于根据癌症中通路表达优先排序药物因为研究特定细胞群体的反应有助于理解药物的敏感性和耐药性。将药物反应预测模型与我们这样的解释模型整合可能对设计新颖的治疗策略非常有益。此外我们推测我们可解释的VAE将有助于理解一般生物学问题。 潜空间的灵活性为分析生物模块的活性铺平了道路如通路、转录调控因子甚至细胞类型特异性。本研究中展示了模块。我们还进一步证明我们的模型可用于共同探索不同细胞亚群的细胞类型和细胞状态无论是否涉及扰动实验。此外解码器连接的权重提供了潜在变量与原始特征之间关系的直接解释性。例如它可以用于对比推断基因调控网络中的相互作用置信度或通过数据驱动的方式对基因在特定生物模块中的重要性进行排名。 当前架构的明显局限在于模型中稀疏的单层解码器。事实上这种架构设计阻碍了普适性和鲁棒性的进一步提升。因此我们的VAE的生成能力受到了削弱。例如理论上我们的模型可以用潜在向量算术进行可解释的响应预测类似于 scGen10但我们方法有限的生成能力牺牲了预测性能以换取可解释性。因此我们的设计选择在潜在空间的生物学解释性上妥协了这些方面。我们相信网络生物学中的先进见解例如能够更全面描述调控机制的多层GRNs能够缓解这些局限。这将为对特定细胞群体进行靶向的、计算机内激活和抑制生物程序的可能性打开了可能以研究其对发育或疾病进展的影响。 另一方面线性解码器的硬编码连接在上下文需要时没有空间纠正生物因素的先验知识这在其他潜在变量模型如f-scLVM 29中存在。事实上从现有数据库如MSigDB获得的先前生物学知识可能不完整或非上下文特定因为额外的未注释基因在某些生物学因素中可能起重要作用。我们认为通过使用标准正则化方法对解码器的权重进行正则化而非掩蔽这一点可以更容易改善。我们把这个问题留给后续探索。方法VEGA的网络结构我们在此介绍用于生物模块分析的变分自编码器VAE架构。我们的VAE是一个深度生成模型旨在最大化单细胞数据集X在生成过程下出现的可能性过程7,10描述为其中θ是神经网络的可学习参数。这里我们选择让潜在变量Z明确表示生物学抽象如通路、基因调控网络或单元类型标记集。为实现这一点我们将VAE架构中的解码器部分修改为单层、掩码、线性解码器。具体来说该层连接潜在节点zj与输入基因特征之间的连接通过掩码M指定M是一个二元矩阵其中训练过程中我们会归零与掩码权重相关的梯度使反向传播仅适用于来自先验注释的权重其余权重保持为0。此外解码器被限制为正权重w ≥ 0以保持对生物模块活动方向性的可解释性。在模型生成部分明确指定了基因与潜在变量之间的联系后我们确保潜在空间代表生物模块的活动及其可解释性。我们选择将潜在变量Z的集合建模为多元正态分布由我们的推理网络参数化参数可学习φ如下这种变分分布的选择很常见并且在以往的单细胞研究中已被证明效果良好 9,10。训练模型时目标是下界证据ELBO类似于许多标准VAE实现 7,10其中变分分布的期望可以通过对一小批数据的蒙特卡洛积分近似且当我们设先验时库尔布莱克-莱布勒散度项的解是闭形式。在采样变分分布时使用了重参数化技巧7以便在训练模型时应用标准反向传播。 为了保留预先注释的生物网络中不存在的基因信息我们在模型的潜在空间中增加了额外的全连通节点。这有两个效果1它能够模拟未注释基因的表达这对于训练期间数据的良好重建至关重要;2它可以帮助捕捉现有生物因素无法解释的数据额外方差显著提升模型的训练效果。贝叶斯微分因子活动研究两组细胞间因子激活的差异及其意义具有重大意义。为此我们从贝叶斯差异基因中汲取灵感表达式程序在9中引入并提出了类似的差分因子分析方法。对于给定的因子k一对细胞xa xb及其对应的基群IDsa sb例如两种不同的处理条件我们有两个互斥的假设这直观上可以看作是测试一个细胞的生物因子激活平均值是否高于其他细胞期望值代表经验频率。我们通过研究定义为对数贝叶斯因子K来评估最可能的假设与9类似假设单元格独立我们可以计算从各组随机抽样的多个单元格对的平均贝叶斯因子。这有助于我们判断某一因子在某一群体中被激活频率更高。通过本文我们认为如果K的绝对值大于3相当于比值比为≈20因子会显著差异激活。数据集与预处理Kang等人数据集Kang等人14的数据集包含两组PBMCs一组对照组一组用干扰素-β刺激组。我们选择使用与 scGen 作者 10 所描述的相同的预处理步骤。简而言之细胞通过与八个原始细胞类型簇中最大相关性进行注释平均取前20个簇基因的平均值。随后对数据进行了筛选去除表达基因少于500个和表达在5个或更少细胞中的基因。然后对每个单元格的计数进行了归一化和对数变换我们选择了顶部6998个高度变异的基因最终数据集为18,868个细胞。原始数据可在GSE96583获取。Zheng 等人数据集Zheng 等人.28数据集包含来自健康供体的3000个PBMCs。过滤后每个细胞的计数被归一化并进行对数转化。然后我们对基因进行了子集使用Kang等人PBMC数据集中的相同6998个基因。最终数据集包含2623个细胞和6998个基因。原始数据可在 https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.1.0/pbmc3k 年获取。MIX-Seq 数据集MIX-Seq17数据集来自 https://figshare.com/s/139f64b495dea9d88c70我们利用实验3的数据收集了足够的单元实现模型的平滑训练。对于5个可用数据集97个分别用DMSO、曲美替尼、达布拉非尼、纳维托克莱克斯和BRD3379处理的细胞系我们切除了表达基因少于200个的细胞以及表达少于3个细胞的基因。然后我们对每个单元格的计数数进行了归一化并对数变换了数据。最后每个药物治疗实验数据集与对照数据集DMSO处理结合提取出前5000个高度变异基因。最终数据集规模为曲美替尼DMSO数据为16732细胞4999基因;达布拉非尼DMSO数据为16942细胞5000基因;NavitoclaxDMSO数据为14507细胞5000基因;BRD3379DMSO数据为15304细胞5000基因。Darmanis 等人数据集Darmanis等人提供的原始胶质母细胞瘤数据25来自 http://www.gbmseq.org/ 并预处理如下我们移除表达基因少于200个的细胞以及表达在3个或更少细胞中的基因。每个单元格的计数被归一化然后对数变换数据。最后我们将转录组限制为前6999个高度变异基因。最终数据集共计3566个细胞。原始数据可在GSE84465获取。Field 等人数据集Field等人27的皮层类器官数据处理方式与胶质母细胞瘤数据集类似。经过归一化和高度变异基因选择后数据集共有4378个细胞6999个基因。原始数据GSE106245可用。VEGA潜空间基因注释的选择在初始化模型的潜在空间时我们选择使用分子签名数据库MSigDB12中的预注释基因集。特别是我们选择使用标志基因集注释50个基因集或Reactome数据库674个基因集。该研究的MIX-Seq分析部分采用了反应组MIX-Seq分析部分采用了标志性H。对于胶质母细胞瘤细胞的基因调控网络分析我们从胶质母细胞瘤的整体RNA-Seq样本中推导出ARACNe13,30网络。具体来说该网络来自一篇先前发表的论文31并重新用于研究胶质母细胞瘤单细胞转录组学谱。对于皮层类器官分析中的细胞类型标记基因我们联系了作者以获取用于注释这些细胞类型的最相关基因。标记列表作为补充表提供。可视化降维方法在可视化数据集时我们使用了UMAP Python软件包中实现的UMAP算法32。我们使用默认参数除了min_dist参数设置为0.5。我们还使用了 sklearn python 包中实现的 t-SNE 33默认参数。与GSEA的比较我们用Python中的gseapy包中的预秩函数运行了基因集富集分析12GSEA。简而言之我们使用Wilcoxon秩和检验计算了对照组和处理组每个基因的差异表达得分该检验体现在Scanpy包34的 rank_genes_groups 功能中。我们根据基因的测试统计数据对基因进行了排名并以以下设置运行了GSEA预排序最小基因集大小为5最大基因集大小为1000以及1000个排列。我们根据FDR对基因组进行了排名并在FDR≤0.05时考虑了显著的命中。当GSEA返回的FDR为0时我们用1e-5替代以避免对数计算时出现数学错误。评价指标我们计算了轮廓分数以评估模型潜伏空间中的生物分离情况。我们利用潜空间中的欧几里得距离计算每个单元格i的轮廓系数定义为其中ai和bi分别是单元i的平均簇内距离和单元格的最近簇平均距离。我们使用了Kang等人14中的刺激标记或细胞类型标记来评估模型潜伏空间的生物学相关性见补充说明1。sklearn package35 silhouette_score实现用于计算。
