RNA速率实战：用scVelo倒放单细胞命运轨迹

干这行久了，处理过的单细胞数据越多，越会觉得手里攥着的东西有点“不对劲”。我们花了大价钱做实验、跑测序，最后拿到的其实是一张静态快照：细胞死了，RNA被抽出来了，那一刻的转录状态被定格。但细胞是活的，它一直在变化，从一个状态走向另一个状态。尤其是当你盯着UMAP图，看到一簇未分化的细胞和一簇终末分化的细胞离得很远，中间连着一长串过渡态的时候，心里那个问题就绕不开了：它们到底是怎么走过去的？中间经历了什么？有没有办法让这些点“动起来”，甚至把命运演化的过程“倒放”一遍？RNA速率（RNA velocity）就是我目前找到的最接近这个答案的工具。这篇文章我会把它的原理、实操、以及“倒放”细胞命运的具体玩法，完整地拆开讲一遍。

1. 项目的缘起：为什么我们需要“倒放”细胞命运

做单细胞转录组（scRNA-seq）分析时间久了，你会发现一个很有意思的矛盾：我们最想研究的恰恰是数据里最缺失的那个维度——时间。测序的那一刻，每个细胞都被固定在了它生命历程中的一个特定转录状态上，但我们对它的过去和未来一无所知。伪时间分析（如Monocle、PAGA）能根据表达相似性给细胞排个序，但那本质上是在“静态数据里找渐变规律”，排序的方向需要人为指定，带有很强的先验假设。而RNA速率不一样，它是在单个细胞内直接读取“时间箭头”。

这个思路的起点相当朴素：细胞表达基因的时候，并不是一下子就把成熟的mRNA丢到细胞质里。新转录出来的RNA先是带着内含子的前体（pre-mRNA），经过剪接加工变成成熟的mRNA，最后再被降解。所以，同一个基因在同一时刻会存在“未剪接”和“已剪接”两种状态。如果某个基因的未剪接mRNA占比很高、已剪接mRNA相对偏少，说明这个基因正在加速转录，生产线刚开起来；反过来，如果未剪接的库存很少，已剪接的mRNA却很多，说明这条基因的表达在往下走，正在停产清库存。把几千个基因的这种“生产状态”综合起来，就能判断这个细胞此刻正在往哪个方向移动。

这正是“倒放录像”的第一个关键点：普通表达矩阵丢掉的信息，在RNA速率里被重新捡了回来。用一个不严谨但很好懂的说法——传统单细胞数据是电影的某一帧，而RNA速率给出了这一帧里每个像素的“速度矢量”。知道了每个细胞的速度方向，你就能顺着时间轴正推，看它下一步会变成什么；也可以逆着方向回溯，追查一个分化终末态的细胞是从哪条路走过来的。

这篇文章适合所有正在做单细胞数据分析、手里有scRNA-seq数据但只做了聚类和注释的朋友。如果你只做过Seurat或Scanpy的标准流程，这篇文章能帮你把RNA速率的原理和实操串起来；如果你已经试过velocyto或scVelo但结果不理想，后半部分的排查技巧和“倒放”玩法应该会对你有直接帮助。当然，RNA速率不是万能钥匙，它有很多坑，很多结果需要谨慎解读。我的目标是用尽量通俗的语言，把这些坑和经验一并讲清楚。

2. 核心原理拆解：RNA速率到底在算什么

2.1 剪接状态：细胞给自己打的“时间标记”

理解RNA速率，绕不开一个分子生物学基础：基因表达不是瞬间完成的，它是一条流水线——转录、剪接、降解。在单细胞测序的文库制备中，如果你比对到的是含内含子的reads（或者通过其他方式识别未剪接reads的来源），就可以把它归为“未剪接”（unspliced）；比对到全长成熟转录本的，就是“已剪接”（spliced）。于是每个基因在每个细胞中有了两个计数，而这两个计数的比例关系，就隐含了表达趋势。

