trajectory-analysis

軌跡解析

概要

単一細胞 RNA-seq から、連続的な細胞状態の遷移を再構築します。以下の4つのレイヤーがあります。

1. 拡散疑似時間（scanpy.tl.dpt） — ユーザーが指定した起点から、拡散多様体に沿って細胞を順序付けます。
2. PAGA — クラスター間の粗視化された接続性を示し、力学ベースのレイアウトを初期化します。
3. RNA velocity（scVelo） — 未スプライス／スプライス済み mRNA 比から方向性のある流れを推定します。注：scVelo はメンテナンスモードです（新機能の追加はありません）。API は安定しています。
4. CellRank v2 — velocity、疑似時間、実時間、または CytoTRACE のシグナルを組み合わせ、マルコフ連鎖を用いてマクロ状態、終末状態、運命確率を推定する統合カーネル API です。CytoTRACEKernel を使うと、velocity を使わずに軌跡を推定できます。

入力には、前処理済みの AnnData（正規化済み、対数変換済み、高変動遺伝子（HVG）抽出済み、近傍計算済み）を想定しています。Velocity にはさらに、velocyto または STARsolo から得た spliced／unspliced レイヤーが必要です。

使用方法

自然言語で解析内容を指定します。

• 「HSC クラスターを起点として拡散疑似時間を計算してください。」
• 「leiden クラスターに PAGA を実行し、力学ベースのレイアウトで再埋め込みしてください。」
• 「動力学 scVelo を実行し、velocity のストリームをプロットしてください。」
• 「CellRank と CytoTRACEKernel を使って、velocity データなしで終末状態を見つけてください。」
• 「VelocityKernel と ConnectivityKernel を組み合わせて、運命確率を計算してください。」

このスキルは、実行可能な Python スクリプトまたはノートブックのセルを生成します。

応答形式

1. 手法の選択 — ユーザーのデータに適した手法を示します（以下の判断フローを使用）
2. 前提条件の確認 — 必要なレイヤー／アノテーションが存在することを確認します
3. 実行可能なコード — コピーしてすぐに実行できる完全なスクリプト
4. 主要なパラメーター — 分かりにくい選択項目について説明します
5. 注意点 — 選択した手法に特有のよくある失敗を警告します

判断フロー — 軌跡解析手法の選択

Have spliced/unspliced layers?
├─ YES → RNA velocity (scVelo dynamical) → CellRank VelocityKernel
│         Best for: steady-state or dynamical kinetics inference
│
└─ NO
    ├─ Have experimental time points?
    │   └─ YES → CellRank RealTimeKernel
    │             Best for: time-course experiments (reprogramming, perturbation)
    │
    └─ NO time points
        ├─ Want directionality (root → terminal)?
        │   ├─ CytoTRACEKernel (uses gene-count gradient as differentiation proxy)
        │   │   Best for: differentiation without velocity; no extra input needed
        │   └─ PseudotimeKernel (uses precomputed DPT as direction)
        │       Best for: when you trust your DPT root choice
        │
        └─ Just want cluster connectivity?
            └─ PAGA (no directionality, shows which clusters connect)

基本概念

1. 前提条件 — spliced／unspliced カウント

Velocity には spliced と unspliced のレイヤーが必要です。velocyto または STARsolo を使って生成します。

# Option A: velocyto run10x
velocyto run10x \
    --samtools-threads 8 \
    /path/to/cellranger_outs/sample_id \
    /path/to/genes.gtf

# Option B: STARsolo (Velocyto-compatible matrices)
STAR --runThreadN 16 \
     --genomeDir /path/to/star_index \
     --readFilesIn R2.fq.gz R1.fq.gz \
     --readFilesCommand zcat \
     --soloType CB_UMI_Simple \
     --soloCBallowlist 3M-february-2018.txt \
     --soloFeatures Gene Velocyto \
     --soloCBstart 1 --soloCBlen 16 --soloUMIstart 17 --soloUMIlen 12 \
     --outSAMtype BAM SortedByCoordinate

