编译 | 刘悦
单细胞多组学数据旨在相同细胞中捕获细胞多分子模式，提供了探索细胞表观基因组和转录组关系的机会。作者开发了MultiVelo，这是一种基因表达的机制模型。它纳入表观基因组数据的同时扩展了RNA速率框架。MultiVelo使用概率隐变量模型从单细胞数据估计染色质可及性和基因表达的转换时间和速率参数，提供表观基因组和转录组变化之间时间关系的定量总结。与仅从RNA估算的方法相比，结合染色质可及性数据显著提高了细胞命运预测的准确性。在来自大脑、皮肤和血细胞的单细胞多组学数据集上拟合MultiVelo揭示了两类不同的基因。该模型还确定了四种细胞状态——表观基因组和转录组耦合的两种状态和两种不同的解耦状态。MultiVelo推断的参数量化了基因占据四种状态的时间长度，根据转录组和表观基因组之间的耦合程度对基因进行排序。最后，作者确定了转录因子表达和结合位点可及性之间以及疾病相关SNP可及性和连锁基因表达之间的时滞。
1
简介
从DNA到RNA再到蛋白质的基因表达调控是控制细胞命运的关键过程。协调、逐步的基因表达变化——基因按一定顺序打开和关闭——是细胞发育过程的基础。高通量单细胞测序技术越来越多地被应用于揭示这些逐步的基因表达变化。由于实验测量会破坏细胞，因此只有时间快照测量可用，这并不能观察同一单个细胞随时间的变化。
计算方法可以利用单细胞快照推断发育过程中连续的基因表达变化。例如，细胞轨迹推断算法使用成对细胞相似性将细胞映射到与预测的发育进程相对应的“伪时间”轴上。然而，基于相似性的轨迹推断无法预测细胞转变的方向或相对速率。推断RNA速率的方法通过拟合微分方程系统来解决这些限制，该系统使用剪接和非剪接转录计数来描述转录变化的方向和速率。最初的RNA速率方法依赖于稳态假设来拟合模型参数。后来的工作开发了一个动态模型，除了稳态外，还显式地拟合了基因表达的诱导和抑制阶段。这种RNA速率的动力学模型还推断出每个细胞的潜在时间值，为重建细胞分化过程中基因表达变化的顺序提供了方法。最近的论文进一步扩展了RNA速率框架，但这些方法没有考虑基因表达。
单细胞多组学测量提供了将表观基因组数据纳入转录机制模型的机会。例如，SNARE-seq、SHARE-seq和10X Genomics Multiome可以量化同一细胞中的RNA和染色质可及性。表观基因组和转录组在细胞分化过程中都会发生变化，因此单细胞多组数据集中的时间快照可能揭示了这些分子层之间的相互作用。例如，如果表观基因组谱系启动发生在特定的基因组位点，单细胞多组学数据可以揭示基因染色质重塑与其转录之间的显著时间间隔。同样，观察转录因子表达的动态变化和假定结合位点的染色质可及性可以揭示它们的时间关系。
现有的RNA速率模型假设基因的转录速率在基因表达的整个诱导阶段是一致的。然而，表观基因组变化在调节基因表达中起着关键作用。例如，从常染色质到异染色质的转变显著降低了该位点的转录速率。因此，更现实的模型将反映增强子和启动子染色质可及性对转录率的影响。
作者提出了MultiVelo，一种从单细胞多组数据集推断基因表达及表观基因组调控的计算方法。作者扩展了动态RNA速率模型，将多组测量纳入其中，以更准确地预测每个细胞的过去和未来状态，联合推断每个模式的瞬时诱导率或抑制率，并确定模式之间的耦合程度或时滞。MultiVelo使用概率潜变量模型估计基因调控的转换时间和速率参数，提供表观基因组和转录组变化之间时间关系的定量总结。
作者证明MultiVelo能够准确地恢复细胞谱系，并量化染色质可及性和基因表达暂时不同步的启动和解耦间隔的长度。作者的微分方程模型精确地拟合了来自胚胎小鼠大脑、胚胎人脑和新生成的来自人类造血干细胞和祖细胞的单细胞多组数据集。此外，作者的模型预测了通过染色质可及性调节基因表达的两种不同机制。最后，作者使用MultiVelo来推断转录因子与其结合位点之间以及GWAS单核苷酸多态性与其连锁基因之间的时间关系。总之，MultiVelo为表观基因组变化在细胞命运转变期间的调节机制提供了基本的见解。
2
结果
2.1.MultiVelo区分了胚胎小鼠大脑中两种基因表达调控模型
作者首先将MultiVelo应用于胚胎小鼠大脑（E18）的10X Multiome数据。MultiVelo准确地拟合了观察到的染色质可及性、未切片的前体mRNA和剪接mRNA在整个脑细胞群体中的计数，识别出426个模式高度符合模型的基因。MultiVelo推断出的速率矢量和潜伏期值准确地恢复了哺乳动物皮层发育的已知轨迹。具体而言，位于室下外侧区（OSVZ）的放射状胶质（RG）细胞产生神经元、星形胶质细胞和少突胶质细胞。在神经元迁移过程中，皮质层以由内而外的方式形成，新生细胞移向上层，老细胞留在深层。RG细胞可以分裂成中间祖细胞（IPCs），作为神经干细胞，并进一步在不同层生成各种成熟的兴奋性神经元
与纯RNA模型（如scVelo）相比，结合染色质可及性和基因表达可提高速率估计的准确性（图1A）。特别是，仅RNA模型预测上层神经元内生物学上不可信的回流（图1B）。MultiVelo和scVelo使用类似的参数设置和估计算法，表明表观基因组数据提供了关于细胞过去和未来状态的重要附加信息，而不仅仅是转录组数据。
图1：MultiVelo揭示了两种不同的基因调控机制。
2.2.MultiVelo鉴定胚胎小鼠大脑中的表观基因启动和解耦
MultiVelo能够量化染色质可及性与分化细胞内基因表达之间的不一致和一致性。具体来说，MultiVelo推断出开关时间参数，这些参数识别每个基因处于四种可能状态之一（启动，耦合，解耦和耦合关闭）。接下来，作者研究了这些推断的状态和时间间隔是否可以准确地捕获胚胎小鼠脑细胞中表观基因组学和转录组学变化之间的相互作用。
MultiVelo 在 10X Multiome 数据中识别出四种状态中每种状态的明确示例（图 2A）。例如，Grin2b是一个仅诱导基因，其表达向神经元命运增加，因此只有诱导状态 - 启动和耦合 - 被预测到该基因（图2A，左）。1型基因Nfix的相位具有完整的轨迹形状，并被标记为所有四种状态（图2A，中间）。相反，Epha5是一个2型基因，它的可及性在整个时间范围内继续上升，因此它只占据耦合的开启和去耦状态（图2A，右）。
状态赋值可以通过在 UMAP 坐标上绘制可访问性（c）和表达式（u和s）并排检查它们来定性地确认（图 2B）。在视觉上，作者观察到当状态赋值耦合打开或耦合关闭时，c和u UMAP 图的颜色匹配，并且当赋值状态启动或解耦时，颜色差异发生。例如，Robo2 RNA表达和染色质可及性之间最大的差异发生在圆圈区域，预计该区域处于解耦状态（图2B，顶部）。Robo2是一个1型基因;染色质关闭开始后，表达保持在相对较高的水平，即使其可及性已经向成熟的神经元下降。同样，Gria2的可及性与解耦状态下的RNA不同（图2B，中间）。2型基因Gria2的染色质可及性在转录诱导阶段之后继续增加。此外，基因Grin2b显示了染色质启动阶段的一个明显例子，在此期间染色质在RNA产生之前打开（图2B，底部）。
图2：MultiVelo量化胚胎小鼠大脑中的表观基因组启动和解耦
2.3.MultiVelo在来自小鼠毛囊的SHARE-seq数据中量化表观基因组启动
最近的一项研究使用SHARE-seq研究毛囊组织中转运扩增细胞（TAC）的快速增殖，这些细胞产生几个成熟的效应细胞，包括内根鞘（IRS）和毛干层：角质层，皮质层和延髓层.当应用于该数据集时，MultiVelo正确识别了从TAC到IRS和毛干细胞的分化方向（图3A）。潜伏时间预测是根细胞单独使用RNA的速率分析未能捕获毛干分化方向（图3B）。与小鼠大脑相比，作者在该数据集中观察到更多的仅诱导基因和更少的Model2基因（图3C).
原始SHARE-seq论文的主要结果之一是鉴定了启动子和增强子染色质可及性预示着基因表达的基因，作者将这种现象称为"染色质电位"。这种现象最明显的例子是Wnt3，它编码一个旁分泌信号分子，在控制头发生长方面很重要。事实上，由可访问性着色的UMAP图以及未剪接和剪接的mRNA表达显示出跨模态的明显时间延迟（图3D）。
图3：MultiVelo量化小鼠皮肤中的表观基因组启动。
2.4.MultiVelo揭示了人类造血干细胞和祖细胞的早期表观基因组和转录组学变化
造血祖细胞由干细胞样细胞群组成，这些细胞群快速连续地分化成各种中间和成熟的血细胞类型，随着它们进入更多的谱系限制状态，自我更新潜力逐渐降低.
作者将纯化的人CD34 细胞培养7天，然后使用10X Multiome平台对其进行测序。作者用单核RNA-seq和ATAC-seq数据获得了11，605个高质量的细胞。使用先前描述的标记基因，作者发现了类似于许多早期血液发育人群的集群，包括 HSC、多效祖细胞 （MPP）、淋巴样引发的多能祖细胞 （LMPP）、粒细胞-巨噬细胞祖细胞 （GMP） 和巨核细胞-红细胞祖细胞 （MEP）。作者还确定了类似于早期粒细胞，红细胞，树突状细胞（DC）和血小板的簇。
血细胞分化是一个具有挑战性的系统，需要用RNA速率进行建模，但作者发现，结合染色质信息可显著提高预测细胞方向的局部一致性和生物学准确性（图4A）。相比之下，仅从RNA推断的速率方法并不能准确反映HSPC的已知分化层次结构。与小鼠大脑一样，MultiVelo预测Model 1基因在此数据集中比Model 2基因更常见;仅诱导是第三大最常见的基因类别（图4B）。观察到的引爆和去耦间隔的中位数长度短于耦合的长度（图4C）。这些模式与作者在小鼠大脑数据集中观察到的一致，这表明可能存在一个共同的潜在生物学机制。
与小鼠大脑数据集一样，HSPC数据集中的Model 2基因对于与细胞周期相关的GO项显着富集。"有丝分裂细胞周期的调节"，"有丝分裂中期/反相过渡的调节"和"有丝分裂姐妹染色单体分离的调节"都富集在FDR <0.002的Model 2基因中。如果作者检查骨髓，红系和血小板谱系的单独轨迹，许多G2 / M相标记基因显示清晰的 Model 2 模式，在表达开始下降后具有最高的染色质可及性（图4D).
图4：MultiVelo鉴定造血干细胞中的启动。
2.5.MultiVelo与发育中人脑中的转录因子、多态性位点和基因表达有关
接下来，作者将MultiVelo应用于最近发布的10X Multiome数据集，该数据集来自人类皮层。与胚胎小鼠大脑数据集一样，MultiVelo推断的速率载体与已知的脑细胞发育模式一致（图5A）。MultiVelo正确地推断出桡神经胶质细胞附近的循环细胞群是潜伏时间最早的细胞类型。相反，在没有染色质信息的情况下推断的速率载体预测了中间祖细胞和上层兴奋性神经元中不协调的回流（图5B).
与小鼠大脑数据集一样，作者确定了Model1和Model2基因的明显例子（图5C），尽管预测在人类数据集中遵循Model2的基因较少（图5D）。有趣的是，仅RNA模型预测的MEF2C是一种Model2基因，具有大部分抑制期。 然而，染色质信息的添加允许正确预测基因具有诱导和抑制阶段。
图5：MultiVelo推断胚胎人脑中的表观基因组和转录组动力学。
3
总结
总之，MultiVelo准确地恢复了细胞谱系，并量化了染色质可及性和基因表达暂时不同步的启动和解耦间隔的长度。作者的模型准确地拟合了来自胚胎小鼠大脑、小鼠背侧皮肤、胚胎人脑和人类造血干细胞的单细胞多组学数据集。此外，作者的模型确定了两类基因，它们在染色质闭合和转录抑制的相对顺序上有所不同，作者在研究的所有组织中发现了这两种机制的明显例子。MultiVelo将提供对一系列生物学环境中基因表达的表观基因组调控的见解，包括正常的细胞分化，重编程和疾病。
作者在PyPI和GitHub上提供MultiVelo的开源Python实(https://github.com/welch-lab/MultiVelo).
参考资料
Single-cell multi-omic velocity infers dynamic and decoupled gene regulation. Chen Li, Maria Virgilio, Kathleen Collins, Joshua D Welch，bioRxiv2021.12.13.472472;
doi: https://doi.org/10.1101/2021.12.13.472472