我把这个过程类比成一家工厂的生产线。未剪接mRNA是刚下线、还没包装的半成品，已剪接mRNA是贴上标签准备发货的成品。如果半成品堆积如山、成品区空空荡荡，那说明这台机器刚刚开足马力，很快就会有一大批成品涌出来——这个基因的表达量正在上升。反之，半成品寥寥无几、成品越积越多，说明机器已经停了，接下来成品会被降解清掉——表达量正在下降。而在稳态下，流水线的产量恰好等于出货加损耗，半成品和成品比例基本恒定。

这个类比里的“半成品率”就是RNA速率计算的原料。现实中我们并不真的关心每个基因的绝对转录速率，而是关心“相对趋势”。只要知道每个基因在当前表达水平下的剪接状态是偏高还是偏低，就能推断这个基因是往上调还是往下调。所有基因综合起来，就构成了细胞的转录速度向量。

很多初学者会忽略一个重要细节：想要区分未剪接和已剪接reads，在测序数据的预处理阶段就得做好准备。普通的STAR比对模式会丢掉或混淆内含子信息，你需要用专门的方式（比如10x数据在Cell Ranger内置的splice-aware比对，或velocyto的run10x等流程）生成一个包含未剪接/已剪接计数的矩阵。这决定了后续分析的原料质量，如果这一步没做好，后面怎么调参都是事倍功半。

2.2 稳态模型与动态模型：两种“算速度”的路径

有了未剪接/已剪接计数，接下来怎么算速度？目前主流有两条技术路线。

第一条是早期velocyto提出来的稳态模型（steady-state model），也被称为“残差模型”。它的逻辑是：对于一个基因，在长时间量纲下，细胞会达到一个剪接动力学稳态，此时已剪接mRNA数量（s）与未剪接mRNA数量（u）之间存在一个近似线性的关系，斜率由剪接速率与降解速率的比值决定。用所有细胞的数据拟合这条线性回归线，然后看每个细胞落在这条线的哪一侧：如果某个细胞的u相对s偏高，说明它高于稳态线，即正在激活这个基因；反之则在抑制。

稳态模型简单、计算快，但有一个致命前提：所有细胞必须处于稳态。这在真实组织中经常不成立，比如发育早期、应激状态、细胞周期活跃的组织，很多基因根本没到稳态就被扰动打断了。结果就是拟合线极不靠谱，速度方向全被带偏。这也是后来动态模型出现的原因。

第二条是scVelo引入的动态模型（dynamical model），核心思路是不再假设稳态，而是直接对每个基因的剪接动力学全过程建模。它用一组包含转录、剪接、降解过程的微分方程描述基因表达随时间的变化，然后通过期望最大化（EM）算法去拟合每个基因的动力学参数。

这两条路径的核心差异可以用一句话概括：稳态模型是在“拍一张快照，看谁偏离了平均线”；动态模型是在“完整建模从开启到关闭的过程，找到每个细胞处于过程中的哪个时间点”。动态模型需要估算的参数更多、计算量更大，但能同时给出“基因阶段”和“细胞潜在时间”（latent time），对复杂组织数据的适应性远比稳态模型好。实际操作中，除非你的样本是纯化得很干净的稳定细胞系或简单细胞群体，否则我建议直接跳过稳态模型，用动态模式。

2.3 关键公式与参数直觉：从表达量到速度向量

理论上，RNA速率的核心动力学可以用一个很朴素的常微分方程组来描述。为了不用吓跑读者，我用最简化的形式说明：

对于单个基因，未剪接RNA的数量u和已剪接RNA的数量s随时间的变化为：

du/dt = α - β·u

ds/dt = β·u - γ·s

这里α是转录速率，β是剪接速率，γ是降解速率。当系统达到稳态时，du/dt和ds/dt都为零，可以得到s/u = β/γ（稳态比）。但动态模型不会停留在稳态，它会拟合α、β、γ这些参数，然后把每个细胞（u, s）坐标映射到整个动力学轨迹上，算出该细胞瞬间的“净生产速率”。

