alterlab-scvelo

scVelo — RNA 速度分析

概述

scVelo 是用于单细胞 RNA-seq 数据 RNA 速度分析的领先 Python 包。它通过模拟 mRNA 剪接动力学来推断细胞状态转变——利用未剪接（前体 mRNA）与已剪接（成熟 mRNA）的丰度比，判断每个细胞中的基因是上调还是下调。借此，无需时间序列数据即可重建发育轨迹并识别细胞命运决定。

安装： uv pip install "scvelo==0.3.4" — 0.3.4（上传于 2026-02-24）是 当前版本，该软件包目前仅得到最低限度的维护：已发布的版本说明止于 0.2.5，上游 CI 仍以 Python 3.9/3.10 为目标，且其依赖声明（numpy>=1.17、scanpy>=1.5、anndata>=0.7.5）都没有版本上限，因此 pip 会轻易构建出一个 scVelo 从未经过测试的环境。

实际操作中，应在独立环境中安装，并固定 NumPy 版本，而不要将其装进当前的 scverse 技术栈。如果在 NumPy 2 下遇到 dtype 或标量赋值错误——最常见于 mode="stochastic"，其二阶矩回归会将一个单元素最小二乘结果赋给标量——请固定 numpy<2（例如 numpy==1.26.4）并重新运行。检查速度结果是否存在静默故障，不要想当然地认为结果正常：比较随机模型与动力学模型，并确认流线与已知生物学规律一致。

关键资源：

• 文档：https://scvelo.readthedocs.io/
• GitHub：https://github.com/theislab/scvelo
• 论文：Bergen 等（2020），Nature Biotechnology。PMID: 32747759

何时使用此技能

在以下情况下使用 scVelo：

• 根据快照数据推断轨迹：确定细胞分化的方向
• 预测细胞命运：识别祖细胞及其下游命运
• 识别驱动基因：找出最能解释观测轨迹动态的基因
• 发育生物学：模拟造血、神经发生、上皮间质转化
• 估计潜在时间：根据剪接动态推导的伪时间对细胞排序
• 作为 Scanpy 的补充：为 UMAP 嵌入添加方向信息

不适用的情况

场景	改用
QC、归一化、PCA/UMAP、Leiden 聚类、标记基因（不涉及速度分析）	alterlab-scanpy
读取/处理 .h5ad 对象、图层和合并数据	alterlab-anndata
概率整合 / 批次校正 / 基于模型的差异表达（scVI、scANVI、veloVI）	alterlab-scvi-tools
从 FASTQ 生成剪接/未剪接计数图层	alterlab-rnaseq-quant
根据表达量推断基因调控网络	alterlab-arboreto

前置条件

scVelo 需要未剪接和已剪接 RNA 的计数矩阵。这些矩阵可通过以下工具生成：

1. STARsolo 或 kallisto|bustools 的 lamanno 模式
2. velocyto CLI：velocyto run10x / velocyto run
3. 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

# Configure settings
scv.settings.verbosity = 3       # Show computation steps
scv.settings.presenter_view = True
scv.settings.set_figure_params('scvelo')

# Load data (AnnData with spliced/unspliced layers)
# Option A: Load from loom (velocyto output)
adata = scv.read("cellranger_output.loom", cache=True)

# Option B: Merge velocyto loom with Scanpy-processed AnnData
adata_processed = sc.read_h5ad("processed.h5ad")  # Has UMAP, clusters
adata_velocity = scv.read("velocyto.loom")
adata = scv.utils.merge(adata_processed, adata_velocity)

# Verify layers
print(adata)
# obs × var: N × G
# layers: 'spliced', 'unspliced' (required)
# obsm['X_umap'] (required for visualization)

2. 预处理

# Filter and normalize (follows Scanpy conventions)
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20,   # Minimum counts in spliced+unspliced
    n_top_genes=2000        # Top highly variable genes
)

# Compute first and second order moments (means and variances).
# scv.pp.moments runs PCA + a kNN graph internally if they're absent, so
# you do NOT need a separate sc.pp.neighbors call here. Calling moments with
# n_pcs/n_neighbors and ALSO running sc.pp.neighbors first just recomputes the
# graph with possibly mismatched params — let moments own it on a fresh object.
scv.pp.moments(
    adata,
    n_pcs=30,
    n_neighbors=30
)

3. 速度估计——随机模型

随机模型速度快，适合探索性分析：

# Stochastic velocity (faster, less accurate)
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)

# Visualize
scv.pl.velocity_embedding_stream(
    adata,
    basis='umap',
    color='leiden',
    title="RNA Velocity (Stochastic)"
)

4. 速度估计——动力学模型（推荐）

动力学模型拟合完整的剪接动力学，准确度更高：

