scvelo

scVelo — RNA 速度分析

概述

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

安装： 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 的计数矩阵。这些矩阵可通过以下工具生成：

1. 使用 lamanno 模式的 STARsolo 或 kallisto|bustools
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)
# 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. 估算速度 — 随机模型

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

# 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 速度
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	慢	高	适用于发表级成果；可识别驱动基因

最佳实践

• 探索时从随机模式开始；最终分析时切换到动力学模式
• 需要良好的未剪接 reads 覆盖度：短 reads（< 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 等（2020），Nature Biotechnology。PMID：32747759
• velocyto（预处理）：http://velocyto.org/
• CellRank（命运预测，可扩展 scVelo）：https://cellrank.readthedocs.io/
• dynamo（代谢标记替代方案）：https://dynamo-release.readthedocs.io/