这里有一个特别重要的直觉：RNA速率不是简单地等于“表达量的变化率”。因为同一时刻测得的u和s是动态过程的一个快照，RNA速率向量给出的是“细胞在当前转录状态下，沿着分化轨迹移动的倾向”。它描述的是一个在分化轨迹上运动的细胞的瞬时速度，而不是轨迹本身。因此，速率向量里的每个分量都有明确的生物学含义：正值表示该基因表达正在上升，负值表示正在下降。把所有基因的分量汇集起来做降维展示，就得到了在UMAP或t-SNE上叠加的“箭头流场”。

我建议在实操中不要把注意力全放在漂亮的箭头图上，先去检查几个关键指标：模型是否收敛（scVelo会输出收敛信息）、latent time是否沿着分化轨迹单调递增、每个基因的动力学拟合是否合理。参数不是越多越好，动态模型的很多参数需要根据数据规模调整，比如n_top_genes（用于拟合速度的高变基因数）和min_shared_counts（基因在最小细胞数中的表达量阈值）。过于激进的过滤会丢掉稀有细胞类型中的关键基因，过滤不足则会引入大量低表达基因的噪声，导致速度计算被“嘈杂基因”主导。初次跑通流程后，我每次都会重新审视这两个参数，它们对结果的影响往往比很多人想象的更大。

3. 实操全流程：从原始数据到细胞命运录像

3.1 输入数据：你需要一张“带时间戳”的表达矩阵

RNA速率分析的前提是和普通单细胞分析不一样的输入：你必须同时拿到已剪接和未剪接两个版本的标准表达矩阵（spliced和unspliced counts）。如果是10x Genomics的Chromium数据，velocyto提供了专门的脚本（比如velocyto run10x）来处理Cell Ranger的BAM文件，比对到参考基因组后区分出两类reads。如果用的是dropEst、STARsolo等工具，也能输出包含未剪接/已剪接信息的结果。

拿到数据之后，我习惯先把两种counts合并成AnnData对象，并保留原始计数。有一个很容易踩的坑：Scanpy的很多预处理函数默认会对X矩阵做归一化，但RNA速率分析用的是单独的层（layers），千万不要在spliced/unspliced上随便跑total counts归一化。你可以在之后用scVelo自带的过滤和归一化函数，它会同时处理两层数据并保持它们的比例关系。手动混用两套预处理逻辑的结果基本就是：速度方向全部扭曲，箭头图看着花哨但毫无生物学意义。

代码层面的准备工作大致是这样：

# 通过velocyto生成loom文件（示例命令，具体参数按实际环境调整） velocyto run10x -m repeat_mask.gtf \ -o ./velocyto_output \ ./cellranger_output/possorted_genome_bam.bam \ ./cellranger_output/genes.gtf

或者直接用已经做好的h5ad文件（包含spliced和unspliced层）：

import scanpy as sc import scvelo as scv adata = sc.read_h5ad("data_with_velocity.h5ad") # 查看数据结构 print(adata) print(adata.layers.keys()) # 应该包括 'spliced' 和 'unspliced'

如果不能生成包含spliced/unspliced的数据，后续的所有分析都无法进行。这是整个RNA速率分析的第一个硬性门槛。

3.2 预处理与质量控制：比普通scRNA-seq更挑剔

RNA速率分析的预处理整体上沿用Scanpy的框架，但细节上比标准流程更挑剔。我的习惯做法是：

先做一个基础的细胞过滤（比如低质量细胞、双细胞），然后在AnnData上运行scVelo的过滤和归一化：

scv.pp.filter_and_normalize( adata, min_shared_counts=20, n_top_genes=2000, min_counts=2, min_cells=5 ) scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

这里面有几个细节值得展开。

第一，min_shared_counts这个参数决定了一个基因至少要在多少个细胞中同时检测到未剪接和已剪接信号。设得太低（比如5），大量基因本身表达量极低、剪接比例不稳定，会把噪声带进后续拟合；设得太高（比如50以上），会丢失掉只在少数细胞中特异性表达的谱系关键基因。我一般从20开始，根据数据稀疏程度上下调整。如果数据测序深度很深（每细胞UMI数很高），可以适当提高；如果数据很稀疏，反而要降一点。

