RegVelo: 基因调控启发的单细胞动态.

• paper：RegVelo: Gene-regulatory-informed dynamics of single cells

0. TL; DR

细胞命运的转变是由复杂的基因调控网络（GRN）驱动的，然而，经典的RNA velocity模型在推断细胞动态时并未显式地考虑基因间的调控关系，这限制了我们对底层机理的理解。反之，传统的GRN推断方法大多忽略了生物系统的动态特性。为了解决这一概念上的脱节，本文提出了RegVelo，一个自下而上、可操作且可解释的深度学习框架，它能够同时对RNA剪接动力学和基因调控相互作用进行建模。

在多种生物系统中，RegVelo都展示了其在预测细胞终末状态、基因相互作用和扰动模拟方面的可靠能力。通过将RegVelo应用于斑马鱼神经嵴发育的全长Smart-seq3数据以及配对的基因表达和染色质可及性数据，作者描绘了控制命运决定的调控程序。在in silico（计算机模拟）扰动预测的指导下，并通过CRISPR-Cas9敲除和单细胞Perturb-seq实验的验证，作者确定了tfec是色素细胞命运的早期驱动因子，并揭示了elf1是一个新的调控因子。RegVelo为连接基因调控和细胞命运决定建立了一个定量的框架，使得研究人员能够模拟调控扰动并生成可检验的假设。

1. 背景介绍

单细胞技术以高分辨率揭示了细胞分化的过程，而Waddington的能量景观理论将表型状态形象地比喻为“山谷”，分化路径则如同山脊。然而，构建一个能同时捕捉GRN和细胞在这一景观上动态演化的模型，一直是计算单细胞生物学中的一个重大挑战。这类模型对于洞察生物学机理、预测基因敲除的动态效应以及指导实验至关重要。

目前主流的计算方法存在明显的“二分法”：

• 轨迹推断方法：如伪时间（pseudotime）和RNA velocity，它们专注于重构细胞的分化路径。伪时间方法能排列细胞顺序，但缺乏模拟动态的动力学信息。RNA velocity通过对mRNA剪接过程建模，能够估计细胞状态的瞬时变化，但它依赖于严格的假设，如基因间相互独立和转录速率恒定，从而忽略了细胞分化中至关重要的转录调控。
• GRN推断方法：这类方法利用单细胞RNA测序（scRNA-seq）或多组学数据（如scATAC-seq）来推断基因间的调控关系。这些推断出的GRN常被用于量化基因重要性或预测基因扰动的影响。然而，它们通常不直接建模细胞的动态过程，或者只是在事后将调控与动态联系起来。

因此，当前迫切需要一种能直接将基因调控与细胞动态联系起来的模型，以揭示GRN如何编排动态的细胞过程。

为了填补这一空白，作者提出了RegVelo（regulatory velocity）。RegVelo是一个深度生成模型，它将剪接动力学与基因调控耦合起来，通过一个受先验GRN指导的神经网络，将转录过程建模为一个依赖于调控因子和时间的动态过程。这使得RegVelo能够构建一个非线性的、全基因组范围的动态微分方程系统，并提供一个连续的速度向量场。更重要的是，它能够进行in silico的调控扰动模拟，从而在机理层面预测调控变化如何驱动细胞命运的决定。

2. RegVelo 方法

RegVelo是一个贝叶斯深度生成模型，它以未剪接（unspliced）和已剪接（spliced）的RNA丰度为输入，旨在推断基因和细胞特异性的潜变量，包括潜伏时间（latent time）、后验速度分布以及代表转录调控的GRN。

2.1 耦合基因调控的动力学方程

与传统的RNA velocity模型为每个基因建立独立的动力学方程不同，RegVelo通过一个高维常微分方程（ODE）系统来描述所有基因的耦合动态。对于每个基因$g$，其未剪接RNA $u_g(t)$和已剪接RNA $s_g(t)$的动态由以下方程组描述：

这里的关键创新在于转录速率$\alpha_g$不再是一个常数，而是被建模为一个受基因调控影响的、随时间变化的函数：

其中$\mathbf{s}(t)$是所有基因在时间$t$的已剪接RNA丰度向量。$\mathbf{W}$是一个GRN权重矩阵，其中$W_{gj}$代表基因$j$对基因$g$的调控效应（正为激活，负为抑制）。$\mathbf{b}$是基础转录速率向量。$h$是一个非线性激活函数（如softplus），确保转录速率为非负。

通过这种方式，RegVelo将所有基因的动力学通过GRN矩阵$\mathbf{W}$耦合在了一起。

2.2 模型架构与推断过程

