生信程序 | VeloCycle：使用流形约束的RNA速度模型进行统计推断揭示了细胞周期速度的调制

生信菜鸟团

发布于 2025-01-02 14:44:29

1830

文章被收录于专栏：生信菜鸟团生信菜鸟团

Basic Information

英文标题： Statistical inference with a manifold-constrained RNA velocity model uncovers cell cycle speed modulations
中文标题：使用流形约束的RNA速度模型进行统计推断揭示了细胞周期速度的调制
发表日期：31 October 2024
文章类型：Article
所属期刊：Nature Methods
文章作者：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速度场。
b, 由标准方法描述的无约束速度估计的示意图。
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上获得的细胞阶段条件。
b, 未剪接-剪接延迟（以弧度计）的直方图。
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 相位峰值基因的基因表达散点图。垂直线对应于剪接（蓝色）和未剪接（红色）计数的峰值相位。
h，RPE1 细胞周期速度的后验估计图。
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。
m，细胞周期进展期间累积 EdU 标记的示意图。
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）相应地标记。
E14-E15 区域域的速度估计的小提琴图。
RG 各阶段的区域比例条形图。
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速度统计评估的核心。
最早引入贝叶斯变分推断用于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

V(x) 是描述低维潜在空间动力学的向量场。

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

假设我们现在有多个条件（或重复实验）c，并且寿命 τ_g = γ_g^(-1) 是条件独立的，我们观察到这种关系

Gene sets and quality control filtering

基因集和质量控制过滤

Para_01

为了选择预期会在细胞周期中表现出周期性行为的基因进行速度分析，我们应用了三个不同大小的基于文献的循环基因集：‘小’包含 97 个基因，‘中’包含 218 个基因，以及‘大’包含 1,918 个基因。
VeloCycle 默认使用‘中’基因集，以最小化噪声或低表达基因对流形和速度估计的影响；然而，在测序深度和数据集质量特别高的情况下，我们也使用了‘大’基因集。
函数‘velocycle.utils.get_cycling_gene_set’可用于访问这些人类和小鼠基因集。
还根据剪接和未剪接计数的平均检测进行了额外的基因过滤，具体如以下部分所述。

Categorical and continuous cell cycle phase assignment

分类和连续的细胞周期阶段分配

Para_01

使用 scanpy 函数 ‘sc.tl.score_genes_cell_cycle’ 进行了细胞周期阶段（G1、S 和 G2/M）的分类分配，如前所述。
使用从 ‘scvelo.pp.moments’ 获得的速度信息和原始出版物中描述的标准参数，通过 DeepCycle 在模拟和真实数据集上实现了连续的细胞周期阶段分配。

Inference of the unspliced–spliced delay

未剪接-已剪接延迟的推断

Para_01

为了从 VeloCycle 的结果中计算未剪接-已剪接的延迟，我们基于每个基因计算了未剪接和已剪接 UMI 峰值表达相位之间的差异（以弧度为单位），使用未剪接（ElogU）和已剪接（ElogS）计数的估计期望值。

Posterior probability sampling

后验概率抽样

Para_01

除非另有说明，潜在变量及其相关估计不确定性是在模型训练后从 500 个后验样本中收集的，使用了‘pyro.infer.predictive’，可信区间测量在第 5 百分位和第 95 百分位之间。

Para_02

从速度函数 ω(φ) 获得的细胞周期速度估计值按平均降解半衰期进行缩放，即 mean(γg)。
为了推断整个细胞周期的细胞周期期，我们在 20 个相位（从 0 到 2π）的网格上从速度函数取样，并使用 'scipy.integrate.trapz' 计算曲线下面积。
然后使用 'numpy.mean' 和 'numpy.percentile' 计算后验均值、第 5 百分位数和第 95 百分位数。
通过计算第 95 百分位数和第 5 百分位数估计值之间的差异来获得后验估计的完全不确定性范围。