第二，n_top_genes控制用于速度估计的高变基因数量。默认2000听起来够用，但在异质性很强的组织（比如胚胎组织、肿瘤微环境）里，我会提升到3000以上。RNA速率的方向是由大量基因的“投票”决定的，高变基因太少会丧失对弱信号的捕捉能力。但也不要无脑加到5000——动态模型对每个高变基因都要做一轮动力学拟合，基因数越多计算时间越长，而且边际收益会递减。

第三，moments函数本质上是基于邻域图对表达量做平滑，降低单细胞噪声。n_neighbors和n_pcs在这里并不只是可视化的参数，它们会直接影响速度向量的计算。尤其是n_pcs，降维保留的主成分太少会损失稀有细胞群体的信息，太多又会把技术噪声引进来。对于几千个细胞的小数据集，n_pcs取15-20足够；几万个细胞的大数据集，30-50都可以试。我通常会跑一个小网格搜索，用速度一致性（速度向量与转录方向的一致性程度）作为粗略的评估指标。

预处理阶段还有一个经常被忽视的问题：批次效应。如果样本来自多个批次、多个样本，一定要在RNA速率分析前先做批次整合（比如Harmony、scVI等），否则速度场会被批次差异完全主导。我在一个模拟项目X中遇到过这样的情况：两个批次的细胞各自朝着自己的方向“跑”，箭头图看起来像是两股对吹的风——问题根本不是生物学差异，而是批次差异。先整合、再算速度，结果才恢复正常。

3.3 速率估计：用动态模型拟合每个基因的“生命周期”

预处理完成之后，核心的计算就在scVelo里。我强烈建议直接用动态模型而不是默认的稳态模型，尤其是当你研究的是发育过程、组织再生、肿瘤异质性这类高度动态的系统时。具体操作：

# 恢复动力学（拟合每个基因的转录、剪接、降解参数） scv.tl.recover_dynamics(adata, n_jobs=8) # 用恢复的参数计算RNA速率 scv.tl.velocity(adata, mode="dynamical") # 计算速度图（用于向量场和转移矩阵） scv.tl.velocity_graph(adata)

recover_dynamics这一步是动态模型的核心，它会为每个基因独立拟合出α、β、γ参数。n_jobs可以按CPU核心数调整，这一步是最耗时的，但值得等。拟合完成后，scVelo会为每个基因给出一个“基因阶段排序”（gene expression phase），本质上就是该基因在分化过程中开启到关闭的时间顺序。这个输出非常有用，后面讲“倒放”玩法时我们还会用到它。

需要注意scVelo会输出大量warnings，别被它们吓到。常见的一种是“基因在极少细胞中表达，跳过动力学拟合”之类——这通常不是错误，只是提醒你某些基因没有足够的信息构建可靠动力学模型。对应的基因会被自动排除在速度计算之外，这是一种合理的保守策略。

另一种需要警惕的情况是：recover_dynamics的拟合过程不收敛。导致不收敛的原因很多，最常见的是数据稀疏度过高，或者基因本身缺乏明显的剪接动态模式。如果大量基因（比如超过20%）拟合失败，不要急着调参，先回去检查输入数据——是否有批次效应、是否包含了过多的低质量细胞、是否在预处理阶段误用了对spliced/unspliced比例有影响的归一化步骤。

3.4 向量场可视化与轨迹推断

有了速度向量，接下来就可以把它投影到UMAP或其他降维坐标上，得到经典的“流线图”和“箭头图”。

# 将速率向量投影到UMAP坐标上 scv.tl.velocity_embedding(adata, basis="umap") # 画流线图（stream图，看起来就像一条顺着分化方向流动的河） scv.pl.velocity_embedding_stream(adata, basis="umap", color="cell_type") # 画箭头图（每个细胞一个箭头，方向直观但噪声较大） scv.pl.velocity_embedding_grid(adata, basis="umap", color="cell_type")

看流线图的时候，有几个判断技巧。

如果结果理想，你会看到箭头流场大致沿着细胞分化的方向——比如从干/祖细胞群指向终末分化群，而且它们之间没有明显的“对流”或“对撞”。如果流场的箭头毫无方向感，像一团乱毛线，要么是数据质量问题，要么是这个组织里同时存在大量互相独立的分化过程，UMAP把它们压在一起后箭头彼此抵消。

