python单细胞流程-RNA velocity

为什么要进行 RNA velocity 分析？

单细胞 RNA-seq 会在测序时破坏细胞，因此只能得到每个细胞在某一时刻的转录组快照，不能连续追踪同一个细胞之后会发生什么变化。

伪时序分析可以根据细胞之间的表达相似性推断相对顺序，但通常需要人为指定起点，也不能直接提供转录状态变化的方向。RNA velocity 进一步利用 unspliced RNA 和 spliced RNA 之间的关系，估计细胞当前转录状态的局部变化方向。

简单理解就是：

unspliced RNA 的变化先出现 → spliced RNA 的变化随后出现 → 利用二者之间的时间差估计表达变化方向

RNA velocity 适合尝试回答：

• 某段连续状态变化大致朝哪个方向发展。
• 哪些细胞更接近起始状态或终末状态。
• RNA velocity 推断的方向是否支持已有的伪时序结果。
• 哪些基因表现出较清晰的诱导或抑制动力学。

它不能直接证明：

• 一个细胞真实变成了另一个细胞。
• 两个细胞之间存在真实祖先-后代关系。
• UMAP 上箭头指向的 cluster 一定是未来细胞类型。
• 某个 velocity gene 是驱动分化的因果基因。

RNA velocity 仍然是基于 snapshot 数据和模型假设得到的推断，不能替代 lineage tracing、真实时间序列或活细胞追踪实验。

RNA velocity 适合什么数据？

RNA velocity 最适合存在活跃、连续并且被充分采样的状态转换过程，例如胚胎发育、祖细胞分化、组织再生和短时间尺度的刺激响应。

使用前要重点判断两个条件。

第一个条件是生物过程的时间尺度。目标过程的变化速度应当与 RNA 分子的剪接和降解时间尺度大致相当。胰腺内分泌细胞发生适合用于 RNA velocity；持续数月或数年的慢性疾病进程通常不能直接用当前的 splicing-based RNA velocity 解释。

第二个条件是数据中是否真的存在状态转换。成熟 PBMC 数据如果主要由稳定的 T cells、B cells、NK cells 和 monocytes 组成，并且这些细胞之间没有正在发生的转换，就不适合强行计算 RNA velocity。

不适合直接做 RNA velocity 的常见情况包括：

• 对象中没有原始 spliced 和 unspliced count。
• 样本主要由彼此无直接转换关系的成熟细胞类型组成。
• 数据只捕获到起点和终点，中间状态严重缺失。
• 状态变化主要由批次、测序深度或低质量细胞驱动。
• 同一个基因在不同 lineage 中存在明显不同的 kinetics，却被放在一起拟合。

如果一个基因在不同 lineage 中表现出多套明显动力学，应先按生物学相关的 lineage 进行 subset，再分别分析。

RNA velocity 和伪时序有什么区别？

伪时序主要根据细胞在表达空间中的位置和连通关系计算相对顺序；RNA velocity 则利用 unspliced 和 spliced RNA 估计局部变化方向。

• Pseudotime：细胞在轨迹上的相对位置。
• RNA velocity：细胞当前转录状态的局部变化方向。

二者可以互相验证，但不能互相替代。如果伪时序和 RNA velocity 的方向冲突，应检查数据、模型假设、邻居图和生物学先验，而不是选择图更好看的结果。

RNA velocity 的基本原理

DNA 转录后先产生包含 intron 的 unspliced pre-mRNA，经过剪接后形成 spliced mature mRNA。对某个基因，可以用下面的简化模型描述这个过程：

其中：

•  表示基因  的 unspliced RNA abundance。
•  表示基因  的 spliced RNA abundance。
•  表示转录速率。
•  表示剪接速率。
•  表示成熟 RNA 的降解速率。
•  表示 spliced RNA 的变化，也就是基因  的 RNA velocity。