RegVelo的整体架构是一个深度生成模型，其推断过程基于变分推断：编码器神经网络将每个细胞的$u$和$s$丰度编码为一个低维的潜变量$z^{(n)}$。潜伏时间解码器（Latent Time Decoder）以潜变量$z^{(n)}$为输入，为每个细胞$n$和每个基因$g$输出一个特异性的潜伏时间$t_{ng}$。

转录速率预测器（Transcription Predictor）参数化了GRN权重矩阵$\mathbf{W}$。它接收所有基因的$s$丰度作为输入，并根据公式计算出每个细胞每个基因的转录速率$\alpha_{gn}$。这个网络在训练过程中会受到一个先验GRN的指导和约束。

结合学习到的转录速率$\alpha_{gn}$、剪接速率$\beta_g$和降解速率$\gamma_g$，模型使用一个可并行的ODE求解器来积分上述动力学方程，从而预测在潜伏时间$t_{ng}$下的$u$和$s$的丰度。

模型假设观测到的$u$和$s$丰度服从以预测值为均值的高斯分布。通过最大化证据下界（ELBO），模型可以同时优化所有网络参数和动力学参数（$\mathbf{W}, \mathbf{b}, \beta, \gamma$等）。

2.3 对GRN的建模与约束

RegVelo将GRN表示为一个权重矩阵$\mathbf{W}$，并通过两种策略来结合先验知识并处理其稀疏性：

• 硬约束（Hard constraints）：只允许先验GRN中存在的边拥有非零权重，其他边的权重被强制设为0。
• 软约束（Soft constraints）：允许模型学习先验GRN中不存在的新连接，但通过一个正则化项$R_{prior}(\mathbf{W}; \mathbf{G}) = |\mathbf{W} \odot (1 - \mathbf{G})|_2$来惩罚这些新连接，使其倾向于被学习为零。这为校正不完整的先验GRN提供了灵活性。

通过对ODE系统的雅可比矩阵（Jacobian matrix）进行L1正则化，可以促进学习到的GRN的稀疏性，这符合生物系统中调控网络的普遍特征。

2.4 下游分析与扰动模拟

训练完成后，RegVelo能够提供丰富的信息：

• 后验速度分布：作为一个贝叶斯模型，RegVelo可以量化速度估计的不确定性（uncertainty）。
• In silico扰动模拟：这是RegVelo的核心功能。通过在训练好的模型中修改GRN矩阵$\mathbf{W}$（例如，将某个转录因子对应的列全部设为零，模拟该TF的敲除），模型可以重新计算速度场，并预测扰动对细胞动态和最终命运的影响。作者通过计算原始速度场和扰动后速度场之间的余弦相异度来量化局部扰动效应，并通过CellRank框架来评估扰动对全局细胞命运概率的影响。

3. 实验分析

作者通过一系列模拟和真实的生物学数据集，系统地验证了RegVelo的性能和应用价值。

实验一：在模拟和细胞周期数据集上的基准测试

在一个具有明确方向的U2OS细胞周期数据上，RegVelo的推断结果与实验观测高度一致。

• 图C展示了基于FUCCI荧光标记的细胞周期真实进程。
• 图D（左）显示RegVelo推断的潜伏时间与细胞周期进程正相关（Spearman相关系数=0.683），速度流（velocity stream）也准确地指向了细胞周期的前进方向。其跨界正确性（cross-boundary correctness, CBC）得分高达0.864。
• 图D（右）展示了推断出的GRN，其中识别出的核心调控因子（如TGIF1）及其靶基因（如BUB1, TOP2A）均是已知的细胞周期相关基因。

实验二：胰腺内分泌细胞生成的调控动态解析

作者将RegVelo应用于小鼠胰腺内分泌细胞生成的数据集，并结合CellRank框架进行下游分析。

• 图B：RegVelo的速度场与CellRank结合，成功识别了四个已知的终末细胞状态（alpha, beta, delta, epsilon）。
• 图C：在识别终末状态的任务中，RegVelo能够正确识别所有终末状态，其性能（TSI score）优于不考虑GRN的scVelo和veloVI。
• 图D, E (In silico扰动)：作者模拟了Neurod2和Rfx6这两个转录因子的调控网络敲除。
  • 图D显示，扰动Neurod2的调控网络主要影响epsilon细胞的命运，其中Neurod2-Rfx6这条调控轴的贡献最大，这与近期的实验研究相符。
  • 图E显示，扰动Rfx6的调控网络则同时影响beta和epsilon细胞的命运，这反映了Rfx6在内分泌谱系决定中的多效性。

实验三：人类造血分化和肌生成中的调控基序恢复