我还习惯同时画一下latent time：

scv.pl.scatter(adata, color="latent_time", basis="umap")

latent time就是动态模型推断出的每个细胞“内部时钟”读数。在理想的数据中，latent time应该沿着分化轨迹单调变化。如果你发现latent time在某个细胞群里突变或跳跃，说明模型的动力学假设和该细胞群的真实生物学状态不匹配——比如细胞周期效应混了进来，或者这个群里核心的驱动基因没有被选入高变基因。

可视化之后，另一个常用的下游分析是构建细胞状态间的转移关系。scVelo自带的velocity_graph本质上已经用速度向量构建了细胞之间的马尔可夫转移矩阵。基于这个矩阵，你可以做多类分析：从指定起点模拟分化“流向”，鉴定命运转归概率，甚至回溯细胞来源。这就是下一节“倒放”玩法的基础。

4. “倒放”的三种玩法：追溯、排序、预测

标题里提到的“倒放”，对很多人来说可能只理解成一个比喻——但实际操作上，它确实意味着可以把时间方向反过来用。下面是我在真实分析中使用最多的三种“倒放”思路。

4.1 回溯祖细胞状态：从终末细胞反推来路

最常见的“倒放”需求是：我已经知道某群细胞是终末分化状态（比如成熟的胰岛β细胞、终末分化的神经元），我想知道它们是从哪些祖细胞来的、经过了哪些中间状态。单细胞聚类能给你一群形态相近的细胞，但聚类本身不会告诉你它们的“来路”。

这时候，速度向量场可以反向使用：把每个细胞看作流场中的一个点，把时间方向取反，让“水流”从终点倒流回起点。实际操作上，基于velocity_graph转移矩阵，可以反复执行“反向转移”，找到终末细胞群在速度场上最可能的来源区域。

如果你手头有CellRank，它会把这个过程做得更系统。CellRank的核心思路是：用RNA速率信息构造一个细胞间的转移核（transition kernel），然后它不仅能算出从任意起点细胞出发最终落到哪个“吸收态”（命运终点）的概率，还能算出“反向吸收概率”——即任意细胞在过去最可能来自哪个起始状态。这个反向概率几乎就是“倒放细胞命运”的数学化身。

4.2 转录因子层级的时间排序

第二种“倒放”玩法更精细：追溯发育过程中转录调控网络的建立顺序。发育生物学里有一类问题特别常见：一个转录因子家族里好几个成员都在终末细胞中高表达，但谁先启动、谁后启动，哪个是上游的“开启者”，哪个只是下游的“执行者”？普通表达数据只能看到最终结果，而RNA速率能给出每个基因在表达过程中的“激活-关闭”相位，配合scVelo的输出就能推断调控层级。

具体做法是：先选定一组感兴趣的转录因子（TF），看它们在动态模型拟合结果中的基因相位阶段，以及它们在速率向量中的方向性。如果一个TF在某群细胞中未剪接信号很强，意味着它正在被快速激活，很可能处于上游启动阶段；另一个TF的在相同区域中已剪接信号为主但未剪接信号较少，说明它已经开始关闭，可能在时间轴上更靠前或更靠后，根据具体细胞群判断。

更进一步，可以计算TF与其靶基因之间的速度相关性：TF的转录激活状态（由RNA速率近似）与靶基因表达上升是否存在时序上的先后关系。如果TF的未剪接比例上升之后，靶基因的已剪接比例在latent time轴上跟着上升，这就是一个典型的“先启动后执行”的级联证据。我在处理一个模拟的免疫细胞分化数据时，就靠这个思路把一群候选TF分成了“早期启动者”和“晚期维持者”两类，后续用实验数据验证了其中几个关键节点的时序推断，准确度很高。

4.3 命运转归概率与细胞可塑性

第三种玩法可以视为“对未来的预测”——也就是正放录像。但这其实也间接回答了一个“倒放”问题：一个细胞如果还没确定命运，它有多大概率走向A或B？速度向量场可以量化这个“转归概率”。