loom を既存の AnnData とマージします。

import scvelo as scv

adata = sc.read_h5ad("processed.h5ad")
ldata = scv.read("sample_id.loom", cache=True)
# scv.utils.merge inner-joins on barcodes — ensure obs_names match between objects
adata = scv.utils.merge(adata, ldata)
assert "spliced" in adata.layers and "unspliced" in adata.layers

2. 拡散疑似時間（DPT）

import numpy as np
import scanpy as sc

sc.tl.diffmap(adata)

# Root cell: pick one cell in the known starting population
root_cluster = "HSC"
adata.uns["iroot"] = np.flatnonzero(adata.obs["celltype"] == root_cluster)[0]

sc.tl.dpt(adata)  # writes adata.obs['dpt_pseudotime']
sc.pl.umap(adata, color=["celltype", "dpt_pseudotime"])

3. PAGA

sc.tl.paga(adata, groups="leiden")

# threshold controls edge pruning: higher = fewer edges shown
# Default 0.01 shows almost everything; 0.03-0.1 is typical for clean plots
sc.pl.paga(adata, threshold=0.03, show=False)

# PAGA-initialized layout for cleaner trajectory visualization
sc.tl.draw_graph(adata, init_pos="paga")
sc.pl.draw_graph(adata, color=["leiden", "dpt_pseudotime"], legend_loc="on data")

# Gene trends along a chosen path
sc.pl.paga_path(
    adata,
    nodes=["HSC", "MPP", "GMP", "Mono"],
    keys=["Elane", "Mpo", "Gata1"],
)

PAGA のしきい値に関する注意： sc.pl.paga(threshold=...) は可視化（どのエッジを描画するか）にのみ影響し、基になるグラフには影響しません。接続性の重みは adata.uns['paga']['connectivities'] に格納されています。すべてのエッジを表示するには threshold=0 に設定し、表示を整理するには値を上げます。

4. RNA velocity（scVelo）

ステータス： scVelo はメンテナンスモードです。API は安定していますが、新機能の追加は予定されていません。Velocity データのない新規プロジェクトでは、velocity を使わない代替手法として CellRank の CytoTRACEKernel を検討してください。

import scvelo as scv

scv.settings.n_jobs = 8

# Preprocessing — operates on spliced/unspliced layers
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# Dynamical mode (recommended): fits full kinetic model per gene
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)

scv.pl.velocity_embedding_stream(adata, basis="umap", color="celltype")

# Latent time (dynamical mode only)
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", color_map="gnuplot")

5. CellRank v2 — 統合カーネル API

CellRank v2 は、統合カーネルインターフェイスを提供します。すべてのカーネルは次の特性を備えています。

• AnnData を受け取り、.compute_transition_matrix() を呼び出します。
• + と *（加重和）で組み合わせられます。
• マクロ状態の分解には cr.estimators.GPCCA を使用します。

5a. VelocityKernel（scVelo の velocity が必要）

import cellrank as cr
from cellrank.kernels import VelocityKernel, ConnectivityKernel

vk = VelocityKernel(adata).compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * vk + 0.2 * ck  # weighted combination