• 人类造血：在传统RNA velocity方法表现不佳的人类造血数据集上，RegVelo成功识别了所有五个终末谱系，并准确重构了分化层级。
  • 图B（右）显示，在所有被测试的方法中，只有RegVelo能够正确识别所有五个终末状态（TSI score ≈ 0.95）。
  • 图C显示，模型推断的不确定性（uncertainty）在祖细胞中较高，在终末状态细胞中较低，这与祖细胞具有更高可塑性的生物学事实相符。

• 造血驱动因子：
  • 图D（中）显示，RegVelo成功恢复了著名的GATA1-SPI1“拨动开关”调控基序。模拟敲除GATA1或SPI1的调控网络，能够正确预测其对红细胞和单核细胞谱系的相反影响。
  • 图E, F显示，在预测谱系驱动基因的任务中，RegVelo（无论是基于扰动模拟还是CellRank）的性能（AUROC）均显著优于其他方法，甚至包括使用了代谢标记数据的dynamo。

• 人类肌生成(G-L)：作者展示了如何利用旁路信息（如已知的发育阶段）来指导RegVelo的模型选择。
  • 图H显示，通过比较潜伏时间与发育阶段的相关性，作者选择了一个最优的RegVelo模型（Soft-r）。
  • 图I显示，这个选定的模型在预测细胞状态转换（CBC score）上显著优于scVelo和veloVI。
  • 图J, K, L显示，该模型成功识别了已知的肌生成调控因子（如MSC和MYOG），并且模拟敲除MSC调控网络后，能够正确预测下游成熟标志物（如MYOG）的上调。
• 人类后脑发育 (M-P)：在人类后脑发育数据中，RegVelo同样成功恢复了已知的从放射状胶质细胞到不同神经元类型的分化层级。

实验四：斑马鱼神经嵴发育中的调控机制发现与验证

这是本文最核心的应用，展示了RegVelo如何从头发现新的生物学机制并指导实验验证。

首先，作者利用RegVelo分析了自采的斑马鱼神经嵴发育Smart-seq3数据，并结合了从多组学数据推断的先验GRN。

• 图C, D：模型成功恢复了神经嵴分化的终末状态和正确的发育层级。
• 图E, F：在预测谱系驱动因子的任务中，RegVelo的性能（AUROC）优于其他基于相关性的方法。特别是，它能够识别出那些在分化早期表达、但在终末状态下调的“早期驱动因子”（如ets1, sox9b），而传统方法则会忽略它们。
• 图G-J (实验验证)：作者通过CRISPR-Cas9敲除和HCR-FISH实验，验证了RegVelo预测的几个间充质谱系驱动因子（sox9b, sox6, zic2b, hoxa2b）的功能。结果显示，敲除这些基因确实导致了相应谱系的细胞数量显著减少，证实了模型的预测。

为了系统性地验证in silico扰动预测，作者进行了大规模的Perturb-seq实验，并开发了新的评估指标。

• 图B-G：以nr2f2为例，MELD分析和RegVelo的密度变化似然度（density change likelihood）都准确预测了其敲除主要影响面部间充质细胞，而对咽弓2细胞影响不大。这一谱系特异性效应也通过HCR-FISH得到了验证（图H）。
• 图I-K：对另一个驱动因子nr2f5的分析也得出了类似的一致性结论。
• 图L (系统性比较)：在包含11个单基因和4个多基因敲除的Perturb-seq面板上，RegVelo的预测与实验结果的相关性（Spearman correlation ≈ 0.51）远高于dynamo和CellOracle（< 0.32），精确率和召回率也接近后者的两倍。

新生物学发现有力地证明，RegVelo不仅是一个强大的动态推断工具，更是一个能够生成可检验假设、并指导实验发现新生物学机制的in silico细胞：

• 图A：RegVelo将tfec预测为色素细胞谱系的一个新的早期驱动因子。
• 图B, C：模型预测的tfec下游靶基因富集在黑色素细胞分化通路。
• 图D, E：实验验证显示，敲除tfec会导致其靶基因slc45a2阳性的色素祖细胞数量显著减少。
• 图F-H：进一步的分析和实验表明，tfec与其他色素谱系因子（如mitfa）存在功能冗余，同时敲除tfec, mitfa, tfeb会导致严重的色素缺失表型。
• 图I：ChIP-seq实验进一步揭示，tfec可能是由神经嵴关键调控因子Sox10直接调控的。

RegVelo还预测了另一个ETS家族转录因子elf1在色素细胞谱系中的促进作用。这一预测同样通过Perturb-seq和HCR实验得到了证实，并且作者通过CUT&RUN实验识别了elf1的直接靶基因，从而揭示了其完整的调控上下文。