Extension to 1D nonperiodic intervals and 2D manifolds

扩展到一维非周期区间和二维流形

Para_01

我们探索了模型向其他一维或二维流形扩展的可能性。
为了适应非周期性，需要进行的主要更改是用其他平滑函数替代傅里叶级数。
为了确保我们的框架适用于不同的生物相关场景，我们对函数形式的假设限制在其一阶和二阶导数的连续性上。
特别是，我们使用三次B样条基（B_3）而不是傅里叶分量：

Para_02

这里，φi，i ∈ {1,2} 表示流形 φ 的两个维度。我们注意到，在一维情况下，p 阶（在我们的情况下 p = 3）B-样条基和基维度 m 是通过以下方程定义在一个变量 t 上，t 有一系列节点 ti，i ∈ {1, …, m + p + 1}：

Para_03

在二维情况下，拼接计数和速度都是使用二维样条建模的，我们通过常规一维 B 样条的笛卡尔积获得这些样条。
尽管这个向量场可能无法完全代表从势能景观衍生出来的分支过程，其中包含临界点，但它足以说明如何原则上处理二维问题。

Para_04

一维和二维样条模型都在 Pyro 中实现。
估计分为两个步骤进行。
第一步，我们根据选定的 φ 条件拟合拼接计数的样条系数。
第二步，我们基于未拼接计数估计动力学参数和速度。
对于样条系数，我们选择了一个零均值的宽正态先验。
在一维模型中，拼接计数的 B 样条基具有五个维度，速度 ω(φ) 被建模为一个标量值，在分化过程中保持不变。
在二维模型中，拼接计数和速度的样条在每个轴上都有六个基维度。
我们使用平均场变分分布，并将所有变量初始化为先验的均值。
该模型使用 SVI 进行训练，优化器可以选择 Adam 或 ClippedAdam。

Structured data simulations and sensitivity analyses of VeloCycle

VeloCycle 的结构数据模拟和敏感性分析

Para_01

为了正确验证 VeloCycle 在所有流形学习和速度学习潜在变量都有真实值的数据集上的性能，我们采用了一种新的结构化模拟方法来保持速度动力学参数（剪接率 β 和降解率 γ）和基因谐波（ν0、ν1sin 和 ν1cos）之间的关系。
这些关系在真实数据中是预期存在的，并且在模拟中是必要的，以避免不现实的情况，即未剪接与已剪接计数的比例不切实际地过高或过低。
我们预计，包含更多速度信息的基因应该是那些具有较大未剪接-已剪接延迟和较慢剪接及降解速率的基因；动力学过快的基因在单细胞 RNA 测序数据中提供的信号将非常有限。
因此，我们构建了一个生成性的 VeloCycle 模型，该模型在基因谐波和速度动力学速率参数之间施加了相关性结构，唯一目的是采样模拟数据（而不是用于推断本身）。
我们定义的相关性如下：基因谐波系数之间存在弱正相关（r = 0.05；假设每个基因只有一个正弦和余弦项），剪接率与零阶基因谐波系数 ν0 之间存在中等正相关（r = 0.30），以及剪接率与降解率之间存在中等正相关（r = 0.30）。

Para_02

使用这个相关矩阵，通过从‘pyro.dist.MultivariateNormal’中随机抽样生成了模拟数据集，用户可以定义基因和细胞的数量。
这些变量，连同用户定义的细胞周期速度（νω）和从0到2π之间的随机均匀分布中抽取的特定细胞相位（φ），被代入速度方程来计算未剪接（ElogU）和已剪接（ElogS）的期望值。
最后，原始数据（S和U）通过使用期望值和从‘pyro.dist.Gamma’中抽取的噪声参数（shape_inv），从‘pyro.dist.GammaPoisson’中抽样得到。
本研究生成的所有模拟数据均可在Zenodo上获取（参见数据可用性）。
使用'velocycle.utils.simulate_data'函数可以模拟额外的数据集。

