scVelo — RNA 速率分析
概述
scVelo 是用于单细胞 RNA-seq 数据 RNA 速率分析的主要 Python 包。它通过建模 mRNA 剪接动力学来推断细胞状态转换——利用未剪接（前体 mRNA）与已剪接（成熟 mRNA）丰度的比例判断每个细胞中基因是上调还是下调。这使得无需时间序列数据即可重建发育轨迹并识别细胞命运决策。
安装： uv pip install scvelo
关键资源：
- 文档：https://scvelo.readthedocs.io/
- GitHub：https://github.com/theislab/scvelo
- 论文：Bergen et al. (2020) Nature Biotechnology. PMID: 32747759
何时使用本技能
在以下情况下使用 scVelo：
- 从快照数据推断轨迹：确定细胞分化的方向
- 细胞命运预测：识别祖细胞及其下游命运
- 驱动基因识别：找到动力学最能解释观察到的轨迹的基因
- 发育生物学：模拟造血、神经发生、上皮-间质转化
- 潜在时间估计：根据剪接动力学对细胞进行伪时间排序
- 补充 Scanpy：为 UMAP 嵌入添加方向性信息
前提条件
scVelo 需要同时包含未剪接和已剪接 RNA 的计数矩阵。这些由以下工具生成：
- STARsolo 或 kallisto|bustools（使用
lamanno模式） - velocyto CLI：
velocyto run10x/velocyto run - alevin-fry / simpleaf（输出剪接/未剪接结果）
数据存储在 AnnData 对象中，包含 layers["spliced"] 和 layers["unspliced"]。
标准 RNA 速率工作流程
1. 设置和数据加载
import scvelo as scv
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
# 配置设置
scv.settings.verbosity = 3 # 显示计算步骤
scv.settings.presenter_view = True
scv.settings.set_figure_params('scvelo')
# 加载数据（包含剪接/未剪接层的 AnnData）
# 选项 A：从 loom 文件加载（velocyto 输出）
adata = scv.read("cellranger_output.loom", cache=True)
# 选项 B：将 velocyto loom 与已处理的 Scanpy AnnData 合并
adata_processed = sc.read_h5ad("processed.h5ad") # 包含 UMAP、聚类信息
adata_velocity = scv.read("velocyto.loom")
adata = scv.utils.merge(adata_processed, adata_velocity)
# 验证层
print(adata)
# obs × var：N × G
# layers: 'spliced', 'unspliced'（必需）
# obsm['X_umap']（可视化必需）
2. 预处理
# 过滤和归一化。从 scVelo 0.3 起，filter_and_normalize() 仅过滤基因并归一化每个细胞——
# 它不再接受 n_top_genes 参数，也不再进行 log 变换，
# 因此 log 步骤和 HVG 选择需要来自 Scanpy。
scv.pp.filter_and_normalize(
adata,
min_shared_counts=20 # 剪接+未剪接中的最小计数
)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)
# 计算一阶和二阶矩（均值和方差）
# 必须先计算 knn_connectivities
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30)
scv.pp.moments(
adata,
n_pcs=30,
n_neighbors=30
)
3. 速率估计——随机模型
随机模型速度快，适合探索性分析：
# 随机速率（较快，精度略低）
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
# 可视化
scv.pl.velocity_embedding_stream(
adata,
basis='umap',
color='leiden',
title="RNA 速率（随机模型）"
)
4. 速率估计——动力学模型（推荐）
动力学模型拟合完整的剪接动力学，精度更高：
# 恢复动力学（计算密集型；约 10-30 分钟处理 10K 细胞）
scv.tl.recover_dynamics(adata, n_jobs=4)
# 从动力学模型计算速率
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
5. 潜在时间
动力学模型可以计算共享的潜在时间（伪时间）：
# 计算潜在时间
scv.tl.latent_time(adata)
# 在 UMAP 上可视化潜在时间
scv.pl.scatter(
adata,
color='latent_time',
color_map='gnuplot',
size=80,
title='潜在时间'
)
# 识别按潜在时间排序的 top 基因
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
scv.pl.heatmap(
adata,
var_names=top_genes,
sortby='latent_time',
col_color='leiden',
n_convolve=100
)
6. 驱动基因分析
# 识别速率拟合最好的基因
scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
print(df.head(10))
# 速率和一致性
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
adata,
c=['velocity_length', 'velocity_confidence'],
cmap='coolwarm',
perc=[5, 95]
)
# 特定基因的相图
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2'],
ncols=3, figsize=(16, 4))
7. 速率箭头和伪时间
# UMAP 上的箭头图
scv.pl.velocity_embedding(
adata,
arrow_length=3,
arrow_size=2,
color='leiden',
basis='umap'
)
# 流线图（更清晰的可视化）
scv.pl.velocity_embedding_stream(
adata,
basis='umap',
color='leiden',
smooth=0.8,
min_mass=4
)
# 速率伪时间（替代潜在时间）
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
8. PAGA 轨迹图
# 基于速率信息的 PAGA 图
scv.tl.paga(adata, groups='leiden')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
# 绘制带有速率的 PAGA
scv.pl.paga(
adata,
basis='umap',
size=50,
alpha=0.1,
min_edge_width=2,
node_size_scale=1.5
)
完整工作流脚本
import scvelo as scv
import scanpy as sc
def run_rna_velocity(adata, n_top_genes=2000, mode='dynamical', n_jobs=4):
"""
完整的 RNA 速率工作流。
Args:
adata: 包含 'spliced' 和 'unspliced' 层以及 UMAP 在 obsm 中的 AnnData
n_top_genes: 用于速率的 top HVG 数量
mode: 'stochastic'（快速）或 'dynamical'（精确）
n_jobs: 动力学模型的并行任务数
Returns:
添加了速率信息的已处理 AnnData
"""
scv.settings.verbosity = 2
# 1. 预处理（scVelo 0.3 已将 log/HVG 从 filter_and_normalize 中移除）
scv.pp.filter_and_normalize(adata, min_shared_counts=20)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes, subset=True)
if 'neighbors' not in adata.uns:
sc.pp.neighbors(adata, n_neighbors=30)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
# 2. 速率估计
if mode == 'dynamical':
scv.tl.recover_dynamics(adata, n_jobs=n_jobs)
scv.tl.velocity(adata, mode=mode)
scv.tl.velocity_graph(adata)
# 3. 下游分析
if mode == 'dynamical':
scv.tl.latent_time(adata)
scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)
scv.tl.velocity_confidence(adata)
scv.tl.velocity_pseudotime(adata)
return adata
AnnData 中的关键输出字段
运行工作流后，会添加以下字段：
| 位置 | 键 | 描述 |
|---|---|---|
adata.layers | velocity | 每个基因每个细胞的 RNA 速率 |
adata.layers | fit_t | 每个基因每个细胞的拟合潜在时间 |
adata.obsm | velocity_umap | UMAP 上的二维速率向量 |
adata.obs | velocity_pseudotime | 基于速率的伪时间 |
adata.obs | latent_time | 动力学模型的潜在时间 |
adata.obs | velocity_length | 每个细胞的速率 |
adata.obs | velocity_confidence | 每个细胞的置信度分数 |
adata.var | fit_likelihood | 基因级模型拟合质量 |
adata.var | fit_alpha | 转录速率 |
adata.var | fit_beta | 剪接速率 |
adata.var | fit_gamma | 降解速率 |
adata.uns | velocity_graph | 细胞间转移概率矩阵 |
速率模型比较
| 模型 | 速度 | 精度 | 使用时机 |
|---|---|---|---|
stochastic | 快 | 中等 | 探索性分析；大型数据集 |
deterministic | 中 | 中等 | 简单的线性动力学 |
dynamical | 慢 | 高 | 出版级质量；识别驱动基因 |
最佳实践
- 从随机模式开始进行探索；最终分析时切换到动力学模式
- 需要良好的未剪接读数覆盖：短读长（< 100 bp）可能导致内含子覆盖不足
- 最少 2,000 个细胞：细胞数更少时 RNA 速率噪音较大
- 速率应具有一致性：箭头应遵循已知生物学规律；随机性表明存在问题
- k-NN 带宽很重要：邻居太少 → 速率噪音大；太多 → 过度平滑
- 合理性检查：根细胞（祖细胞）的标记基因应具有高的未剪接/已剪接比率
- 动力学模型需要清晰的动力学状态：最适合清晰的分化过程
故障排除
| 问题 | 解决方案 |
|---|---|
| 缺少未剪接层 | 重新运行 velocyto 或使用 STARsolo 的 --soloFeatures Gene Velocyto |
| 极少的速率基因 | 降低 min_shared_counts；检查测序深度 |
| 箭头看起来随机 | 尝试不同的 n_neighbors 或速率模型 |
| 动力学模型内存错误 | 设置 n_jobs=1；减少 n_top_genes |
| 处处为负速率 | 检查剪接/未剪接层是否互换 |
附加资源
- scVelo 文档：https://scvelo.readthedocs.io/
- 教程笔记本：https://scvelo.readthedocs.io/tutorials/
- GitHub：https://github.com/theislab/scvelo
- 论文：Bergen V et al. (2020) Nature Biotechnology. PMID: 32747759
- velocyto（预处理）：http://velocyto.org/
- CellRank（命运预测，扩展了 scVelo）：https://cellrank.readthedocs.io/
- dynamo（代谢标记替代方案）：https://dynamo-release.readthedocs.io/