编译 | 宋思宁，梁瑶瑶，侯同帅 审稿 | 王海云
本文介绍由美国马萨诸塞州怀特黑德生物医学研究所的Xiaojie Xu和Jonathan S. Weissman以及匹兹堡大学计算与系统生物学系的Jianhua Xing共同发表在Cell的研究成果：基于单细胞测序(scRNA-seq)RNA速度和代谢标记预测细胞状态。作者提出了一个分析框架dynamo (https://github.com/aristoteleo/dynamo-release)，推断绝对RNA速度，重建预测细胞命运的连续向量场，利用微分几何提取潜在的规则，最终预测出最佳的重编程路径和扰动结果。进一步分析了dynamo在克服传统基于剪接的RNA速度分析的基本限制方面的能力，表明其能在代谢标记的人类造血scRNA-seq数据集上精确估计速度。此外，微分几何分析揭示了驱动早期巨核细胞出现的机制，并阐明了PU.1-GATA1电路中的不对称调节。利用最小作用路径方法，dynamo可以准确预测驱动无数造血系统的转变，并最终由计算机干扰预测基因微扰引起细胞命运的转变。综上，Dynamo有助于开展细胞状态转变的定量分析和预测。
1
简介
后生动物具备一个受精卵分化成多种细胞类型的能力，同时可以保持相同的基因组。为了说明这一过程，Waddington引入了表观遗传学，用于直观地解释细胞分化、转分化及重编程。然而，该领域的核心目标是超越这种定性的、隐喻性的概念性解释，设计定量的预测模型。
动力系统建模为深入了解基因调控网络(GRNs)如何控制生物过程提供了一个强大的工具。单细胞基因组学的最新进展使细胞状态转变的分析具有前所未有的分辨率。单细胞分析技术的进步推动了计算方法的发展，通过快照测量推断细胞动力学。其中最主要的是伪时间方法，即通过学习基于转录组相似性的单细胞图流形来推断生命进程的顺序。最近，几个研究小组已经将代谢标记的批量RNA测序技术应用于单细胞方法。尽管单细胞分析取得了惊人的进展，但我们充分利用这些数据的能力受到了限制，因为缺乏一个适当的分析框架来解释数据和指导未来的实验。
在本研究中提出一个构建和解释单细胞转录组向量场的框架。该框架提供了四项创新。首先，通过协调RNA代谢标记和内在剪接动力学，建立了一个包容性的表达动力学模型，不仅准确估计全基因组RNA更替率，而且克服了传统的基于剪接的RNA速度的内在限制，以推断绝对速度。其次，开发了一个通用的算法，从离散的、稀疏的和噪声的单细胞测量稳健地重建连续转录组向量场。第三，将基于机器学习的向量场重构方法的可扩展性与微分几何分析(包括雅可比矩阵、加速度、曲率和散度)的可解释性结合起来，以获得进一步的生物学见解。第四，利用直接从scRNA-seq数据集重建的分析向量场，开发了最小作用路径(LAP)和计算机微扰两种原则方法，以对细胞命运转变的最优路径和关键驱动因素以及遗传扰动的结果进行非微不足道的预测。
2
结果
一个包含向量场和微分几何分析的细胞状态转换的基本框架
原则上，速度矢量场提供了基因如何相互调节的完整描述。作为一个简单的例子，考虑一个双基因开关基序，如参与造血的pu1 /SPI1-GATA1调控网络(图1 A1)。这个主题的向量场通常是制定为一组常微分方程(图1 A1),模型自活化和相互抑制涉及PU.1 GATA1,指定一个细胞的瞬时速度在任何给定的表达状态,并预测细胞的进化状态(图1 A2-4)。我们可以进一步刻画该向量场在其基因表达空间中的拓扑特征，用分隔符将空间划分为三个吸引子盆地，每个吸引子包含一个稳定的不动点(吸引子)，对应于一个稳定的表型(图1 A4)。三个从吸引子A1的同一吸引子盆地的不同状态出发的代表细胞，每个细胞都沿着向量场定义的轨迹(流线)传播，最终定居在同一吸引子状态A1(图1 A2-4)。对向量场的分析也有助于产生关于基因如何调节细胞状态的假设(图1B和1C)。
一些额外的微分几何量提供基因调控的互补信息。加速度场(图1D)显示了基因表达子空间，在这些子空间中，速度发生了显著的变化。因此，早在细胞表现出可识别的谱系特异性基因表达差异之前，就有可能检测到在祖先状态中具有较大加速值的基因，对细胞命运的承诺做出关键贡献。一个相关但不同的量是曲率场，显示了速度方向突然改变的基因表达热点。拨动开关基序说明了矢量场和各种微分几何分析在研究调节网络的动力学中的重要性。然而，这些简化的基序被嵌入到一个未知的全基因组调控网络中(图1A)。因此，应用机器学习方法直接从单细胞测量数据重建转录组向量场是可取的(图1E)。
图1 使用速度矢量场和微分几何分析建模单细胞表达动力学
RNA代谢标记和表达动力学的整合模型提供了mRNA动力学参数的全基因组估计
原始的RNA速率方法使用从cscRNA-seq数据中偶然捕获的内含子读码，并假设普遍的剪接速率常数。假设稳态与极端细胞高表达,并使用替换 (和是各自的拼接和降解速率常数)，认为这种限制可以通过tscRNA-seq来放松，它可以以一种可控的、更少偏差的、时间解决的方式来测量RNA周转动态。
作者构建了一个包含模型(图2A)，建立了从cscRNA-seq和tscRNA-seq数据集中提取RNA动力学信息的统一框架，其中考虑了RNA代谢标记(当使用tscRNA-seq数据时)、RNA剪接和降解。当只有cscRNA-seq数据可用时，或者当需要使用tscRNA-seq实验的剪接数据时，可以使用dynamo来估计相对降解速率常数（）相对剪接RNA速度(图2B)。
为了证明提出的方法的有效性，将dynamo应用于两个先前报道的数据集:通过小鼠ESCs的scNT-seq获得的降解数据集和通过RPE-1细胞的scEU-seq获得的动力学数据集 (图2D-2I)。有趣的是，小鼠基因的剪接和降解率与人类同源基因的剪接和降解率相关，但普遍高于人类同源基因(图2D)，这与之前观察到的情况类似。特别是，在动力学实验中，新的RNA和总RNA呈现出预期的强线性关系，且斜率随标记时间的增加而增加(图2E和2G-2I)。
在动力学实验中，在拼接和未拼接RNA的“相平面”上绘制未拼接/拼接速度，在总RNA和新RNA的“相平面”上绘制新/总速度。新的RNA速率总是非负性，因为在短时间的标记实验中，标记RNA的水平普遍增加(图2G)。
图2. 包含 RNA 代谢标记的单细胞表达动力学模型
用Dynamo基于代谢标记RNA克服了传统基于剪接的RNA速度的基本限制
为了证明大规模的、基于UMI的tscRNA-seq数据集可以改善对cscRNA-seq数据集的速度分析，基于对体外培养第四天和第七天进行多系分化后的人类CD34 HSPCs生成的scNT-seq数据集进行了分析。利用dynamo处理数据，以量化每个细胞中每个基因的未剪接、剪接、新RNA和总RNA，仅基于剪接数据(未剪接和剪接的RNA)进行了cscRNA-seq RNA速度分析(图3A)。结果显示，采用传统的基于剪切的RNA速度方法，得到的RNA速度矢量不仅噪音大，还给出了生物上的错误方向（从成熟细胞类型如巨核细胞等指向祖细胞如祖细胞样细胞等），而相比之下，在dynamo的框架下，基于代谢标记的RNA速度则能正确的指向干细胞分化的方向。为了评估dynamo反卷积正交细胞过程的能力，重新分析了细胞-细胞周期进展和糖皮质激素受体(GR)激活探究的sci-fate数据集，对原始研究中检测到的GR反应和细胞周期基因的组合或个体集进行了时间解析的总RNA速度分析。结果显示，在DEX(地塞米松)初始处理后，细胞在0到2,4,6和8 h的时间点发生了平稳的顺序转移，并从细胞周期基因组的分析中确定了一个与细胞周期进程匹配的循环(图3F)。又进一步分析了来自scEU-seq研究的数据集，基于标记的RNA速度分析准确地显示，在scNT-seq mESC数据集上Tet 1/2/3三重敲除(TetTKO)下，罕见的2C样全能型细胞增加。（图3H）
图3.代谢标记实验改进和推广了RNA速度估计
精确、稳健、高效的单个细胞向量场重建
为计算速度测量的噪声和离群值，sparseVFC依靠EM算法迭代优化内层集合，以及每个基函数的优化系数集，进一步提高了向量场重构的鲁棒性。（图4A）
为了探索向量场重构的潜力，首先测试了5,000个随机采样点的模拟数据集上重构的有效性。结果显示，重构向量场的流线估计以及不动点等与解析流线几乎无法区分，此外，还可以精确地恢复通过状态空间的雅可比矩阵，估计的高阶向量微积分量与真实的解析计算量非常匹配，所推导的向量演算的解析公式速度远远快于最先进的同样需要重构向量场的数值方法（图4B.C.D.E）。还在其他数据集上验证了重建的向量场可以几天内的预测单细胞的轨迹。
图4. 绘制矢量场，量化其地形，并转向微分几何分析
微分几何分析揭示了造血的时间和调节机制
接下来设计了一套与矢量场相关的微分几何分析来揭示基因调控的定量信息(图5A)。首先使用tscRNA-seq数据集学习了向量场。UMAP基于空间的向量场中确定的不动点正确地反映了系统的拓扑结构。然后将向量场组织成树状结构，该树状结构正确地总结了造血系统的谱系层次(图5B.5C)。
首先基于向量场的伪时间确认了巨核细胞谱系的早期出现，又发现巨核细胞谱系主调节因子FLI1在造血干细胞状态下就已经表现出高表达。此外，雅可比矩阵分析显示，FLI1和KLF1之间存在相互抑制和FLI1的自激活。这些分析共同表明，FLI1的自激活在造血干细胞状态下能保持其较高的表达，这使得HSPC首先以高速和加速的方式向Meg谱系转移，同时通过抑制KLF1抑制其向红细胞的转移。雅可比矩阵分析显示，CEBPA抑制RUNX1和GATA2，以及RUNX1和GATA2之间的相互激活和它们在祖细胞中的自激活。（图5G）
为了从不同的角度了解关键调控的机制，发展了三个互补的策略:细胞层面、轨迹层面和平面层面的分析。用典型的SPI1(PU.1)-GATA1展示了这些策略。早期研究揭示他们相互拮抗又各自自我激活，而这些调控关系都包含在由测量数据重构的向量场中。对向量场的分析得出GATA-SPI1非对称的调控机制，SPI1对GATA1转录的有效抑制需要达到一定的阈值，但GATA1对SPI1的抑制为低阈值。在功能上， GATA1-SPI1的不对称可能有助于平衡谱系发育。（图5H.I）
图5. 人类造血过程的矢量场和微分几何分析
最小作用路径预测最佳造血细胞生命转换的驱动因素
借助从可用的 scRNA-seq 数据集构建的连续向量场，开发了一种原则性策略，来揭示最佳路径、与驱动相关的关键转录因子(TFs)以及与它们相应的表达动态（图6A）。造血scNT-seq数据集非常适合测试具有许多已知发育、去分化和转分化事件的这种方法（图6B）。
首先，采用 LAP (F-LAP) 来描述每个转换过程。从 HSC 到终端细胞类型的发育 F-LAP 不仅仅是基因表达空间中的最短路径，而是遵循包含表达动力学信息的向量场指定的弯曲流动(图6D)。此外，发育F-LAP 与去分化 LAP 路径不同，并且通常具有更短的过渡时间和更小的动作(图6F)。类似地，从一种细胞类型到另一种细胞类型的转分化 LAP不同于反向转变，这表明细胞是一个不可逆系统这一事实（图S7B）。值得注意的是，HSC 向 Meg 谱系分化的发育 F-LAP 需要最少的时间（约 31 小时），这与Yamamoto 等人在2013年中报道的结果一致，进一步证实了早期的观察Meg谱系的外观（图6E和S7D）。
接着，测试了 LAP 方法优先考虑各种造血细胞命运转变的关键驱动因素的能力(如图6H、6I、S7E和S7F)。对所有报告的正常发育和重编程实验编译了已知的TFs，并根据它们的累积 MSD 对其进行评分。实验确定的所有报告的转分化事件的TFs在LAP分析中排名一样高(大多在前80%)(图6H和图S7F)，所有报告的转移的AUC曲线下总面积得分约为0.85(图6I)。
以上分析揭示了 LAP 方法在预测细胞-命运转变的最佳路径和TF混合物方面具有高准确度的潜力。
图 6. 最小作用路径方法准确预测最佳细胞转化路径
通过计算机干扰实验预测遗传扰动后的细胞命运转移
通过计算机干扰预测每个细胞中每个基因的表达响应和遗传扰动后的细胞命运转移(图7A)。结果显示抑制 GMP 谱系的主调节因子SPI1，将细胞转向巨核细胞和红细胞，而抑制 MEP 谱系的主调节因子 GATA1，将细胞转向单核细胞和中性粒细胞(图7B)。然而，抑制这两个基因会使细胞陷入祖细胞状态。这些预测与Rekhtman et 等人在1999年中发表的预测非常吻合，并揭示了SPI1和GATA1在驱动GMP和MEP谱系中的跷跷板效应调节(图7B, iii)。类似地，抑制HSPC维持基因HLF1会促使细胞远离祖细胞(图7B, v)，最后，GATA1、KLF1和TAL1、已知的红细胞因子，以及用于将成纤维细胞重新编码为红细胞的 TF 的三重激活将大多数其他细胞转移到 Ery 谱系中（图 7B，vi）。
图 7. 硅扰动剖析了遗传扰动下的细胞命运转变
3
讨论
单细胞基因组学的实验进展为了解单个受精卵如何以精确协调的方式产生复杂的有机体这一过程提供了独特而丰富的视角，但缺乏适当的分析框架来利用这些数据。在这项研究中，开发了 dynamo，通过将黑盒机器学习方法与可解释的动态系统方法相结合，从单细胞数据集中获得定量见解，从而填补这一未得到满足的空白。
本研究的分析框架包括四个完整的阶段。首先，从单细胞数据估计全基因组动力学速率常数和 RNA 速度向量。接下来，使用 RNA 丰度和速度向量来重建向量场。然后，应用由解析矢量场实现的微分几何分析，从而获得生物学见解。最后，应用 LAP 方法和计算机干扰实验来预测细胞状态转换的最佳路径和遗传扰动的结果。在动力学参数估计的第一阶段，由于dynamo实现了一个通用的建模系统，它与现有的 RNA 代谢标记策略以及可能开发的新标记协议广泛兼容。此外，收集了一些关于人类造血高质量的tscRNA-seq数据集，并建立了从该数据集和其他使用dynamo的tscRNA-seq数据集估计的总RNA速度，克服了传统RNA速度估计的固有限制。
在第二阶段，将单细胞速度向量样本作为输入来稳健地学习转录组空间中的连续向量场。在伪时间排序、RNA 速度和 sci-fate方面的早期努力下共同推动了动力学推理的重要发展。在这方面关键性发展是现在能够在转录组空间中重建分析和连续向量场。通过重建的连续向量场，可以预测过去或未来较长时间段内的细胞状态，正如对中性粒细胞分化或小鼠造血的顺序转录组分析和克隆命运追踪的分析所证明的那样。Dynamo还能够在计算机上追踪细胞群随时间推移的转录组动力学，这可能为活细胞成像或谱系追踪提供重要补充。
在第三阶段，应用预测动力系统方法和微分几何分析，从矢量场中提取调控信息。因此，Dynamo 可以使用单细胞基因组学数据直接探索调控机制，甚至恢复动力学参数，例如希尔系数、潜在的细胞命运转变。
在第四阶段，提出了两种原则性方法，即最小动作路径(LAPs)和计算机扰动(in silico perturbation)，分别预测最优转变路径和遗传扰动的结果。进行计算机扰动实验有助于通过大量可能的成对和高阶扰动进行搜索，从而发现令人们感兴趣的细胞状态和转变的基因组合。
总之，本文建立了一个分析转录动力学的通用框架，可以应用于许多生物系统。若再结合单细胞方法的实验进展，包括 RNA 代谢标记、谱系追踪、RNA年龄、信号通路记录以及遗传扰动，dynamo使我们转向整体动力学模型和用于细胞图谱项目的整个生物体的理论，以了解复杂的细胞状态如何由有限数量的因素的组合规则产生，并最终解决在任何细胞类型之间转换的最终目标。
代码
Dynamo包下载：
https://github.com/aristoteleo/dynamo-release
论文图片使用教程：
https://github.com/aristoteleo/dynamo-notebook
Dynamo使用教程：
https://github.com/aristoteleo/dynamotutorials
参考资料
Qiu X., Zhang Y., Martin-Rufino J.D. et al. Mapping transcriptomic vector fields of single cells. Cell (2022).
doi: 10.1016/j.cell.2021.12.045.