Para_03

使用了20个数据集来评估流形学习步骤，每个数据集包含3,000个细胞和300个基因，独立模拟的真实速度为0.4。
相同的数据库也用于验证速度学习步骤。
为了对细胞和基因数量进行敏感性分析，生成了四个独立模拟的数据集，每个数据集包含10,000个细胞和1,000个基因；使用数据子集测试模型在不同数量的细胞（流形学习从100到5,000，速度学习从50到10,000）和基因（流形学习从100到1,000，速度学习从50到1,000）上的性能。
为了评估速度学习在具有不同真实速度的数据集上的表现，我们模拟了四个具有共享动力学和基因谐波参数的数据集，但具有16种不同的真实速度，范围从0.0到1.5。

Para_04

使用‘velocycle.utils.circular_corrcoef’计算估计变量和模拟真实变量之间的循环相关性，该方法将输入数据转换为单位圆坐标，并通过计算估计值与真实值的复共轭的乘积的平均值来计算相关性。
为了比较 VeloCycle 获得的相位与 DeepCycle 获得的相位，使用相同的模拟数据集计算速度矩，然后按照原始出版物中描述的默认参数运行 DeepCycle。

VeloCycle estimation across multiple standard scRNA-seq datasets

跨多个标准单细胞RNA测序数据集的VeloCycle估计

Para_01

在这项工作中，我们使用 VeloCycle 在来自不同技术、物种和采样背景的多个已发表数据集上进行了流形几何和细胞周期速度估计。
对于所有数据集，原始原始数据均使用 VeloCycle2 重新处理以获得拼接和未拼接的计数矩阵。
上述已经描述了在这些类型的单细胞 RNA 测序数据上运行 VeloCycle 的一般程序，并由相应 GitHub 页面上的教程支持这些工作。
在这里，我们提供了一个摘要，概述了基于数据集的具体过滤标准和参数。

FACS-sorted mouse embryonic stem cells

FACS 分选的小鼠胚胎干细胞

Para_01

使用标准参数对来自 Smart-seq2 mES 细胞培养的 279 个单细胞进行了 VeloCycle 细胞周期阶段和基因谐波估计。
用于流形学习的基因是‘velocycle.utils’中可用的‘大’基因集（GO；1,918 个基因）中的基因，在过滤掉每个细胞平均未剪接计数≤0.5 或平均剪接计数≤2 的基因后（剩余 1,358 个基因）。
流形学习使用了 3,000 次训练步骤。

Para_02

为了评估使用 VeloCycle 阶段预测细胞周期分类（G1、S 或 G2/M）的能力，使用了来自 'sklearn.tree' 的 'DecisionTreeClassifier'，用 65% 的细胞进行训练，保留 35% 的细胞作为测试集和混淆矩阵的计算。
为了与使用总基因表达矩阵预测细胞周期分类的模型进行比较，使用了来自 'sklearn.linear_model' 的线性 'LogisticRegressionCV' 模型（c.v. = 5），并使用与决策树相同的训练-测试细胞分割进行了训练。

Mouse embryonic stem cells and human fibroblasts

小鼠胚胎干细胞和人成纤维细胞

Para_01

VeloCycle 使用标准参数分别在 5,191 个来自 mES 细胞培养的单细胞和 2,557 个来自人类成纤维细胞培养的单细胞上运行。
根据作者的注释，在分析前过滤掉了非循环细胞。
用于流形学习的基因是从 velocycle.utils 中可用的‘中等’基因集（DeepCycle；218 个基因）中选择的，过滤掉每个细胞平均未剪接 UMIs ≤0.1 或每个细胞平均已剪接 UMIs ≤0.3 的基因（mES 细胞和成纤维细胞分别剩余 189 个和 160 个基因）。
使用 5,000 步训练进行流形学习，并使用‘正常’指南和恒定速度模型进行 10,000 步训练进行速度学习。
与 DeepCycle 阶段的比较使用了原始研究中描述的这些确切数据集的已发布估计值。