使用CellRank时，可以设置一个起点细胞群和若干终点细胞群，然后计算每个起点细胞进入每个终点的概率分布。这会给你一个非常直观的产出：细胞命运偏好图，或者一个“转归概率矩阵”。矩阵里每个非零元素都代表：来自起始群i的细胞有多大比例最终进入终点群j。

这个分析的一个关键注意点是：转归概率是“基于当前转录状态的趋势预测”，不是终审判决。细胞是高度可塑的，尤其是受到外界信号扰动时，转录状态会迅速改变，概率也会随之变化。所以我把转归概率理解为“用转录组快照做的最优推断”，而不是宿命论。它最有效的应用场景其实是找到“犹豫不决”的细胞——那些同时具有高概率进入两个终端的细胞，往往就是最有趣的可塑性细胞或处于关键决策点的细胞。这类细胞往往在普通聚类里会被归为某一群，但转归概率矩阵会告诉你，它们的命运其实尚未锁死。

5. 常见问题与排查技巧实录

RNA速率分析的结果不好，大部分时候不是方法本身不行，而是踩了某些隐蔽的坑。下面按我实际遇到的频率排序，整理成了一张速查表：

除了表格里的这些场景，还有一个我从踩坑中总结出来的习惯：永远不要在拿到结果的第一时间就下结论。我的流程是，先观察全局流线图，再单独检查关键谱系的latent time分布和TF动态，然后才回到全局看转归概率。每一步都要做生物学意义上的交叉验证——比如RNA速率指向某个中间态是某种终末细胞的前体，那么目标终末细胞的特异性marker在这个中间态里应该已经出现上升趋势（未剪接信号增强）。如果连这个基本的一致性都验证不了，说明分析管线里一定有某个环节出了问题。

另外一个容易被忽略的问题是技术平台差异。10x平台的未剪接reads占比通常并不高，尤其在测序深度不够的情况下，未剪接信号的稀疏性非常严重。Smart-seq2这类全长转录组平台虽然能更好地覆盖内含子，但细胞数一般较少，邻域图构建和速度向量计算容易不稳定。处理不同平台数据时，预处理参数和模型选择都不能照搬。

6. 我在反复实操中最深的体会

RNA速率这个工具，说到底是把细胞内部的分子记忆变成了一张动态地图。它对直觉的依赖比传统聚类分析高得多——因为它输出的不是“分组”，而是“方向”。当你第一次看到箭头流场把一个未分化细胞群平滑地导向两个终末分化群时，那种感觉确实很像给静止的细胞命运按下了播放键；当你把转移矩阵反过来跑一遍，看到终末细胞群向几个不同的祖细胞状态发散时，也就真的体会到了什么叫“命运录像倒放”。

但我也要说几句泼冷水的话。

RNA速率的结果本质上仍是转录层面的统计学推断，不是直接观测到的细胞命运。它受数据质量、参考基因组注释、剪接位点识别效率的影响非常大，尤其是未剪接reads的识别，在不同比对工具之间可能差异明显。这意味着不同流程跑出来同一个数据集的“录像”可能不完全一样，如果发现这种不一致，不要急着判断哪个是对的，而是去检查最底层的数据比对和计数矩阵，差异往往在源头。

另一个体会是：RNA速率最好和多种分析手段配合使用，而不是孤军奋战。用Monocle或PAGA确定轨迹主干结构，用RNA速率确定方向并细剖析分支点，用转录因子动态推断调控层级，再结合实验验证才能形成完整证据链。我见过不少报告，单看RNA速率的流线图很漂亮，但细看之下PAGA的拓扑和速度箭头完全冲突——这种时候，更应该相信多组学的交叉共识，而不是某一个软件给出的“精彩画面”。

最后还是那句老话：任何生信工具的结果都是“假说生成器”，RNA速率尤其如此。它能高效地生成关于细胞命运的有根据的假说，但真正的命运，最终还是要靠湿实验去确认。这套“倒放”的分析流程，最好的用法是帮你把有限的实验资源集中到最关键的细胞群体和决策时间点上。