当一个基因刚开始被诱导时，unspliced RNA 通常先增加；当转录被抑制后，unspliced RNA 通常先下降，spliced RNA 再滞后下降。

参数推断模型

根据上面的公式，如果剪接速率和降解速率已知，就可以计算 RNA velocity。但是，单细胞 RNA-seq 只能提供不同细胞在未知时间点的 snapshot：我们测到了每个细胞的 unspliced RNA 和 spliced RNA，却不知道细胞所处的真实时间，也不知道基因的动力学参数。因此，RNA velocity 不能直接从表达量中读出来，必须通过模型进行推断。

对于每个基因，可以把每个细胞的 unspliced 和 spliced abundance 组成一个  点，所有细胞的点共同构成 phase portrait。理想情况下，phase portrait 呈部分或完整的杏仁形：上部主要对应基因诱导，下部主要对应基因抑制。真实数据比较稀疏且噪声较大，因此 scVelo 会先根据邻居图对数据进行局部平滑，再利用 steady-state model 或 EM model 推断动力学关系。

了解模型的目的不是推导复杂公式，而是判断  velocity 是怎样算出来的、应该选择哪种计算方法，以及结果是否可信。不同模型使用不同假设，可能得到不同甚至相反的方向；如果模型假设不符合数据，即使 UMAP 上的流线很清晰，也不能直接进行生物学解释。

常见模型怎么选择？

Steady-state model 和 EM model 是两种可选的 velocity 推断方法，不是必须依次完成的两个步骤。二者可以单独运行，也可以分别运行后比较结果。

Steady-state model：快速估计 RNA velocity

Steady-state model 先根据 phase portrait 中接近稳态的细胞估计 unspliced/spliced steady-state ratio，再计算每个细胞偏离该关系的程度：位于稳态关系上方通常表示基因正在被诱导，位于下方通常表示基因正在被抑制。

下面这行代码会计算每个细胞、每个基因的 RNA velocity，但不会直接生成 UMAP 上的流线或箭头：

scv.tl.velocity(adata, mode="deterministic")

该模型计算快，适合初步检查 velocity 方向。它要求数据中观察到接近稳态的细胞，并假设动力学速率基本保持不变。如果实验没有覆盖足够完整的诱导或抑制过程，推断方向可能错误。

EM model：拟合完整的剪接动力学

scVelo 的 dynamical model 使用 expectation-maximization（EM）框架，为每个基因拟合转录速率、剪接速率和降解速率，同时估计细胞处于诱导、抑制还是稳态阶段。

下面两行代码分别完成两个任务：recover_dynamics() 先拟合动力学参数，velocity(mode="dynamical") 再利用这些参数计算 RNA velocity。

scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode="dynamical")

EM model 不要求实验中一定观察到 steady state，因此更适合连续发育中的 transient cell states，也可以进一步计算 latent time。它的计算量更大，并且仍然要求 phase portrait 大致符合预期动力学；模型更复杂，不代表结果一定更正确。

简单选择原则：

• 只想快速检查 RNA velocity 的大致方向：先使用 steady-state model。
• 研究连续发育或其他 transient states，并且需要更完整的动力学结果：使用 EM model。
• 如果基因没有合理的 phase portrait，或者结果明显违背生物学先验，两种模型的结果都不应直接解释。

RNA velocity 的基本流程

以下使用胰腺内分泌细胞发生数据说明分析流程：

准备原始 spliced/unspliced count → 过滤和归一化 → PCA 和邻居图 → 计算 moments → 选择 steady-state 或 EM model（也可同时计算用于比较）→ 检查 phase portrait 和低维投影

加载库

import os
import scanpy as sc
import scvelo as scv

设置绘图参数：

sc.settings.verbosity = 3
scv.settings.verbosity = 3
scv.settings.set_figure_params("scvelo")

读取示例数据

这里使用 scVelo 内置的 pancreas 数据。该数据包含胰腺内分泌发育过程，并且已经提供 spliced 和  unspliced layers。