Human dermal fibroblasts

人真皮成纤维细胞

Para_01

VeloCycle 在 1,222 个未处理的 dHF 单细胞培养物上运行，使用标准参数；非循环细胞根据作者的注释被排除。
用于流形学习的基因是 ‘velocycle.utils’ 中可用的 ‘大’ 基因集（GO；1,918 个基因），过滤掉每个细胞平均未剪接 UMI ≤0.1 或每个细胞平均剪接 UMI ≤0.3 的基因后，剩余 876 个基因。
流形学习使用 5,000 次训练步骤进行，速度学习使用常数速度模型和周期性速度模型进行 10,000 次训练步骤，使用 LRMN（‘lrmn’）引导。

Para_02

时间序列显微镜数据，包括细胞分割和追踪，对于 dHFs 来自最初发表的研究，并可在 Zenodo 上获取，网址为 https://doi.org/10.5281/zenodo.6245943。
如果一个细胞的测量细胞周期长度小于 8 小时或大于 32 小时，则认为该细胞跟踪不佳并从分析中排除。

PC9 lung adenocarcinoma cell line

PC9肺腺癌细胞系

Para_01

VeloCycle 使用标准参数在 PC9 肺腺癌细胞系接受厄洛替尼治疗前（D0，7,927 个细胞）和后（D3，3,743 个细胞）的数据上运行。
用于流形学习的基因是‘velocycle.utils’中可用的‘大’基因集（GO；1,918 个基因），在过滤掉每个细胞平均未剪接UMI≤0.1或每个细胞平均已剪接UMI≤0.1的基因后得到的。
在初始流形学习步骤之后，仅保留了未剪接和已剪接计数之间的皮尔逊相关性≥0.8且预测的未剪接-已剪接延迟大于≥-0.25的基因。
使用 5,000 次训练步骤进行流形学习，使用‘lrmn’引导和恒定速度模型及周期速度模型进行速度学习，训练步骤为 10,000 次。

Radial glial progenitors from the developing mouse brain

来自发育中鼠脑的放射状胶质祖细胞

Para_01

VeloCycle 在 E10 时间点的所有 RG 前体细胞上联合运行，根据区域身份进行分层（FB，3,293 个细胞；MB，2,388 个细胞；HB，2,012 个细胞），使用标准参数。
用于流形学习的基因是‘velocycle.utils’中可用的‘大’基因集（GO；1,918 个基因）中的基因，在过滤掉每细胞未剪接 UMI 平均值≤0.05 或每细胞已剪接 UMI 平均值≤0.1 的基因后。
在初始流形学习步骤之后，仅保留了未剪接和已剪接计数之间的皮尔逊相关性≥0.8 且预测的未剪接-已剪接延迟大于≥-0.10 的基因。
流形学习使用 5,000 次训练步骤进行，速度学习使用‘lrmn’引导和恒定速度模型进行 10,000 次训练步骤。

Para_02

同样地，VeloCycle 在所有来自 E14 和 E15 时间点的 RG 前体细胞上联合运行，根据区域身份（FB，2460 个细胞；MB，307 个细胞；HB，176 个细胞）进行分层，使用标准参数。
与上述 E10 分析相同的基因过滤步骤，使用了 239 个基因。
使用 5,000 步训练进行了流形学习，使用 'lrmn' 引导和恒定速度模型进行了 10,000 步训练的速度学习。

Para_03

为了在空间上可视化 E10 时间点的 VeloCycle 速度估计，我们运行了 BoneFight 算法，将 scRNA-seq 聚类映射到同一研究中的 HybISS 空间转录组学数据集，然后根据它们的速度估计对相应的聚类进行着色。

Genome-wide Perturb-seq RPE1 cells data

全基因组范围的扰动测序RPE1细胞数据

Para_01

