Basic Information
英文标题： Statistical inference with a manifold-constrained RNA velocity model uncovers cell cycle speed modulations
中文标题：使用流形约束的RNA速度模型进行统计推断揭示了细胞周期速度的调制
文章作者：Alex R. Lederer | Gioele La Manno
文章链接：https://www.nature.com/articles/s41592-024-02471-8
Abstract
Para_01
在生物系统中，细胞经历协调的基因表达变化，导致转录组动态在低维流形内展开。
虽然可以使用 RNA 速度提取低维动态，但这些算法可能脆弱，并依赖于缺乏统计控制的启发式方法。
为了解决这些挑战，我们引入了一种 RNA 速度的贝叶斯模型，该模型在一个重新制定的统一框架中耦合了速度场和流形估计，识别出显式动力系统的参数。
专注于细胞周期，我们实现了 VeloCycle 来研究一维周期流形上的基因调控动态，并通过活体成像验证其推断细胞周期的能力。
我们还应用 VeloCycle 揭示区域定义的祖细胞和 Perturb-seq 基因敲除中的速度差异。
总体而言，VeloCycle 通过一个模块化且统计上一致的 RNA 速度推断框架扩展了单细胞 RNA 测序分析工具箱。
Main
Para_01
单细胞 RNA 测序（scRNA-seq）以破坏性的方式捕获基因表达的静态快照，这使得难以解释生物过程的动力学方面。
为了解决这个问题，已经出现了从 scRNA-seq 数据中重建细胞状态之间时间信息的计算方法。
例如，RNA 速度利用未剪接和已剪接转录本之间的比例来估计一个描述基因表达变化率的向量。
该模型考虑了一组描述 mRNA 生命周期的一阶常微分方程，其关键参数是剪接和降解速率。
Para_02
原始的 RNA 速度框架在 velocyto 中实现，固定了一个跨基因的共同剪接速率，从剪接-未剪接相图中推断出相对的基因依赖降解率。
RNA 速度估计的扩展模型是‘动力学模型’，首次在工具 scvelo 中实现，该模型为每个基因引入了细胞层面的潜在时间，以支持沿伪时间轴变化的动力学参数的估计，使其可以直接识别。
通过利用期望最大化算法，scvelo 估计潜在时间和动力学参数。
其他方法已经利用了这些建模思想或致力于扩展它们；然而，RNA 速度分析仍然对预处理选择高度敏感，并且需要各种启发式方法来获得最终估计。
Para_03
一种普遍但可能危险的启发式方法是用于近似 RNA 计数期望的最近邻平滑；这一过程可能会使信息从某些基因泄露到其他基因，并导致失真。
此外，使用通用非线性降维技术将高维速度向量映射到二维嵌入（例如，均匀流形逼近和投影（UMAP）和 t-分布随机邻居嵌入（t-SNE））可能会引入伪影。
例如，与正交过程（如增殖和分化）相关的速度可能会混合在一起，而相邻但不相关的细胞群可能会影响最终的向量。
其他通常需要关注的算法步骤和特殊情况也已经被指出；单细胞代谢标记测量方法可以解决其中的一些问题，但它们的应用仅限于特定的实验设计和体外环境。
Para_04
然而，一些最早期且至今仍常用的一些 RNA 速度模型的一个很少被讨论的局限性在于它们依赖于基因特异性的动力学参数和速度。
在这种情况下，即使事后寻求全局协调，估计的动力学参数仍然保持独立；这导致了一个物理上和几何上不一致的速度向量，该向量的基因特异性成分处于不同的时间尺度，其最终方向不一定与细胞穿过的低维流形相切。
因此，希望进行联合基因拟合以规范估计值，这是一种由最近的方法引入的策略，其中流形、非恒定的动力学参数和速度都是一个非线性函数（例如，神经网络）的输出，该函数具有共享的潜在表示。
然而，这种非结构化的相互依赖并不能完全控制从数据到估计的信息流动，并使得理解正则化是如何应用变得困难。
Para_05
最后，RNA速度缺乏已建立的真实情况，限制了对新开发方法进行敏感性分析的严格性，使得高级扩展的基准测试环境变得具有挑战性。
特别是，过度参数化成为一个问题，尤其是对于假设较少、包含多个非线性或自由度较高的模型。
此外，提出的‘动态模型’的贝叶斯公式返回一个高维均值场后验分布，这与低秩动力学的假设不一致，并且不适合于速度推断和细胞群体动力学的统计比较。
Para_06
我们通过将 RNA 速度分析重新表述为基于流形约束的概率模型的推理框架来应对这些挑战。
采用这种方法，我们提出了 RNA 速度作为定义在流形坐标上的场的显式参数化。
我们专注于一个称为 VeloCycle 的框架中的一个维度（1D）周期性流形，这使得模型验证和应用于细胞周期动力学成为可能。
细胞周期是生物学中最普遍的周期性过程，对胚胎发育、组织再生和疾病起着根本性的作用。
尽管在单细胞 RNA 测序（scRNA-seq）数据集中普遍存在，但默认的细胞周期分析管道仍仅限于基于少数标记基因的分类阶段分配。
在这项工作中，我们不仅解决了在速度估计过程中保持几何约束的更广泛问题，还在改进 scRNA-seq 数据中的细胞周期分析方面取得了进展，突出了其连续性质，并提供了对实际生物时间尺度的控制。
我们将 VeloCycle 应用于不同的生物学背景，通过延时显微镜测量进行实验基准测试，并展示了执行统计检验的能力。
Results
Manifold-constrained RNA velocity addresses shortcomings
流形约束的RNA速度解决了现有方法的不足
Para_01
我们首先试图通过将流形和速度推断统一到一个单一的概率框架中来重新设计RNA速度估计（图1a，左）。
在我们的模型中，细胞作为所有测量基因空间内嵌入的低维流形x上的点随时间移动。
剪接和未剪接分子仅作为x的函数（s(x)和u(x））。
然后，通过将速度向量场参数化为流形坐标的自主函数V(x)，我们将RNA速度向量限制为与流形相切（图1a，右）。
这与先前的方法相反，在先前的方法中，速度方向不受限制，是基于每个基因的估计结果（图1b）。
我们对预期剪接计数求导，应用链式法则，并代入动力学方程，以获得一个将所有基因的动力学参数和潜在坐标动态联系起来的速度向量场（方法）。
测量的原始读数中的噪声被建模为负二项式，也作为流形的函数。
对于所有其他参数，包括每个基因的剪接（β）和降解（γ）速率，选择生物化学知情的先验（图1c和方法）。
Fig. 1: Statistical inference of RNA velocity with a manifold-constrained framework for the cell cycle.
a, 联合框架的示意图，用于参数化基因表达流形和RNA速度场。
c, 潜在变量和可观测数据之间概率关系的板图。S 从期望、流形坐标和流形几何中采样。U 从流形信息、动力学参数和速度函数中采样。坐标定义每个细胞在潜在空间中的位置，而几何定义沿流形的表达变化。
d, 使用细胞特异性坐标（x）和基因特异性几何家族（f），对剪接计数（s）定义流形公式，通过该公式可以将观测数据直接映射到高维空间（顶部）。底部：未剪接计数（u）的速度公式被定义为速度场函数（V），具有连锁的动力学参数（β, γ）。我们通过对这些实体应用链式法则来获得速度估计，将速度描述为流形 x(t) 的直接函数。
e, 周期性过程的流形约束速度估计示意图。首先，流形学习估计坐标和几何；其次，速度学习估计动力学参数和速度函数。
f, 使用 VeloCycle 可能进行的新类型速度分析的示意图：(i) 多个样本之间的统计可信度测试以及与零假设的对比；(ii) 通过 MCMC 采样对模型参数的后边缘分布分析；(iii) 速度外推至实际生物时间，可通过活体显微镜验证；(iv) 从大型参考集到小型目标数据集的基因流形转移学习。星号表示统计显著性。NS，不显著。
Para_02
这种公式构成了一种潜在变量框架，用于估计基因表达流形和RNA速度。
特定维度、拓扑结构及相关的功能参数化约束其几何形状的选择可以根据具体应用进行调整（图1d；讨论）。
我们提出了两种统计学习程序中的推断方法：(1) 流形学习，联合学习定义基因表达空间几何的参数，并为每个细胞分配一个流形（潜在）坐标；以及 (2) 速度学习，基于流形几何和细胞坐标找到速度场和动力学参数（图1d，e）。
Para_03
我们在这个方案中考虑了一种流形拓扑先验信息很强的情况：细胞周期，一个基因表达平滑变化且可以使用傅里叶级数参数化的1D周期空间。
我们的框架 VeloCycle 构成了一个具有两组潜在变量的生成概率模型，并在 Pyro31 中求解（方法）。
第一组与流形学习有关，定义了低维流形 x，该流形被参数化为细胞周期阶段（φ）和基因特异性傅里叶系数（ν0, ν1sin 和 ν1cos），使用预期的剪接计数作为阶段的函数（图 1e 和扩展数据图 1a,b）。
第二组与从预期未剪接计数中学习速度有关，包括基因特异性降解率（γg）、有效剪接率（βg）和速度谐波系数（νω），这些参数化了一个描述细胞周期速度沿流形（φ）变化的角速度函数 ω(φ)（图 1e 和扩展数据图 1c,d；方法）。
使用随机变分推断（SVI），VeloCycle 返回潜在变量的联合后验概率，可用于（1）进行统计速度显著性测试；（2）表征潜在变量不确定性之间的基础相关性；（3）估计生物相关时间尺度上的细胞周期速度；以及（4）通过迁移学习促进速度在小数据集中的应用（图 1f）。
Sensitivity analysis on simulated data validates VeloCycle
对模拟数据的敏感性分析验证了VeloCycle
Para_01
设计完我们的模型后，我们试图通过模拟数据来评估其性能，因为没有真实的数据集包含有关阶段、速度和 RNA 动力学参数的真实信息。
我们采用了一种旨在保留真实数据中预期的重要关系的模拟方法，并避免生物学上不可能的情景（方法和扩展数据图 2a-c）。
具体来说，我们在剪接速率、降解速率和基线（平均）表达水平之间引入了正相关性（r = 0.30），这有助于确保生物学上的合理性，并避免与经验观察相悖的不现实参数配置（扩展数据图 2a）。
这种结构自然地导致了剪接速率与总剪接计数之间的正相关以及剪接速率与总未剪接计数之间的负相关（扩展数据图 2b,c）。
Para_02
首先，我们在20个单独模拟的数据集上评估了流形学习，每个数据集包含3,000个细胞和300个基因，发现VeloCycle推断的阶段与真实情况非常接近，圆相关系数为rφ = 0.95（图2a，b）。
估计误差始终小于由后验定义的不确定性，99.2%的细胞的真实值落在5-95%的置信区间内（扩展数据图2d）。
我们还验证了特定于基因的傅里叶级数系数紧密跟踪原始真实情况（rν0 = 0.95，rν1sin = 0.98 和 rν1cos = 0.98）（图2c和扩展数据图2e）。
对于这些参数，更宽的置信区间对应于变异系数更大的噪声基因（扩展数据图2f）。
这些结果证实VeloCycle正确识别了流形几何和细胞坐标。
Fig. 2: Sensitivity analysis of VeloCycle on simulated data.
a, 细胞周期阶段分配（估计）与模拟的真实情况（GT）的散点图。
b, 估计和GT相位之间的圆相关系数的箱线图（最小值，0.932；最大值，0.963；中位数，0.957）。
c, 使用a中的数据集，估计和GT值的基因谐波（傅里叶）系数（v0, v1sin, v1cos）的散点图。
d, 使用不同数量的细胞和基因计算的估计和GT相位之间的平均圆相关系数的热图（三次模拟的平均值）。
e, 使用a中的数据集，DeepCycle获得的细胞周期阶段估计与GT的散点图。
f, DeepCycle估计和GT相位之间的圆相关系数的箱线图（最小值，0.416；最大值，0.851；中位数，0.788），跨越b中的数据集。
g, 在20次模拟中，使用VeloCycle（最小值，0.22；最大值，0.29；中位数，0.23）和DeepCycle（最小值，0.43；最大值，1.03；中位数，0.53）进行相位估计的每个细胞MSE的箱线图。
h, 使用VeloCycle（左）和DeepCycle（右）从一个模拟数据集中随机选择的30个细胞的估计（Est.）和模拟GT之间的相位差异的极坐标图。每个点代表一个细胞，线条连接估计的相位分配（浅灰色）到模拟的GT（深灰色）。
i, 估计的动力学比率与模拟GT的散点图，针对300个基因。
j, 估计和GT速度（ω = 0.4）之间百分比误差的箱线图（最小值，2.0；最大值，23.0；中位数，14.5）。
k, 说明恢复的剪接率（logβg）、降解率（logγg）、剪接计数和未剪接计数之间的关系的散点图，针对300个模拟基因。
l, 上：四个模拟中，0.0至1.5 rpmh之间16种不同模拟速度的估计和GT估计的散点图。下：上述模拟对应的后验不确定性区间的标准差点图。每个图中的黑点表示四次模拟的平均后验区间；误差带表示1倍标准差。
m, l中条件下的估计和GT速度之间百分比误差的散点图。
n, l中条件下平均未剪接-剪接表达延迟的散点图。
o, 使用不同数量的细胞和基因对三个独立模拟的速度估计范围的敏感性分析热图。每个框中的文本值表示三个数据集的平均速度，热图强度表示绝对范围。20个单独模拟数据集上的皮尔逊相关系数（r）用红色标出（a,c,e,i,k）。每个绿点代表单个基因（a,e）。每个紫点代表单个基因（c,i,k）。箱线图的边界（b,f,g,j）由四分位距（IQR）定义；须线延伸每个箱子1.5倍IQR。
Para_03
为了评估模型在不同数据集大小上的稳健性，我们进行了敏感性分析，改变了细胞和基因的数量（方法）。
我们发现估计总体上是准确的，使用最少100个细胞或100个基因即可获得大于0.70的循环相关系数（图2d）。
我们进一步将我们的推断与DeepCycle进行了基准测试，DeepCycle是一种最近基于自动编码器的方法。
这种比较显示VeloCycle通常更准确（平均均方误差(MSE)低60%，rφ=0.95），比DeepCycle（rφ=0.73）更优，尽管后者使用速度矩来实现其估计（图2e-h）。
Para_04
接下来，我们对 VeloCycle 进行了模拟相位和基因谐波的条件设置，以评估速度学习。
我们在 20 个单独模拟的数据集中观察到基因特异性动力学参数的准确估计，特别是降解-剪接速率比与真实值的匹配非常接近（rγ/β = 0.997，rβ = 0.918，rγ = 0.617；图 2i 和扩展数据图 2g,h）。
值得注意的是，VeloCycle 能够返回平均角速度的准确估计（误差率在 5.4% 至 22.6% 之间；图 2j）。
VeloCycle 在模型公式中未强加的情况下，恢复了估计的动力学参数和总计数之间的生物学相关结构（图 2k 和扩展数据图 2a,b）。
Para_05
我们进行了敏感性分析，以了解估计值在不同真实速度下的表现。
我们考虑了细胞周期速度的一个大范围，完全涵盖了生物上可能的速度范围（从 0 到 1.5 弧度每平均半衰期（rpmh），每个速度四个模拟）。
结果突显了方法性能的稳定性，估计值与真实值相差 0.2-35.8%（图 2l,m）。
在较慢的速度下，误差增加，动力学参数与真实值之间的皮尔逊相关系数较低（扩展数据图 2i，左侧）。
确实，较慢的速度对应于未剪接和已剪接 RNA 之间更短的延迟（图 2n 和方法部分），这更难以准确表征。
在所有模拟中，降解-剪接率比几乎完美匹配真实值（平均 rγ/β = 0.99）（扩展数据图 2i，右侧）。
最后，我们调查了速度学习性能是否受数据集大小的影响。
我们发现对细胞数和基因数有依赖关系，最大的数据集获得了最高的准确性和最窄的后验范围；然而，使用更多细胞可以弥补较少的基因数量，反之亦然（图 2o 和扩展数据图 2j）。
我们确定了 500 个细胞（至少 50 个基因）或 350 个基因（至少 50 个细胞）作为准确速度估计的最低限度。
Para_06
最后，我们评估了增加模拟非循环细胞比例（即每100个循环细胞中有0到200个非循环细胞）对流形学习的影响，在混合群体中，每100个循环细胞含有高达50个非循环细胞时，获得了大于0.70的环状相关性。
使用高达每100个循环细胞50个非循环细胞时，速度估计也保持在真实值的25%以内。
Manifold learning robustly estimates accurate phases
流形学习稳健地估计准确的相位
Para_01
在模拟数据上验证后，我们将 VeloCycle 部署到使用不同单细胞 RNA 测序化学方法生成的真实数据集上。
我们认为，即使类别（例如 G1、S 和 G2/M）是分类的，访问细胞周期阶段的真实情况也会有助于评估我们的阶段分配。
因此，我们对基于荧光泛素化细胞周期指示器（FUCCI）系统转导的小鼠胚胎干细胞（mES 细胞）进行了流形学习，这些细胞是通过荧光激活细胞排序（FACS）进行索引分选的。
我们使用代表广泛基因本体论（GO）查询的基因集对剪接计数进行了细胞周期阶段拟合，并将结果与 FUCCI-FACS 类别进行了比较。
属于同一类别的细胞被分配到相似的阶段（图 3a,b）；基于两个阈值并在 VeloCycle 阶段上训练的分类器在预测注释方面达到了 82.7% 的准确率，几乎匹配了在所有基因上训练逻辑分类器所获得的 87.8% 准确率（图 3c）。
此外，支持流形学习的基因拟合紧密复制了预期的细胞周期基因顺序模式。
高置信度拟合包括早期峰值组蛋白乙酰转移酶 Hat1，随后是转录因子 Trp53，以及后期促进复合体成员 Ube2c（图 3d）。
当在较小的 209 个基因集合上执行流形学习时，基因顺序和振荡幅度得到了重现，这说明该方法在灵敏度较低的化学方法上也有效（图 3e,f）。
Fig. 3: Manifold learning and gene periodicity on different datasets and technologies.
a，279 个 mES 细胞的相位分配散点图，按 FACS 分选的类别相位着色。
b，FACS 分选标签在 VeloCycle 分配的各相中的密度图。
c，条形图报告了使用仅基于 VeloCycle 相位估计的双阈值决策树获得的类别相位预测器与使用整个表达矩阵训练的逻辑回归分类器之间的比较。
d，基因拟合的代表性散点图。曲线黑线表示通过流形学习获得的基因特异性傅里叶级数。‘峰值’表示细胞周期流形上最大表达的位置（φ）。
e,f，使用小（x 轴）或大（y 轴）基因集进行流形学习时，基因峰值位置（e）和振幅（f）的散点图。
g，左：1,358 个标记基因（最小值 -0.56；最大值 0.84；中位数 0.18）和非标记基因（最小值 -1.5；最大值 0.80；中位数 -0.13）的基因振幅箱线图。右：标记基因（最小值 0.03；最大值 0.25；中位数 0.09）和非标记基因（最小值 0.02；最大值 0.31；中位数 0.15）的谐波系数不确定性箱线图。
h，振幅最大的 200 个周期性基因的类别组成饼图。
i，基因总谐波系数（ν）不确定性和振幅的散点图。基因点被标为标准‘标记’或‘非标记’。红色虚线代表标记的平均值。
j，人类成纤维细胞估计基因谐波的极坐标图。每个点代表一个基因（n = 160）。圆上的位置表示最大表达的相位，从中心的距离表示总振幅。彩色基因（橙色/绿色）是用于使用 scanpy 或 Seurat 计算标准细胞周期评分的基因。
k，早期（CDKN3, CCND1）、中期（CDC6 和 HELLS）和晚期（CDK1 和 TOP2A）细胞周期标记物的基因拟合选定散点图。
l，VeloCycle 估计的相位与 DeepCycle 相比的散点图。红色表示循环相关性。
m，VeloCycle 相位的总原始剪接 UMI 计数的散点图。黑色线条表示分箱的平均 UMI 水平。b,f,g,j 中的箱线图边界由四分位距（IQR）定义；胡须延伸每个箱子 1.5 倍的 IQR。
Para_02
鉴于我们的傅里叶参数化，我们可以通过峰值表达的相位、振荡幅度和估计不确定性对基因进行分类（补充表1）。
检查相位-幅度关系发现，通常用于Seurat和scanpy等软件包评分的标记基因（下文简称‘标准标记’）按相位聚集，这与基于FACS的地面真实情况一致（图3g,h）。
与非标记基因相比，标准标记基因平均具有更高的幅度（平均值0.14对比-0.15）和更低的后验不确定性（标准差0.26对比0.43）（图3g）；然而，在基于幅度排名前200的周期性基因中，大多数（74.5%）并非标准标记（图3h），并且其中许多（n=78）可以同样或更可靠地作为细胞周期预测因子（图3i）。
其中包括钙结合蛋白Calm2、剪接辅因子Son和细胞周期蛋白Ccnb1，它们都在细胞增殖中发挥作用。
Para_03
我们继续使用人类成纤维细胞的 10x Chromium 数据进行流形学习（图 3j,k 和补充表 2）。
为了将 VeloCycle 与其他方法进行比较，我们将它估计的阶段与通过 DeepCycle 获得的阶段进行了比较，发现有很强的对应关系（人类成纤维细胞，r = 0.882；图 3l）。
因此，VeloCycle 在不使用速度的情况下实现了与 DeepCycle 类似的阶段估计，并且与拟合单个基因谐波相结合。
作为进一步验证正确捕获了细胞周期动态的证据，我们观察到总独特分子标识符（UMIs）沿着阶段逐渐增加，随后在胞质分裂期间出现急剧下降（图 3m）。
这些结果表明，流形学习估计了一个生物学上有意义的一维几何空间，该空间跟踪了不同化学物质下的细胞周期。
Unspliced–spliced delays identify cell cycle speeds
未剪接-已剪接延迟确定细胞周期速度
Para_01
我们接下来研究了未剪接分子计数与 VeloCycle 阶段是否足够具有信息量来估计细胞周期速度。
为了在进行全面推断之前直观地探索这一点，可以独立提取未剪接和已剪接UMI的相位和基因谐波，并使用我们推导出的速度近似公式（方法）。
我们将这种方法应用于两种平行培养并在相同条件下生长的人RPE1细胞，以便我们还可以通过重复比较来评估稳健性。
首先，我们通过流形学习为每个数据集提取相位，并测量每个基因峰值未剪接表达与剪接表达之间的延迟（相位差）38（图4a）。
我们观察到基因的延迟一致且为正值（图4b），并且在重复之间具有良好的相关性（r = 0.90；图4c）。
我们将这种相关性解释为数据包含速度信息的首个证据，因此我们继续使用上述近似公式估计一个细胞周期期。
计算得出的周期是平均半衰期的18.5倍，假设实际平均半衰期为1小时，则对应18.5小时（图4d）。
Fig. 4: Analysis of delays, velocity scale, and parameter uncertainties in the choice of variational distribution.
a, 极坐标图显示了通过流形学习分析的两个RPE1细胞重复样本中106个标记基因的峰值未剪接-剪接表达。基因按其在Cyclebase 3.0中的分类注释着色。
未剪接基因拟合是分别推断的，基于运行流形学习于剪接UMI上获得的细胞阶段条件。
c, 来自a的重复样本之间的未剪接-剪接延迟散点图（r = 0.90）。
d, 使用一阶近似点估计（方法）获得的细胞周期期的条形图。
e, 两个重复样本的恒定、缩放的细胞周期速度（rpmh）的后验估计图。黑色虚线表示500个后验预测的平均值，彩色条表示置信区间（第5到第95百分位）。
f, 基因在logγg和νω之间具有不相关（左）和相关（右）后验不确定性情况下的假设情景示意图。蓝色圆圈代表高斯核分布密度；红色线条代表两个任意固定点之间的不确定性区间。
g, 使用VeloCycle的原始SVI模式对培养的人类成纤维细胞推断的速度后验估计图，使用所有基因（左）或随机基因子集（总基因数的50%；右）。
h, 使用SVI、MCMC和LRMV（SVI + LRMN）速度学习模型估计后的缩放速度（rpmh）的小提琴图。
i, 跨160个基因的降解率（logγg）和角速度（νω）后验不确定性之间的皮尔逊相关性的小提琴图。
j, 降解（logγg）和剪接（logβg）不确定性之间的皮尔逊相关性的小提琴图。
k, MCMC与SVI（顶部）或SVI + LRMN（底部）对于TOP2A和RRM2的重叠logγg−νω后验分布的密度表示。Kullback-Leibler（KL）散度得分用红色表示。h-j中的小提琴图基于500个预测样本构建；白色线条表示平均值。
Para_02
除了是一个近似值外，点估计的其他限制还包括它没有基于适当的噪声模型，并且没有与不确定性度量相关联。
为了获得更准确的估计和统计置信度度量，我们在两个RPE1重复样本上学习了完整的贝叶斯模型（速度学习），条件是在流形学习中推断出的随机变量上。
将获得的速度按拟合的平均半衰期缩放后，得到两个重复样本的平均细胞周期分别为20.1小时±0.2小时和20.0小时±0.2小时（均值±95%可信区间）（图4e和扩展数据图3a,b）。
后验分布广泛重叠（71.2%重叠），表明重复样本之间没有可信的速度差异。
为了在真实数据上验证VeloCycle可以估计沿生物学相关动态范围的细胞周期速度，我们对mES细胞进行了速度学习，这是一种快速循环的细胞类型。
对于这个数据集，VeloCycle返回了一个10.5±0.3的平均半衰期估计（扩展数据图3c）。
与RPE1细胞一样，该模型恢复了总UMI计数和基因特异性剪接和降解率之间的预期关系，正如之前在模拟中观察到的那样（扩展数据图3d和图2i）。
综合起来，这些发现证实了VeloCycle可以估计细胞周期速度并采样信息后验分布。
A structured distribution accurately models uncertainty
结构化分布准确地建模了不确定性
Para_01
尽管我们已经展示了我们的变分公式在模拟和真实数据中使用 SVI 能够恢复细胞周期阶段和速度的准确估计，但合理地质疑简化均值场变分族在表示潜在变量之间联合不确定性结构方面的局限性是合理的。
我们假设这种参数化选择可能导致对估计速度后验的过度自信，因为这些潜在变量的不确定性可能是内在相关的（图 4f）。
一个证据是观察到随机基因子集的估计值落在所有基因拟合的后验置信区间之外（图 4g）。
为了消除低估速度不确定性的偏差，我们决定通过使用马尔可夫链蒙特卡洛（MCMC；方法）对其进行采样来表征模型的联合后验。
使用无 U 转弯采样器，我们用 MCMC 研究了人类成纤维细胞的速度后验，揭示出与均值场 SVI 相比，不确定性宽了五倍（0.10 rpmh 对比 0.02 rpmh；图 4h）。
Para_02
与我们的假设一致，这个更宽的置信区间伴随着相关联的联合后验分布，捕捉了不同潜在变量之间的不确定性依赖关系。
检查后验分布时，我们发现某些基因的角速度（νω）和降解率（logγg）样本表现出相关结构（平均 r = 0.26；图 4i）。
此外，对于每个基因，我们注意到剪接（logβg）和降解（logγg）率的后验样本之间存在强烈相关性（平均 r = 0.96）（图 4j）。
Para_03
这些发现提倡重新设计我们的变分分布，以适应由MCMC推断出的后验分布的典型特征，从而保持推理准确性但避免耗时的采样过程。
我们重新制定了变分分布，将logγg和νω建模为低秩多元正态（LRMN），并将logβg建模为对应logγg的条件正态（方法部分）。
在重新训练这个新的SVI + LRMN模型后，我们获得了具有更大不确定性范围（0.08 rpmh）的速度估计，比使用平均场SVI获得的估计大（图4h-j）。
此外，我们在SVI + LRMN后验样本中检测到logγg和νω之间的相关性，这与MCMC结果中的一部分基因重叠；这导致SVI + LRMN和MCMC后验之间的Kullback-Leibler（KL）散度小于SVI和MCMC后验之间的KL散度（图4k和扩展数据图3f）。
Para_04
值得注意的是，在 SVI + LRMN 和 MCMC 中，具有高 logγg 和 νω 不确定性相关性的特定基因之间存在对应关系（扩展数据图 3g）。
logγg 和 νω 之间相关性较大的基因往往是那些未剪接-剪接延迟较大的基因（扩展数据图 3h）。
我们推测，一个基因的 logγg 和 νω 之间的依赖程度与其对速度估计的贡献程度有关。
这一推测得到了一个留一法实验的支持，其中降解率较小的单个基因对速度估计的影响最大（扩展数据图 3i）。
当将 SVI + LRMN 应用于 mES 细胞时，logγg 和 νω 后验不确定性之间的相关性也是可重复的（扩展数据图 3j,k）。
总体而言，这些实现变化导致生成了一个更稳健的模型，可以自信地用于推理，同时保留了真实后验的潜在相关结构。
Cell tracking and labeling validate the inferred velocities
细胞追踪和标记验证了推断的速度
Para_01
使用 VeloCycle 估计受流形约束的细胞周期速度最方便地用平均半衰期单位表示（方法）。
由于许多细胞类型的半衰期平均值通常是已知的，因此可以实时获得并验证整个周期中的 RNA 速度估计值。
在这方面，我们认为延时显微镜提供了一种有力的方法，将 VeloCycle 估计值与真实情况进行比较。
Para_02
为了将我们的速度估计框架与实验确定的细胞周期期进行基准测试，我们检查了一个通过延时显微镜监测的人真皮成纤维细胞（dHFs）的数据集，并收集了单细胞 RNA 测序数据（方法）。
我们的 SVI + LRMN 模型推断出一个恒定的细胞周期期为 15.3 ± 1.2 小时，假设建模转录本的平均半衰期为 1 小时（图 5a 和扩展数据图 4a-d）。
接下来，我们使用 VeloCycle 推断非恒定（周期性）细胞周期速度，并获得了类似的估计持续时间为 16.5 ± 2.1 小时，最大速度接近有丝分裂（φ 在大约 3π/2 到 2π 之间）（图 5b）。
然后，我们使用 CellPose 和 TrackMate 重建了 268 个单独细胞的细胞周期期，这些细胞随后进行了延时成像（图 5c）。
从这些数据中，我们恢复了一个中位数为 15.8 小时的细胞周期（标准差 3.1 小时），这与 VeloCycle 估计的后验可信区间重叠（图 5a、b、d）。
当使用较小的循环基因集时，也得到了类似的结果（扩展数据图 4e），并且包含不同比例的非循环 G0 细胞的敏感性分析结果与模拟观察到的适度稳健性一致（扩展数据图 4f、g）。
综上所述，这些结果表明，可以从活体成像和 VeloCycle 获得相当的细胞周期速度估计。
Fig. 5: Validation of computationally inferred velocities by cell tracking and labeling experiments.
dHFs40 细胞周期速度的后验估计图，包括常数（a）和周期性（b）。
c，上：时间延时显微镜追踪连续细胞分裂的示意图。下：在多个时间点拍摄的示例图像，用于说明通过两次分裂追踪单个分段的 dHF（粉红色）。母亲细胞分裂后（16:40 小时），一个子细胞（用白色箭头指示）被追踪了 15 小时，直到它自己分裂（31:40 小时）。
d，282 个通过活细胞成像追踪的 dHFs 的细胞周期期直方图。
e，根据分类相位分配（G1，514 个细胞；S，383 个细胞；G2M，325 个细胞）绘制的 dHF 细胞周期速度的小提琴图。中位速度由黑线表示（G1，0.35；S，0.37；G2M，0.48）。
f，未剪接–已剪接（U–S）表达延迟（左）与速度（右）之间的对应关系的双轴图。左：基因按相位分组为 20 个等分箱，以计算未剪接–已剪接延迟。实红线表示分箱平均延迟；红带表示一个标准差。右：从 b 得到的速度估计值。底部：分类相位分配概率。
g，S 和 M 相位峰值基因的基因表达散点图。垂直线对应于剪接（蓝色）和未剪接（红色）计数的峰值相位。
i，从出生（3:20 小时）到随后的分裂（20:00 小时）追踪单个 RPE1 细胞的图像。
j，337 个通过活细胞成像追踪的 RPE1 细胞的细胞周期期直方图。
k，累积 EdU/p21 实验的示意图。细胞连续暴露于 EdU，固定在不同的时间点，并进行 EdU 检测和 p21 免疫染色。
l，左：累计 EdU 标记 2 小时、8 小时和 36 小时后的 p21（绿色）、4,6-二氨基-2-苯基吲哚（DAPI）（青色）和 EdU（洋红色）染色图像（来自三次实验重复之一）。比例尺，100 μm。右：不同染色组合的单个细胞图像。比例尺，10 μm。
n，点图表示不同时间点 p21+ 细胞的平均百分比。黑色水平线表示平均值（最小值，0.24；最大值，0.76；平均值，0.52），误差条表示标准差；每个点代表单个重复中 p21+ 细胞的百分比（n = 29；总共三个重复，分别有十个、十个和九个时间点）。
o，上：13 个时间点（从 30 分钟到 73 小时）EdU+ 细胞的比例线图。数据显示三次重复的平均值（除 2 小时外，其余均来自两次重复），误差条表示标准差。下：静止细胞（p21+）中 EdU 阳性细胞的比例随时间变化的线图。A，生长和平台期交点的 x 值；B，线性拟合的 y 轴截距；GF，平台期的 y 值。
p，scEU-seq 实验设计的插图，该设计生成 24 个表，供 Dynamo 使用以产生细胞周期期的标准估计。RFP，红色荧光蛋白；GFP，绿色荧光蛋白；RPEs，视网膜色素上皮。
q，左：VeloCycle、Dynamo（无代谢标记信息）和 Dynamo-Metabolic（有代谢标记信息）模型采用的不同实验测量、流形和细胞路径推断方法的示意图。右：显示 VeloCycle、Dynamo 和 Dynamo-Metabolic 模型获得的估计细胞周期期的图。小提琴图显示 VeloCycle 输出的后验分布，圆圈是个别 LAP 评估；红色星号表示平均值。红色虚线表示 d 和 j 中的中位数。白色虚线表示 500 次后验预测的平均值，黑色条表示可信区间（第 5 到第 95 百分位）在 a、b、f（右）和 h 中。NS，不显著；**P < 0.01。
Para_03
我们接下来根据独立的细胞周期阶段对速度进行了分类，以便更细致地评估这些评价和模型行为。
我们观察到，在G2/M期（平均缩放速度为0.47 rpmh）细胞周期的进展比G1期（0.37 rpmh）和S期（0.36 rpmh）更快（图5e）。
动力学参数及其后验不确定性在恒定速度模型和周期性速度模型之间高度相关（扩展数据图4h，i）。
值得注意的是，当估计不同细胞周期阶段峰值基因的平均未剪接-剪接延迟时，我们发现具有较大平均延迟的细胞周期阶段对应于速度较快的区域（图5f）。
延迟较大的基因也是那些剪接和降解速率较小的基因，这符合近似模型（扩展数据图4j，k和方法部分）。
在检查了未剪接-剪接延迟和角速度与降解率之间的低秩基因特异性后验相关性后，我们可以识别出最强烈贡献于潜在速度估计的具体基因（图5g）。
Para_04
为了进一步审查 VeloCycle 推断的细胞周期持续时间与实验获得的时间匹配的程度，我们对同一培养的 RPE1 细胞进行了延时显微镜检查和单细胞 RNA 测序。
使用 VeloCycle 获得的速度约为 17.7 ± 2.1 小时（图 5h 和方法）。
与 dHFs 一样，这一计算估计值与通过延时成像追踪分裂细胞获得的平均细胞周期持续时间 17.7 小时（标准差 3.4 小时）重叠（338 个细胞）（图 5i,j）。
接下来，我们试图将从延时显微镜和 VeloCycle 获得的细胞周期持续时间测量结果与使用正交实验技术获得的结果进行比较。
因此，我们进行了连续的 5-乙炔基-2′-脱氧尿苷 (EdU) 标记，以独立估算细胞周期长度（图 5k,l）。
在 S 期进行 DNA 复制时，循环细胞会吸收 EdU；因此，EdU 脉冲的持续时间与 EdU 阳性细胞的比例直接相关。
在 72 小时内监测了 13 个时间点的 EdU 水平（图 5m），我们使用 p21 (CDKN1A) 染色来考虑处于 G0 期的细胞，并确定平均细胞周期长度为 16.8 小时（图 5n,o 和方法）。
综合起来，这些发现验证了在细胞周期背景下计算的 RNA 速度估计值。
据我们所知，这是首次通过实验方法直接验证 RNA 速度估计的例子，这证明了可以将 VeloCycle 输出用作实际（而非伪）时间单位。
Para_05
另一种进一步验证细胞周期期估计的方法是考虑代谢标记的单细胞 5-乙炔基尿嘧啶（EU）测序（scEU-seq），这提供了显著更丰富的信息内容。
使用 Dynamo 的最小作用路径估计（LAP）程序，特别设计用于获得代谢建模的速度，我们同时处理了最近一次多脉冲代谢标记实验获得的 24 个数据矩阵（包括剪接和未剪接的数据，标记和未标记的，跨越六个脉冲长度）。
值得注意的是，当我们比较这些 LAP 结果与仅基于剪接-未剪接矩阵的 VeloCycle 估计时，VeloCycle 的后验概率与 Dynamo 估计密切重叠，后者使用了更具信息量的标记数据集，并进一步利用 FUCCI 染色信息作为真实嵌入（图 5p,q 和扩展数据图 5）。
这些发现进一步证实了 VeloCycle 在标准 scRNA-seq 数据上使用金标准估计对细胞周期速度的估计。
VeloCycle benchmarks competitively across multiple datasets
VeloCycle 在多个数据集上具有竞争力的表现
Para_01
尽管将 VeloCycle 与具有不同范围和假设的方法进行比较存在概念上的挑战，我们还是在四个独立的数据集（模拟数据（图 2），dHFs（图 5a），代谢标记的 A549 细胞和 RPE1 细胞（图 4h））上对 VeloCycle 进行了性能基准测试，对比了四种 RNA 速度估计方法（scvelo、cellDancer、Dynamo 和 VeloVAE）。
VeloCycle 在多个数据集中相对于现有方法实现了显著提高的跨边界方向正确性（CBDir；扩展数据图 6a,b）和速度一致性得分（扩展数据图 6c），即使使用它们的嵌入来生成用于评估的真实簇也是如此（扩展数据图 6b）。
这种基准测试分析还揭示了 VeloCycle 在剪接和未剪接拟合上的均方误差（MSE）更高，这在我们的预期之中，因为它反映了我们选择优先考虑正则化而非误差最小化的策略（扩展数据图 6d–g）。
VeloCycle can test for drug treatment effects on velocity
VeloCycle 可以测试药物治疗对速度的影响
Para_01
现有的 RNA 速度框架没有提出一种方法来测试所获得估计值的统计显著性，这可能是因为在基因层面的速度参数化下这一任务具有挑战性。
例如，目前无法确定接近零的 RNA 速度估计值是否仅仅应该被解释为噪声。
此外，两个样本之间速度估计值的直接比较无法通过置信度量来支持。
借助 VeloCycle，首次可以对速度进行统计推断，既可以针对特定的零假设，也可以用于细胞群体之间的差异速度显著性分析。
Para_02
为了说明我们的模型如何用于实际中的统计速度测试，我们在使用药物厄洛替尼治疗前（D0）和治疗后（D3）对PC9腺癌细胞系进行了RNA速度分析（扩展数据图7a-e）。
在贝叶斯设置中进行统计检验可以通过从后验分布计算可信区间来实现。
首先，我们考虑了D0细胞的速度后验分布，以询问是否存在非零速度的统计支持。
鉴于可信区间与零没有重叠，我们可以得出结论，数据包含细胞周期进展的统计显著证据（图6a，左）。
然后，我们将处理样本（D3）与对照组（D0）进行了比较。
我们发现时间点之间存在显著的速度差异，在D3检测到较慢的有丝分裂细胞周期速度（图6b）。
例如，我们按相位间隔分层，并检查了后验样本，确认在D3相比D0期间G2/M期速度降低，但在G1和S期则没有（图6b和扩展数据图7f）。
厄洛替尼治疗后有丝分裂细胞数量减少的进一步证据来自被分配M期坐标的D3细胞密度较低（扩展数据图7a，底部）。
Fig. 6: Statistical velocity inference across diverse biological contexts and with transfer learning in genome-wide perturbation screens.
PC9 肺腺癌细胞（D0）相对于零速度对照（红色）的缩放速度的事后估计图。
厄洛替尼治疗前（D0）和治疗后（D3）的缩放速度的事后估计图。区间不重叠的区域表明存在统计学上显著的速度差异。底部：分类相位分配概率。
D0 和 D3 样本之间 273 个基因的平均未剪接-剪接表达延迟的散点图。基因点按峰值表达相位着色。
小鼠 FB、MB 和 HB RG 前体细胞在发育阶段 E10 获得的缩放速度估计的小提琴图。
使用 BoneFight 将单细胞簇的空间投影到四个用 HybISS 分析的参考 E11 胚胎部分上，按速度估计着色。区域域（FB、MB 和 HB）和脑室区（VZ）相应地标记。
E10（h）和 E14-15（i）沿细胞周期流形的细胞分布的核密度估计（KDE）图，按区域域着色。
CRISPR 诱导单基因敲除后 7 天的 RPE1 细胞的细胞周期速度的事后估计图，按 NT 对照（绿色）和细胞周期敲除（米色）条件分层。流形学习使用大（顶部）或小（底部）基因集进行。
来自 j 的 NT 和细胞周期敲除（CC-KD）样本的连续相位分布的核密度图。
所采用的迁移学习方法的示意图。基因谐波系数是在 NT 对照（许多细胞）上获得的，并应用于特定基因敲除条件下具有少量或分布不均的细胞的相位分配。
167,119 个 RPE1 细胞中 986 个单独基因敲除（Δ）条件的速度学习事后估计和标准差的散点图。垂直线对应于 NT（绿色）、细胞周期标记（浅褐色）和其他基因敲除的标准速度估计。
NT、MCM6Δ 和 DBR1Δ 条件下的 KDE（顶部）和未剪接-剪接延迟（底部）图。深绿色线表示平均延迟；浅绿色线表示标准差。
使用小基因集和大基因集为 m 中条件获得的缩放细胞周期速度估计的散点图。
m 中每个条件的总细胞数和事后速度标准差的散点图。在 a、b 和 j 中，黑色虚线代表 500 次事后预测的平均估计值；条形图表示可信区间（第 5 百分位至第 95 百分位）。在 d 和 f 中，黑线表示 500 次事后预测的平均值。在 o 和 p 中，皮尔逊相关系数以红色显示。
Para_03
由于未剪接-剪接延迟与细胞周期速度相关，我们假设D0和D3时间点之间会有不同的延迟，特别是对于在有丝分裂期间达到峰值的基因。
在计算厄洛替尼治疗前后的基因特异性未剪接-剪接延迟后，我们确实注意到在有丝分裂期间表达高峰的一组基因在D0时的相位延迟比D3时更大（图6c）；这包括后期促进复合体成员CDC27（差异延迟(dd) = 0.11弧度），细胞周期依赖性激酶抑制剂CDKN3（dd = 0.10）和中心体支架因子ODF2（dd = 0.09）（扩展数据图7g）。
M期细胞周期速度的降低与厄洛替尼（一种EGF阻断剂，抑制向G1期的进展）预期的效果一致。
这一结果也与证据相符，即PC9细胞系不应观察到完全停滞，据报道该细胞系对完全阻断具有一定的抗性。
Cell cycle speed varies spatiotemporally in radial glia
细胞周期速度在径向胶质细胞中时空变化
Para_01
调节增殖速率以及脑室区径向胶质细胞（RG）的对称和不对称分裂在大脑前-后轴上的发育时间控制中起着关键作用。
为了阐明在小鼠神经发育过程中，占据不同空间区域的祖细胞之间的细胞周期速度是否存在差异，我们在胚胎第10天（E10）阶段对前脑（FB）、中脑（MB）和后脑（HB）的RG细胞进行了VeloCycle估计。
细胞周期速度沿前脑-中脑-后脑轴变化，后部（HB）的祖细胞比前部（FB）的分裂更快（图6d）。
通过使用原位杂交空间转录组学（HybISS）数据和BoneFight算法将这些细胞推断出的细胞周期速度映射到相应位置，可以更细致地观察这一梯度（方法）。
我们观察到靠近脑室区的位置存在快速分裂的RG细胞，这表明细胞增殖沿着脑室区发生，并暗示该区域的不同部分以不同的速率增殖（图6e）。
相比之下，在E14和E15时，来自所有三个区域的RG细胞稳定在同一增殖速度，没有可信的速度差异（图6f）。
在这些较晚的时间点，大多数MB和HB区域的RG细胞积累在非增殖状态；而大多数存在的RG细胞来自FB，这在E10时发展得更慢（图6g-i）。
这些结果与最近的研究一致，显示HB更快地分化为非增殖、分化的细胞类型；因此，在早期发育阶段可能需要更高的增殖能力。
此外，后期的减速是预期中的，并且与EdU追踪研究中的报道一致。
Speed modulation screening is achieved by transfer learning
通过迁移学习实现速度调制筛选
Para_01
先前的 RNA 速度框架在使用仅有限数量细胞的细胞类型或条件时难以获得可靠的估计。
随着最近设计用于筛选数百种遗传、环境或药物扰动的单细胞技术的发展，使用少量细胞评估细胞动力学变化的需求日益增长。
VeloCycle 可以探索以前具有挑战性的 RNA 速度背景：通过将流形学习模型基于从大型参考数据集中预先推断出的基因谐波系数进行条件化，可以使用更少的细胞或属于单一细胞周期阶段的细胞进行速度推断（方法）。
Para_02
为了证明这一点，我们研究了一个大规模的全基因组 Perturb-seq 数据集，在该数据集中，通过靶向的 CRISPR 干扰文库将数百个单个基因敲低引入了 RPE1 细胞系，并在培养 7 天后进行了单细胞 RNA 测序。
首先，我们在非靶向（NT）对照细胞和一组对应于细胞周期（CC-KD）的特征性标记基因的基因敲低条件下运行了 VeloCycle。
NT 细胞的细胞周期期为 25.6 ± 1.3 小时，而 CC-KD 细胞的细胞周期期为 30.9 ± 1.3 小时（图 6j），使用了两个不同大小的基因集合（扩展数据图 7h-l）。
当根据通常被认为是 S 和 G2/M 标记的基因对 CC-KD 条件进行分层时，与 NT 细胞相比，观察到 G1 期空间中的细胞积累（图 6k 和扩展数据图 7m）。
这表明，某些个别与细胞周期相关的基因的功能丧失会扰乱细胞周期进程，要么是在某些阶段减缓增殖速度，要么是在特定进入检查点前完全停止进展。
Para_03
为了研究单个基因敲低条件对细胞周期速度的影响，我们采用了一种迁移学习方法，该方法基于从 NT 和 CC-KD 数据子集中先前推断出的基因谐波进行流形学习，为167,119个细胞和986个个体条件分配了阶段，其中一些条件的细胞数量少至75个（图6l）。
与数据的粗略分层一致，我们观察到，在个别与细胞周期相关的基因敲低条件下，细胞周期速度显著下降，这与NT对照细胞和非细胞周期相关基因敲低的细胞相比尤为明显（图6m）。
几个受影响最严重的细胞周期速度出现在DNA复制（MCM3Δ和MCM6Δ）和翻译起始（E1F3BΔ、EIF2B3Δ和EIF3CLΔ）高度表征基因的敲低中。
有趣的是，几种剪接和mRNA处理基因的敲低条件要么显著降低了估计的细胞周期速度，要么增加了速度，包括内含子套索剪接因子DBR1Δ（与NT条件相比减少了11.7倍），PRPF3Δ（增加了1.2倍）和PRPF31Δ（增加了1.3倍）（图6m,n）。
鉴于RNA速度估计依赖于RNA代谢生命周期的控制微分方程，这一结果表明，与RNA代谢相关的生物扰动破坏了速度框架的生物物理参数化。
此外，每个条件下数据集中的细胞数量直接影响速度不确定性的估计，这意味着更多的细胞，以及因此对于某个条件聚集的稀疏性较少，提高了VeloCycle模型对其获得的速度估计的信心（图6o,p）。
最终，这些分析表明，通过迁移学习方法，速度可以作为评估基因敲低对生物过程动态影响的度量，在大规模扰动背景下应用。
Statistical velocity generalizes across manifold geometries
统计速度在流形几何中具有普遍性
Para_01
VeloCycle 是一种概率速度模型，设计用于一维周期流形（图 1）；然而，我们新的流形约束框架也可以用于生成更高维度和各种几何形状的公式，包括一维（非周期）区间和二维情况（方法和讨论）。
为了探索这种可能性，我们制定了两个新模型，并进行了关键测试：一个是一维区间（非周期）模型，旨在研究分化速度（扩展数据图 8 和 9），另一个是适用于检查更复杂设置的二维模型（扩展数据图 10）。
Para_02
对于一维区间模型，我们使用独立估计的伪时间（扩散伪时间）定义了流形坐标，并使用B样条基函数定义了流形几何，而不是用于VeloCycle的周期性傅里叶级数基（扩展数据图8a）。
首先，我们通过十个独立模拟的数据集验证了基于伪时间的模型性能，正如对VeloCycle所做的那样（图2），准确地恢复了模拟的真实速度和动力学参数（扩展数据图8b-h）。
接下来，我们展示了这一原理验证扩展在胰腺内分泌发生过程中推断速度的能力，以及在小鼠齿状回中的应用（扩展数据图8i-r和方法部分）。
这证明了该框架的潜力，可以分别推断在同一样本内同时发生的两种不同生物过程的速度，即细胞周期和β细胞分化（扩展数据图8l,o）。
最后，我们在模拟数据上评估了一个更复杂的二维案例（扩展数据图10a），成功地同时恢复了描述两个过程（分化和分歧）的真实速度和动力学参数（扩展数据图10b-g）。
总体而言，我们证明了我们的流形约束RNA速度框架可以适应于除周期性情况之外的其他模型；然而，为了提供一个更强大和独立的工具，还需要像对VeloCycle所做的那样进一步验证和表征这些模型。
Discussion
Para_01
在这项工作中，我们通过设计一个框架来解决当前RNA速度方法的几个限制，该框架将流形和速度推断统一到一个单一的概率生成模型中。
我们注意到，为了获得更平滑、更不易过拟合的向量场，已经有人提出并在低维嵌入上应用了速度场的投影和平滑方法，这特别适用于可视化和数据探索目的。
在这里，我们提出了、测试并实验验证了一种显式的RNA速度参数化方法，作为定义在流形坐标上的向量场，从一开始就考虑了切线性作为其核心假设之一。
另一个最近提出的方法graph-dynamo也认识到了拥有切线于流形结构的速度估计的重要性，其切线空间投影在后处理中用于校正向量场以进行可视化和解释。
VeloCycle进一步强调了速度切线性的重要性，通过将流形约束直接整合到估计过程中，为参数识别和推断提供了可能。
Para_02
我们的框架使用变分推断直接从原始数据推断我们生成的RNA速度模型的后验参数，并适当地对数据中的噪声进行建模，而不是使用最近邻平滑等启发式方法。
VeloCycle 返回不确定性估计，使得可以直接评估对估计结果的信心以及样本之间的细胞周期速度比较。
这些能力在不同的生物学环境中都很重要，例如在癌症生物学中，需要使用快照单细胞数据来仔细检查细胞周期进展的改变。
RNA速度以前曾被用于说明细胞周期进展，但这些方法需要多个启发式步骤，并且仅具有探索性价值，因为无法从推理过程中得出结论。
因此，VeloCycle 可能会为疾病进展提供新的生物学见解，例如通过表征不同微环境或患者之间肿瘤增殖率的差异。
Para_03
最早引入贝叶斯变分推断用于RNA速度建模的方法，VeloVAE、VeloVI和Pyro-Velocity，通过简化变分分布的方式限制了估计联合后验的实用性，特别是在未缩放的基因特异性速度参数化的情况下。
更广泛地说，具有高自由度的模型和独立性假设可能会过度拟合噪声并高估速度的置信度。
我们预计，未来的速度方法将进一步探索我们在共享速度函数下约束剪切-未剪切拟合以及控制不确定性依赖性的策略。
Para_04
关于与其他方法的比较，总体上对 VeloCycle 相当有利，但这些基准不应被过度解读。
例如，VeloCycle 获得的最佳速度一致性和跨边界方向正确性（CBDir）分数是我们自洽公式预期的结果。
类似地，我们将基因表达拟合中的较高均方误差（MSE）值解释为获得推理能力所必需的成本。
总体而言，我们主张根据方法的预期应用、已证明的有效性和基本目标来选择方法，而不是基于粗略的指标。
Para_05
VeloCycle 通过将流形的约束直接纳入 RNA 速度估计过程中，丰富了速度分析。
这通过统一动力学参数估计和流形的内在几何结构来实现结构化正则化，确保估计量的一致性。
此外，我们的方法通过提供一个从原始数据拟合的端到端统一模型而脱颖而出，避免了启发式方法，促进了可解释性，并具有严格的推理能力。
对这一模型的主要潜在变量（流形几何、速度和动力学参数）进行推理，避免了与多个基因特异性速度相关的陷阱。
具体来说，它不会导致因将每个基因视为独立而忽视它们不确定性相关性所引起的过度自信。
Para_06
虽然我们用于 RNA 速度估计的模型提供了明显的好处，但仍有许多进一步发展的途径。
首先，VeloCycle 集中于 1D 周期流形的情况，但可以自然地扩展到不同维度和拓扑的潜在空间。
我们还展示了其在 2D 情况下的适用性；然而，如我们对 VeloCycle 所展示的那样，这些模型还需要进一步的实验验证和特征描述。
Para_07
其次，定义维度的问题与基因选择问题相交；由不同基因定义的不同子空间暴露出被不同比例的细胞穿越的不同的流形。
针对这一问题开发的方法最近已被提出，并通过适当的修改，这些方法可以整合到RNA速度估计方法中，以自动选择拓扑和基因集。
在这方面，考虑具有不同拓扑结构的多个流形的框架，这些流形由不同子空间中的细胞跨越，同时还将特定的细胞和基因分配给这些特征，将显著增强流形一致的RNA速度估计的一般适用性和实用性。
这可能涉及使用无偏或GO信息的数据因子化模型来提取和研究沿这些模块的速度。
这种类型的基因模块建模方法已通过cell2fate方法成功演示。
Para_08
第三，我们的模型假设基因特异性剪接和降解速率是恒定的；实际上，对于某些基因，在细胞周期的不同阶段，这些速率可能会发生变化。
未来对 VeloCycle 的扩展可以通过参数化函数定义动力学参数来解决这一限制。
然而，在这些设置下保持模型的良好条件可能并非易事。
Para_09
第四，未来对 VeloCycle 的改进可能涉及添加结构化先验和其他约束：对于细胞周期，可以实施辅助损失以支持细胞因子分裂后总 UMI 下降的配置。
更广泛地说，引入熵正则化可以鼓励细胞在流形上的均匀分布，施加低秩约束于基系数可以改善高维流形的学习。
同样，假设特定基因序列激活或利用基因调控网络知识的先验可以注入有价值的生物学信息。
将基因组特征作为动力学参数的预测因子是另一种有前景的模型正则化策略。
最后，我们讨论了自主动力系统 V(x) 的情况；原则上，由于实验设计使用外部干预或扰动，非自主情况也可以在这个框架内设想，但公式和实现可能需要逐案、非通用的发展。
广泛使用的标准分析流程依赖于一小群标记基因来为单个细胞分配一个分类相位，即使细胞周期进展是一个连续的过程。
最近的方法推断连续相位分配，在基于评分的方法上代表了有价值的改进。
VeloCycle 的流形学习在这方面取得了进展，同时推断个体基因周期模式，提供后验不确定性，并获得与其他方法相比更有利的结果。
值得注意的是，流形学习步骤是灵活的，促进了迁移学习：可以在更大或更高质量的数据集上估计流形几何结构，作为较小数据集的先验。
这增强了速度学习在不同实验条件下的稳健性和适用性。
鉴于条形码策略在单细胞水平筛选中的使用增加，这一点尤其相关。
我们预计此类模型将在药物筛选和评估遗传变化对异质细胞池的影响方面得到未来应用。
验证 RNA 速度向量场整体一致性的方法是将基于启发法估计的人群间转换概率与先前关于它们谱系关系的知识进行相关性分析；然而，这是相关性和间接的。
在这里，我们转而直接将估计值与过程的实际速度进行比较。
通过特定的生物学原理先验，VeloCycle 获得的速度可以解释为增殖速度，这可以在不同的组织位置、发育的不同时刻或由于核心基因调控网络的扰动而有所不同。
尽管我们提倡使用代谢标记技术，这些技术通过更富有信息量的实验设计、靶向化学和结构化数据，可能允许更好地估计速度动力学参数，但这种做法的应用受到限制。
因此，设计实验验证或流形约束的 RNA 速度方法（如 VeloCycle）是一项重要的努力，具有社区普遍应用的目的。
最终，我们的框架代表了从单细胞数据进行动力学估计的严谨性的进步。
将 RNA 速度定制到单一过程中的有希望的结果，倡导开发新的模型，将单细胞数据的高维度分解为具有相应和可解释的 RNA 速度场的独立生物轴。
Methods
Model specifications for manifold-constrained RNA velocity
流形约束RNA速度的模型规范
The manifold
歧管
Measurements and noise model
测量和噪声模型
RNA velocity and chemical kinetics
RNA速度和化学动力学
Latent-space dynamics
潜在空间动力学
Para_03
Manifold-constrained RNA velocity
流形约束的RNA速度
Para_01
我们现在可以将方程 (1)–(4) 连接起来以获得
Para_03
方程(5)为我们方法提供了基础，因为它将左侧低维流形的拓扑结构与右侧的生物学联系起来。
值得注意的是，控制基因动力学的参数（β和γ）原则上也可能依赖于x。
Geometric interpretation
几何解释
u(x) and inference
u(x) 和推理
Duration of biological processes
生物过程的持续时间
Para_01
这种表述的好处在于，可以从轨迹和 V(x) 估计生物过程的实际持续时间：
Manifolds with S 1 topology: the cell cycle
具有S^1拓扑的流形：细胞周期
Likelihoods
可能性
Bayesian model formulation for VeloCycle
VeloCycle 的贝叶斯模型构建
Para_01
我们的模型包括生物定义的先验与从数据中确定的经验贝叶斯风格的先验。
我们的目标是基于上述似然表达式来估计联合后验概率分布的近似值。
Para_03
Para_05
通过找到最大化 corr(Φ_c, ∑_g S_gc) 的全局相位偏移来获得旋转不变性（例如，第一个单元 c0 的任意性，使得 Φ_c0 = 0）。
投影正态分布的集中参数 ε 默认设置为 5，但可以根据对数据质量的整体信心进行调整。
Variational distribution: SVI
变分分布：SVI
Para_01
我们在基础模型中使用的变分分布是平均场分布，其边缘分布可以是正态分布或狄拉克δ分布。
Para_02
Variational distribution: LRMN
变分分布：LRMN
Para_02
Model implementation
模型实现
Para_02
对于流形学习，我们估计每个细胞沿圆形细胞周期流形（φ）的位置以及用于建模剪接计数对数期望（ElogS）的每个基因的傅里叶级数系数（ν），这些模型是从实际数据和负二项分布（NB）建立的。
我们将所有变量初始化为其先验的均值，该均值是使用前两个主成分（φ）或每个基因剪接表达的均值和标准差（ν）确定的。
为了允许不同数据集或批次之间平均表达水平的差异，我们还定义了第一个基因谐波系数的偏移项（Δν）。
Para_03
对于速度学习，我们推断角速度（νω）的傅里叶系数以及速度动力学参数（γ和β），条件是流形学习期间获得的参数后验估计的均值。这些变量用于建模未剪接计数对数的期望值（ElogU），而这些值本身是从真实数据和负二项分布（NB）建模而来。我们将所有变量初始化为先验的均值，角速度的先验均值为零（假设细胞周期速度为零）。为了在学习过程中确保方程（10）中的正数（ω(φ)∑fνgf ∂φζf(φ) + γg），我们使用了ReLU函数。
Overview of VeloCycle latent variables
VeloCycle潜在变量概述
[div_table]
Para_01
给定数据，我们使用 SVI 解决 VeloCycle 模型，并应用 ClippedAdam 优化器和 ELBO 损失函数，学习率从第一个训练迭代的 0.03 衰减到最后一个训练迭代的 0.005。
通常，我们进行 5,000 次训练迭代用于流形学习和 10,000 次训练迭代用于速度学习；然而，提供了一个提前终止训练的选项，如果前 100 次迭代的平均损失与前 10 次迭代的平均损失相差不到五个单位，则不再执行进一步的迭代。
Para_02
在执行 MCMC 时，我们使用一个 No-U-Turn（NUTS）内核，从通过 SVI 首先获得的后验均值估计开始。我们通常使用一条链，2000 步热身采样和 500 步实际采样。
Para_03
VeloCycle 可以使用本地 CPU 或 GPU 在几分钟内运行，使用 GPU 时运行速度显著提高，特别是在使用大量细胞（>30,000 个细胞）或基因（>300 个基因）时。
由于沿基因维度有更多的参数，增加基因数量比增加细胞数量更能快速减少运行时间。
Biological constraints on parameters
生物参数的约束
Para_01
速度动力学参数 β 和 γ 受生物学限制。特别是，
Para_02
此外，每个基因的谐波系数先验是根据数据中所有细胞的表达水平均值和方差确定的。
对于速度谐波系数，我们假设先验均值为无速度（即0），标准差较宽（3.0）。
Para_03
所有先验都可以使用 'velocycle.preprocessing' 函数套件轻松修改，并通过元参数（'mp'）项提供给 Pyro 模型对象。
Approximate point estimate for constant cell cycle velocity
恒定细胞周期速度的近似点估计
Para_02