adata = scv.datasets.pancreas()
adata

指定细胞类型列：

cell_type_key = "clusters"

检查对象：

print(adata.obs.columns)
print(adata.layers.keys())
print(adata.obsm.keys())
print(adata.obsp.keys())

查看 spliced 和 unspliced 的总体比例：

scv.pl.proportions(
    adata,
    groupby=cell_type_key
)

不同实验平台和组织中的 unspliced 比例可以不同，不能机械使用固定比例作为合格标准。还需要结合 count 数量、phase portrait 和样本一致性判断。

选择相关 lineage

如果数据包含多个无关细胞类型，建议先选择一个生物学上合理的 lineage。下面以 pancreas 数据为例，需根据实际细胞注释修改细胞类型。

lineage_cell_types = [
    "Ductal",
    "Ngn3 low EP",
    "Ngn3 high EP",
    "Pre-endocrine",
    "Alpha",
    "Beta",
    "Delta"


    
,
    "Epsilon",
]

available_cell_types = set(adata.obs[cell_type_key].astype(str))
lineage_cell_types = [
    cell_type
    for cell_type in lineage_cell_types
    if cell_type in available_cell_types
]

if len(lineage_cell_types) == 0:
    raise ValueError("没有找到指定的 lineage cell types，请修改 lineage_cell_types。")

adata = adata[
    adata.obs[cell_type_key].astype(str).isin(lineage_cell_types)
].copy()

如果数据本身已经只包含目标 lineage，可以跳过这一步。

数据预处理

RNA velocity 应从未处理的 whole counts 开始。首先过滤在 spliced 和 unspliced RNA 中共同检测不足的基因，然后进行 cell-size normalization、对 adata.X 做 log1p，并选择高变基因。

先过滤低表达基因并进行归一化：

scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20
)

对 adata.X 进行 log1p 转换：

sc.pp.log1p(adata)

选择 2000 个高变基因：

sc.pp.highly_variable_genes(
    adata,
    n_top_genes=


    
2000,
    subset=True
)

重新计算 PCA 和邻居图：

sc.tl.pca(
    adata,
    n_comps=50,
    svd_solver="arpack",
    random_state=0
)

sc.pp.neighbors(
    adata,
    n_pcs=30,
    n_neighbors=30,
    random_state=0
)

如果前面进行了 subset，建议重新计算 UMAP：

sc.tl.umap(
    adata,
    random_state=0
)

根据现有邻居图计算 moments：

scv.pp.moments(
    adata,
    n_pcs=None,
    n_neighbors=None
)

moments() 会根据邻居图对 spliced 和 unspliced abundance 进行平滑，并生成：

• adata.layers["Ms"]：spliced moments。
• adata.layers["Mu"]：unspliced moments。

查看预处理后的细胞分布：

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

使用 steady-state model

首先计算 deterministic RNA velocity：

scv.tl.velocity(
    adata,
    mode="deterministic"
)

构建 velocity graph 并投影到 UMAP：

scv.tl.velocity_graph(
    adata,
    n_jobs=8
)

scv.pl.velocity_embedding_stream(
    adata,
    basis="umap",
    color=cell_type_key
)

这张图适合观察整体趋势。高维 velocity vectors 投影到二维空间后，结果会受到使用的基因、UMAP 结构和绘图参数影响，尤其不应过度解释 embedding 边界处的流线。

使用 EM model

先拟合每个基因的 splicing kinetics：

scv.tl.recover_dynamics(
    adata,
    n_jobs=8
)

拟合参数会写入 adata.var，主要包括 fit_alpha、fit_beta、fit_gamma 和 fit_likelihood。

根据 likelihood 选择拟合较好的基因：

top_genes = (
    adata.var["fit_likelihood"]
    .sort_values(ascending=False)
    .index[:5]
)