# Recover dynamics (computationally intensive; ~10-30 min for 10K cells)
scv.tl.recover_dynamics(adata, n_jobs=4)

# Compute velocity from dynamical model
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

5. 潜在时间

动力学模型可用于计算共享的潜在时间（伪时间）：

# Compute latent time
scv.tl.latent_time(adata)

# Visualize latent time on UMAP
scv.pl.scatter(
    adata,
    color='latent_time',
    color_map='gnuplot',
    size=80,
    title='Latent time'
)

# Identify top genes ordered by latent time
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. 驱动基因分析

# Identify genes with highest velocity fit
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))

# Speed and coherence
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
    adata,
    c=['velocity_length', 'velocity_confidence'],
    cmap='coolwarm',
    perc=[5, 95]
)

# Phase portraits for specific genes
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2'],
               ncols=3, figsize=(16, 4))

7. 速度箭头与伪时间

# Arrow plot on UMAP
scv.pl.velocity_embedding(
    adata,
    arrow_length=3,
    arrow_size=2,
    color='leiden',
    basis='umap'
)

# Stream plot (cleaner visualization)
scv.pl.velocity_embedding_stream(
    adata,
    basis='umap',
    color='leiden',
    smooth=0.8,
    min_mass=4
)

# Velocity pseudotime (alternative to latent time)
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')

8. PAGA 轨迹图

# PAGA graph with velocity-informed transitions
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}')

# Plot PAGA with velocity
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):
    """
    Complete RNA velocity workflow.

    Args:
        adata: AnnData with 'spliced' and 'unspliced' layers, UMAP in obsm
        n_top_genes: Number of top HVGs for velocity
        mode: 'stochastic' (fast) or 'dynamical' (accurate)
        n_jobs: Parallel jobs for dynamical model

    Returns:
        Processed AnnData with velocity information
    """
    scv.settings.verbosity = 2

    # 1. Preprocessing
    scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=n_top_genes)

    if 'neighbors' not in adata.uns:
        sc.pp.neighbors(adata, n_neighbors=30)

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

    # 2. Velocity estimation
    if mode == 'dynamical':
        scv.tl.recover_dynamics(adata, n_jobs=n_jobs)

    scv.tl.velocity(adata, mode=mode)
    scv.tl.velocity_graph(adata)

    # 3. Downstream analyses
    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 velocity
adata.layers	fit_t	每个细胞中每个基因的拟合潜在时间
adata.obsm	velocity_umap	UMAP 上的二维 velocity 向量
adata.obs	velocity_pseudotime	根据 velocity 得出的伪时间
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	细胞间转移概率矩阵

Velocity 模型比较

模型	速度	准确度	适用场景
stochastic	快	中等	探索性分析；大型数据集
deterministic	中等	中等	简单的线性动力学
dynamical	慢	高	适用于发表级别的分析；可识别驱动基因

最佳实践

• 从 stochastic 模式开始进行探索；最终分析时切换到 dynamical
• 需要保证未剪接 reads 的覆盖度良好：较短的 reads（< 100 bp）可能无法覆盖内含子
• 至少需要 2,000 个细胞：细胞数量较少时，RNA velocity 的噪声较大
• Velocity 应具有连贯性：箭头应符合已知的生物学规律；随机分布则表明存在问题
• k-NN 带宽很重要：邻居过少 → velocity 噪声较大；邻居过多 → 过度平滑
• 合理性检查：根细胞（祖细胞）的标记基因应具有较高的未剪接／已剪接比例
• Dynamical 模型要求存在不同的动力学状态：对于分化过程明确的情况，效果最佳

故障排除

问题	解决方案
缺少未剪接层	重新运行 velocyto，或使用 STARsolo 并添加 --soloFeatures Gene Velocyto
Velocity 基因非常少	降低 min_shared_counts；检查测序深度
箭头看起来随机	尝试不同的 n_neighbors 或 velocity 模型
Dynamical 出现内存错误	设置 n_jobs=1；减少 n_top_genes
到处都是负 velocity	检查已剪接／未剪接层是否颠倒

更多资源

• scVelo 文档：https://scvelo.readthedocs.io/
• 教程笔记本：https://scvelo.readthedocs.io/en/stable/VelocityBasics.html
• GitHub：https://github.com/theislab/scvelo
• 论文：Bergen V 等（2020），Nature Biotechnology。PMID：32747759
• velocyto（预处理）：http://velocyto.org/
• CellRank 2（命运预测；使用 velocity kernel，但也可通过伪时间／实时／代谢标记 kernel，在不使用 RNA velocity 的情况下运行——这是目前维护活跃的命运概率分析方案）：https://cellrank.readthedocs.io/
• dynamo（代谢标记替代方案）：https://dynamo-release.readthedocs.io/