scvelo

scVelo — RNA Velocity 分析

概述

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

安装： uv pip install scvelo

关键资源：

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

何时使用此技能

以下情况可使用 scVelo：

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

前提条件

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

1. 使用 lamanno 模式的 STARsolo 或 kallisto|bustools
2. velocyto 命令行工具：velocyto run10x / velocyto run
3. 输出已剪接/未剪接数据的 alevin-fry / simpleaf

数据存储在包含 layers["spliced"] 和 layers["unspliced"] 的 AnnData 对象中。

标准 RNA Velocity 工作流

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. As of scVelo 0.3, filter_and_normalize() only filters
# genes and normalizes per cell -- it no longer takes n_top_genes and no longer
# log-transforms, so the log step and HVG selection come from Scanpy.
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20    # Minimum counts in spliced+unspliced
)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)

# Compute first and second order moments (means and variances)
# knn_connectivities must be computed first
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30)
scv.pp.moments(
    adata,
    n_pcs=30,
    n_neighbors=30
)

3. Velocity 估算——随机模型

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

# 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. Velocity 估算——动力学模型（推荐）

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

# 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. Velocity 箭头与伪时间

# 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 (scVelo 0.3 dropped log/HVG from 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. 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	慢	高	适用于发表级分析；可识别驱动基因

最佳实践

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

故障排除

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

更多资源

• 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/

相关技能

• scanpy：标准的单细胞 RNA-seq 分析流程。
• openpiv：使用 OpenPIV 进行粒子图像测速（PIV）分析。