g = cr.estimators.GPCCA(kernel)
g.compute_schur(n_components=20)
g.compute_macrostates(n_states=5, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()
g.plot_fate_probabilities(same_plot=False)

# Lineage driver genes
drivers = g.compute_lineage_drivers(lineages=["Mono"], use_raw=False)

5b. CytoTRACEKernel（velocity 不要、v2 で追加）

遺伝子数の勾配を、分化ポテンシャルの代替指標として使用します。Velocity データは不要です。

from cellrank.kernels import CytoTRACEKernel, ConnectivityKernel

ctk = CytoTRACEKernel(adata).compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * ctk + 0.2 * ck

g = cr.estimators.GPCCA(kernel)
g.compute_schur(n_components=15)
g.compute_macrostates(n_states=4, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()

5c. RealTimeKernel（実験の時点）

細胞を離散的な時点ごとに採取する時系列実験向けです。

from cellrank.kernels import RealTimeKernel

# adata.obs["day"] must contain numeric time labels (e.g., 0, 2, 4, 7)
rtk = RealTimeKernel(adata, time_key="day").compute_transition_matrix()

g = cr.estimators.GPCCA(rtk)
g.compute_macrostates(n_states=6, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()

5d. PseudotimeKernel（DPT による方向付け）

事前に計算した疑似時間を使って方向性を定義します。

from cellrank.kernels import PseudotimeKernel, ConnectivityKernel

# Requires adata.obs["dpt_pseudotime"] from sc.tl.dpt
pk = PseudotimeKernel(adata, time_key="dpt_pseudotime").compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * pk + 0.2 * ck

g = cr.estimators.GPCCA(kernel)
g.compute_macrostates(n_states=4, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()

5e. カーネルの組み合わせルール

# All kernels support arithmetic combination after .compute_transition_matrix()
combined = 0.6 * vk + 0.2 * ck + 0.2 * pk  # weights must sum to 1

# Inspect kernel quality
combined.plot_projection(basis="umap")  # visualize transition arrows

6. 空間近傍（squidpy）

空間トランスクリプトミクスでは、kNN を空間グラフに置き換えます。

import squidpy as sq

sq.gr.spatial_neighbors(adata, coord_type="generic", n_neighs=6)
# Feed spatial connectivities into CellRank ConnectivityKernel

Quick Reference

Task	Call	Output
Diffusion map	sc.tl.diffmap(adata)	.obsm['X_diffmap']
Set root	adata.uns['iroot'] = idx	—
Pseudotime	sc.tl.dpt(adata)	.obs['dpt_pseudotime']
PAGA	sc.tl.paga(adata, groups='leiden')	.uns['paga']
PAGA plot	sc.pl.paga(adata, threshold=0.03)	figure (threshold=visual only)
PAGA-init layout	sc.tl.draw_graph(adata, init_pos='paga')	.obsm['X_draw_graph_fa']
Merge loom	scv.utils.merge(adata, ldata)	combined AnnData (inner-join on barcodes)
Velocity preproc	scv.pp.filter_and_normalize; scv.pp.moments	moment layers
Dynamics fit	scv.tl.recover_dynamics(adata)	fit params
Velocity	scv.tl.velocity(adata, mode='dynamical')	.layers['velocity']
Velocity graph	scv.tl.velocity_graph(adata)	.uns['velocity_graph']
Stream plot	scv.pl.velocity_embedding_stream(adata, basis='umap')	figure
Latent time	scv.tl.latent_time(adata)	.obs['latent_time']
VelocityKernel	VelocityKernel(adata).compute_transition_matrix()	kernel
CytoTRACEKernel	CytoTRACEKernel(adata).compute_transition_matrix()	kernel (no velocity needed)
RealTimeKernel	RealTimeKernel(adata, time_key="day").compute_transition_matrix()	kernel
PseudotimeKernel	PseudotimeKernel(adata, time_key="dpt_pseudotime").compute_transition_matrix()	kernel
Combine kernels	0.8 * vk + 0.2 * ck	combined kernel
Estimator	cr.estimators.GPCCA(kernel)	estimator
Macrostates	g.compute_macrostates(n_states=5, cluster_key="celltype")	.macrostates
Terminal states	g.predict_terminal_states()	.terminal_states
Fate probs	g.compute_fate_probabilities()	.fate_probabilities
Lineage drivers	g.compute_lineage_drivers(lineages=["Mono"], use_raw=False)	DataFrame
Spatial graph	sq.gr.spatial_neighbors(adata)	.obsp['spatial_*']

よくある間違い

• 誤り： adata.uns['iroot'] を設定せずに sc.tl.dpt を呼び出す 正しい方法： DPT を計算する前にルート細胞のインデックスを設定します：adata.uns['iroot'] = np.flatnonzero(adata.obs['celltype'] == root_cluster)[0] 理由： ルートがないと擬似時間に方向性がなく、結果に意味がありません

• 誤り： AnnData に spliced／unspliced レイヤーがない状態で scVelo の速度解析を実行する 正しい方法： scv.utils.merge の実行後に adata.layers.keys() を確認し、バーコードが一致する両方のレイヤーが存在することを確かめます 理由： 速度計算には両方のレイヤーが必要です。レイヤーがないと、エラーや通知のない失敗が発生します

• 誤り： scv.tl.recover_dynamics を先に実行せずに、scv.tl.velocity で mode='dynamical' を使用する 正しい方法： 必ず scv.tl.velocity(adata, mode='dynamical') の前に scv.tl.recover_dynamics(adata) を実行します 理由： ダイナミカルモデルには遺伝子ごとの速度パラメーターが必要で、これらは recover_dynamics によってのみ計算されます

• 誤り： stochastic モードで速度解析を実行した後に scv.tl.latent_time を計算する 正しい方法： 潜在時間が必要な場合は、速度解析で mode='dynamical' を使用します 理由： 潜在時間はダイナミカルモデルでのみ定義されます。stochastic モードでは必要なパラメーターが生成されません

• 誤り： 先に sc.pp.moments を実行せずに scv.tl.velocity を呼び出す 正しい方法： 速度計算の前に sc.pp.moments(adata, n_pcs=30, n_neighbors=30) を実行します 理由： モーメント（Ms、Mu レイヤー）は速度推定に必要な入力です

• 誤り： クラスターラベルを計算する前に sc.tl.paga を実行する 正しい方法： 先に sc.tl.leiden（または同等のクラスタリング）を実行し、groups= でクラスターのキーを渡します 理由： PAGA は既存のグループ間の接続性を計算します。クラスターがなければ、接続する対象がありません

• 誤り： sc.pl.paga() の threshold パラメーターが基礎となるグラフに影響すると考える 正しい方法： threshold は表示にのみ関係し、描画するエッジを制御するものであって、計算するエッジを制御するものではないと理解します 理由： これを誤解すると、クラスターの接続性について誤った主張につながります

• 誤り： 速度グラフまたは CellRank の遷移行列を計算した後に sc.pp.neighbors を再実行する 正しい方法： 近傍を変更する必要がある場合は、その後に velocity_graph と CellRank カーネルを再計算します 理由： 速度グラフと遷移行列は近傍グラフに依存します。近傍グラフを変更すると、それらは通知なく無効になります

• 誤り： CellRank v2 の compute_macrostates() で n_cells パラメーターを使用する 正しい方法： v2 API ではこのパラメーター名が変更されているため、n_states を使用します 理由： n_cells は v1 の構文であり、v2 ではエラーになるか、通知なく無視されます

• 誤り： マクロステートを先に確認せずに predict_terminal_states() を呼び出す 正しい方法： 終末状態を予測する前に、必ず g.plot_macrostates() を実行し、粗視化された遷移を確認します 理由： 終末状態の予測はマクロステートの品質に依存します。質の悪いマクロステートからは、質の悪い運命マップが生成されます

• 誤り： 速度の矢印がランダムまたは一貫性を欠いて見える場合に、速度解析を無理に使用する 正しい方法： 速度シグナルがない場合は、CytoTRACEKernel または PseudotimeKernel に切り替えます 理由： 速度の矢印がランダムな場合、データに動態シグナルが欠けていることを示しています。速度解析を無理に適用すると、意味のない軌跡が生成されます

参考文献

• CellRank 2：Weiler et al. Nat Methods 2024, https://doi.org/10.1038/s41592-023-02088-9
• scVelo：Bergen et al. Nat Biotechnol 2020, https://doi.org/10.1038/s41587-020-0591-3
• CellRank：Lange et al. Nat Methods 2022, https://doi.org/10.1038/s41592-021-01346-6
• PAGA：Wolf et al. Genome Biol 2019, https://doi.org/10.1186/s13059-019-1663-x
• velocyto：La Manno et al. Nature 2018, https://doi.org/10.1038/s41586-018-0414-6