为了确保分析仅在具有目标基因完全敲除的 RPE1 细胞上进行，我们过滤掉了含有非零未剪接或已剪接 UMI 读数的目标基因的细胞。
VeloCycle 最初在两种条件下的一组数据子集上运行：(1) 控制组 NT 细胞（11,485 个细胞）和一组针对‘小’细胞周期标志物列表中的基因进行敲除的细胞（CC-KO，6,275 个细胞）。
用于流形学习的基因是从‘velocycle.utils’中可用的‘中等’基因集（DeepCycle；218 个基因）中选择的，经过筛选去除了每个细胞平均未剪接 UMI ≤0.1 或每个细胞平均已剪接 UMI ≤0.2 的基因。
在初步的流形学习步骤之后，只保留了未剪接和已剪接计数之间的皮尔逊相关系数 ≥0.7 且预测的未剪接-已剪接延迟大于 ≥−0.5 的基因（剩余 120 个基因）。
流形学习使用 5,000 步训练完成，速度学习则使用‘lrmn’指导和恒定速度模型进行了 10,000 步训练。

Para_02

对条件无关的周期性傅里叶级数分量进行估计尤其具有挑战性，特别是在 Perturb-seq 敲低条件下，这些条件要么（1）细胞数量非常少，要么（2）细胞仅属于细胞周期的一个阶段。
因此，为了推断这些细胞的准确细胞周期阶段，我们首先进行了流形学习，经过 5,000 次训练步骤，以估计较大数量的 NT 和 CC-KO 细胞的基因谐波系数（ν0、ν1sin 和 ν1cos），这些细胞在细胞周期的各个阶段分布更为均匀。
之后，我们再次运行了流形学习，但这次是在包含 167,119 个细胞和 986 个基因敲低条件的整个 Perturb-seq 数据集上进行了 5,000 次训练步骤。
这一次，我们将 VeloCycle 条件化为第一步中学到的基因谐波系数。
这使得每个分层条件下的细胞可以被分配到细胞周期流形上的一个位置，但限制了这些分配，使其基于从更大和更有信息量的数据集中学到的基因表达模式（术语 Δν 允许批次效应表达差异）。
最后，我们在整个数据集上进行了 10,000 次训练步骤的速度学习，为每个基因敲低条件估计了一个单独的恒定速度。

Data availability

Para_01

为本研究新生成的 RPE1 细胞系的原始和处理过的单细胞 RNA 测序数据可在 Gene Expression Omnibus 存取号 GSE250148 处获得。
本研究中使用的其他所有单细胞 RNA 测序数据均来自先前发表的作品，并依赖于原作者提供的细胞类型注释。
这些作品在整个文章中适当引用。
用于重现本研究中所示结果的 Jupyter 笔记本和其他相关文件可通过我们的 GitHub 页面获取，网址为 https://github.com/lamanno-epfl/velocycle/。
所有已发布数据的处理版本（包括拼接-未拼接计数矩阵）也可在上述链接处获得。
模拟的单细胞 RNA 测序数据集、新的 RPE1 样本处理后的单细胞 RNA 测序元数据、活体显微镜下的细胞追踪数据和累积 EdU 染色实验数据也可在上述链接处获得。

Code availability

Para_01

VeloCycle 是用 Python 实现的，并作为开源包在 GitHub 上提供，网址为 https://github.com/lamanno-epfl/velocycle。
可以使用命令 pip install VeloCycle 从 PyPi 安装 VeloCycle，或者通过 GitHub 页面直接安装，使用命令 pip install ‘git+https://github.com/lamanno-epfl/velocycle.git@main’。
需要 Python 3.8 或更高版本。
源代码、安装说明、教程和包含所有必需包依赖项的文件也在 GitHub 上提供。
用于重现本研究结果的附加代码和笔记本可在 Zenodo 上获得，网址为 https://doi.org/10.5281/zenodo.12517650

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-01-02，如有侵权请联系 cloudcommunity@tencent.com 删除

程序