scvelo

scVelo — RNA ベロシティ解析

概要

scVelo は、シングルセル RNA-seq データの RNA ベロシティ解析における主要な Python パッケージです。mRNA スプライシングの速度論をモデル化し、未スプライス（pre-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 と スプライス済み RNA の両方のカウント行列が必要です。これらは以下を使用して生成します。

1. lamanno モードの STARsolo または kallisto|bustools
2. velocyto CLI：velocyto run10x / velocyto run
3. スプライス済み／未スプライスの出力を備えた alevin-fry / simpleaf

データは、layers["spliced"] と layers["unspliced"] を持つ AnnData オブジェクトに保存します。

標準的な 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. 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. ベロシティ推定 — 確率的モデル

確率的モデルは高速で、探索的解析に適しています。

# 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 (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 ベロシティ
adata.layers	fit_t	遺伝子ごと・細胞ごとに適合させた潜在時間
adata.obsm	velocity_umap	UMAP 上の 2D ベロシティベクトル
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 velocity の結果にばらつきが生じます
• Velocity に一貫性がある必要があります：矢印は既知の生物学的知見に沿うべきです。ランダムに見える場合は問題があることを示します
• k-NN のバンド幅が重要です：近傍数が少なすぎると velocity にノイズが生じ、多すぎると過度に平滑化されます
• 妥当性を確認します：ルート細胞（前駆細胞）では、マーカー遺伝子のスプライシング前／後の比率が高くなるはずです
• 力学的モデルには、異なる速度論的状態が必要です：明確な分化プロセスで最も効果を発揮します

トラブルシューティング

問題	解決策
スプライシング前のレイヤーがない	velocyto を再実行するか、--soloFeatures Gene Velocyto を指定して STARsolo を使用します
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）解析。