查看 phase portrait 和拟合轨迹：

scv.pl.scatter(
    adata,
    basis=top_genes,
    color=cell_type_key,
    frameon=False
)

理想 phase portrait 应呈现部分或完整的杏仁形，并能看到较清晰的诱导或抑制过程。如果相图只是一团噪声、出现多套明显 kinetics，或主要由少量离群细胞决定，就不应继续强解释该基因的 velocity。

计算 dynamical RNA velocity：

scv.tl.velocity(
    adata,
    mode="dynamical"
)

scv.tl.velocity_graph(
    adata,
    n_jobs=8
)

scv.pl.velocity_embedding_stream(
    adata,
    basis="umap",
    color=cell_type_key
)

EM model 会覆盖前面保存在 adata.layers["velocity"] 和 adata.uns["velocity_graph"] 中的 steady-state 结果。因此，如果需要比较两张图，应在运行 EM model 前先保存 steady-state 图。

如果两个模型给出的主要方向明显冲突，应回到 phase portrait、lineage subset、邻居图和已知生物学进行检查，不能默认更复杂的模型一定正确。

检查 velocity confidence

可以用邻域中 velocity vectors 的一致性进行辅助检查：

scv.tl.velocity_confidence(adata)

scv.pl.scatter(
    adata,
    basis="umap",
    color=["velocity_confidence", "velocity_length"]
)

velocity_confidence 较低的区域应谨慎解释。velocity_length 是表达空间中的向量长度，不是经过真实时间标定的分化速度。

confidence 不能替代 phase portrait，也不是完整的统计置信区间。

计算 latent time（可选）

Dynamical model 可以将不同基因的 latent time 整合成 gene-shared latent time：

scv.tl.latent_time(adata)

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

latent_time 是模型推断的相对内部时间。

RNA velocity 结果应该怎么验证？

RNA velocity 不能只看一张 UMAP stream。至少需要检查：

• 研究对象是否确实存在活跃、连续的状态转换。
• spliced 和 unspliced 是否来自未处理的 whole counts。
• 关键基因的 phase portrait 是否符合模型预期。
• 多套明显 kinetics 是否需要按 lineage 分开拟合。
• 主要方向是否符合 marker、采样时间点和已有生物学知识。
• UMAP 投影是否在边界处出现不合理回流。
• steady-state model 和 EM model 的主要结论是否严重冲突。

如果只有 UMAP 流线看起来合理，而 phase portrait 和独立证据不支持，应把结果视为探索性提示，不能作为确定的细胞命运结论。

对于需要定量比较终末状态或 fate probability 的项目，可以进一步考虑 CellRank。但这属于 RNA velocity 的下游分析，不在本节展开。

保存

保存结果：

out = (
    "/home/data/t090639/project/example/data/processed/"
    "10_example_rna_velocity.h5ad"
)

adata.write_h5ad(
    out,
    compression="gzip"
)

os.path.exists(out)

常见问题

只有普通 h5ad，能不能直接做 RNA velocity？

只有在对象中已经包含原始 spliced 和  unspliced layers 时才可以。普通 count matrix、Log(CP10k+1) 矩阵或 HVG 表达矩阵不能自动恢复剪接状态。

能不能直接对 PBMC 做 RNA velocity？

如果主要是稳定成熟细胞类型，并不存在被充分采样的连续转换，通常不适合。不能因为能够运行代码，就认为结果具有生物学意义。

EM model 是否一定比 steady-state model 更好？

不一定。EM model 放宽了必须观察到 steady state 等限制，但仍然依赖动力学模型和 phase portrait。模型假设不成立时，复杂模型也可能得到错误方向。

UMAP 上的箭头是否代表真实细胞移动路径？

不是。它们是高维 velocity vectors 的二维投影，会受到基因集合、UMAP 结构和绘图参数影响。

latent time 能不能解释成真实时间？

不能。它是模型推断的相对内部时间。