alterlab-scvelo

scVelo — RNA Velocity 解析

概要

scVelo は、シングルセル RNA-seq データの RNA Velocity 解析を行う主要な Python パッケージです。mRNA スプライシングの動態をモデル化し、未スプライシング（pre-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 で一度もテストされていない環境を問題なく構築してしまいます。

実用上は、現行の scverse スタックに組み込むのではなく、NumPy のバージョンを固定して独立した環境に インストールしてください。NumPy 2 で dtype やスカラー代入のエラーが発生した場合、特に mode="stochastic" で 発生しやすいです。このモードでは、2 次モーメントの回帰処理で、最小二乗法の結果として返される 1 要素の値を スカラーに代入します。numpy<2（例：numpy==1.26.4）に固定して、再実行してください。結果に問題がないと 決めつけず、速度推定結果がひそかに壊れていないか確認してください。stochastic モードと dynamical モードを 比較し、ストリームラインが既知の生物学的知見と一致していることを確かめてください。

主なリソース：

• ドキュメント：https://scvelo.readthedocs.io/
• GitHub：https://github.com/theislab/scvelo
• 論文：Bergen et al.（2020）Nature Biotechnology。PMID: 32747759

このスキルを使う場面

次のような場合は scVelo を使用してください。

• スナップショットデータからの軌跡推定：細胞が分化する方向を判定する
• 細胞運命の予測：前駆細胞と、その下流の細胞運命を特定する
• ドライバー遺伝子の特定：観測された軌跡を最もよく説明する動態を持つ遺伝子を見つける
• 発生生物学：造血、神経発生、上皮間葉転換をモデル化する
• 潜在時間の推定：スプライシングの動態から導出した擬似時間に沿って細胞を並べる
• Scanpy の補完：UMAP 埋め込みに方向情報を追加する

使用しない場面

シナリオ	代わりに使用するもの
QC、正規化、PCA/UMAP、Leiden クラスタリング、マーカー遺伝子（Velocity を使用しない場合）	alterlab-scanpy
.h5ad オブジェクト、レイヤー、結合の読み取り／データ加工	alterlab-anndata
確率的な統合／バッチ補正／モデルベースの DE（scVI、scANVI、veloVI）	alterlab-scvi-tools
FASTQ からスプライシング済み／未スプライシングのカウントレイヤーを作成すること	alterlab-rnaseq-quant
発現から遺伝子制御ネットワークを推定すること	alterlab-arboreto

前提条件

scVelo では、RNA の未スプライシングとスプライシング済みの両方のカウント行列が必要です。これらは次の方法で生成します。

1. lamanno モードの STARsolo または kallisto|bustools
2. velocyto CLI：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 (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. 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
    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 上の 2D 速度ベクトル
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 に切り替えます
• 未スプライスリードを十分にカバーする必要があります：短いリード（< 100 bp）ではイントロンのカバレッジが不足することがあります
• 最低 2,000 細胞が必要です：細胞数がこれより少ないと、RNA velocity はノイズが多くなります
• Velocity に一貫性があることを確認します：矢印は既知の生物学的な動態に沿う必要があります。ランダムな方向を示す場合は問題があることを示しています
• k-NN の帯域幅が重要です：近傍数が少なすぎると velocity にノイズが増え、多すぎると過度に平滑化されます
• 妥当性を確認します：ルート細胞（前駆細胞）では、マーカー遺伝子の未スプライス／スプライス比が高くなるはずです
• Dynamical モデルには明確に異なる動態状態が必要です：明確な分化過程で最も効果を発揮します

トラブルシューティング

問題	解決策
未スプライスレイヤーが見つからない	velocyto を再実行するか、--soloFeatures Gene Velocyto を指定して STARsolo を使用します
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 et al. (2020) Nature Biotechnology. PMID: 32747759
• velocyto（前処理）：http://velocyto.org/
• CellRank 2（運命予測。velocity カーネルを使用しますが、RNA velocity がなくても、擬似時間／実時間／代謝標識カーネルを介して動作します。運命確率を得るための、現在積極的にメンテナンスされている方法です）：https://cellrank.readthedocs.io/
• dynamo（代謝標識を用いる代替手段）：https://dynamo-release.readthedocs.io/