trajectory-analysis

轨迹分析

概述

从单细胞 RNA-seq 数据重建连续的细胞状态转变。包含四个层面：

1. 扩散伪时间（scanpy.tl.dpt）——从用户指定的根节点出发，沿扩散流形对细胞排序。
2. PAGA ——粗粒度的聚类连接关系；用于初始化力导向布局。
3. RNA velocity（scVelo）——根据未剪接 / 已剪接 mRNA 的比例推断方向性流动。注意：scVelo 处于维护模式（不再添加新功能）；API 保持稳定。
4. CellRank v2 ——统一的核 API，将 velocity、伪时间、实时或 CytoTRACE 信号组合成马尔可夫链，用于分析宏观状态、终末状态和命运概率。CytoTRACEKernel 无需 velocity 即可推断轨迹。

输入要求为预处理后的 AnnData（已归一化、对数转换、筛选高变基因并计算邻居）。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. 前置条件——已剪接 / 未剪接计数

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 数据的新项目，可考虑使用 CellRank 的 CytoTRACEKernel 作为无需 velocity 的替代方案。

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_*']

常见错误

• 错误： 调用 sc.tl.dpt 前未设置 adata.uns['iroot'] 正确： 计算 DPT 前设置根细胞索引：adata.uns['iroot'] = np.flatnonzero(adata.obs['celltype'] == root_cluster)[0] 原因： 未设置根时，伪时间没有方向，结果毫无意义

• 错误： 在 AnnData 中没有 spliced/unspliced 层的情况下运行 scVelo velocity 正确： 执行 scv.utils.merge 后检查 adata.layers.keys()；确保这两个层都存在，且条形码匹配 原因： 计算 velocity 需要这两个层；缺少层会导致错误或静默失败

• 错误： 未先运行 scv.tl.recover_dynamics，就在 scv.tl.velocity 中使用 mode='dynamical' 正确： 始终先运行 scv.tl.recover_dynamics(adata)，再运行 scv.tl.velocity(adata, mode='dynamical') 原因： 动力学模型需要每个基因的速率参数，这些参数只能通过 recover_dynamics 计算

• 错误： 以 stochastic 模式运行 velocity 后再计算 scv.tl.latent_time 正确： 如果需要 latent time，请使用 mode='dynamical' 计算 velocity 原因： latent time 仅适用于动力学模型；stochastic 模式不会生成所需参数

• 错误： 未先运行 scv.pp.moments 就调用 scv.tl.velocity 正确： 计算 velocity 前运行 scv.pp.moments(adata, n_pcs=30, n_neighbors=30) 原因： Moments（Ms、Mu 层）是估算 velocity 所需的输入

• 错误： 计算聚类标签前运行 sc.tl.paga 正确： 先运行 sc.tl.leiden（或等效聚类方法）；通过 groups= 传入聚类键 原因： PAGA 计算的是现有群组之间的连通性；没有聚类就没有可连接的对象

• 错误： 认为 sc.pl.paga() 中的 threshold 参数会影响底层图 正确： 请注意，threshold 仅用于可视化——它控制绘制哪些边，而不影响计算结果 原因： 误解这一点会导致对聚类连通性做出错误判断

• 错误： 计算 velocity graph 或 CellRank 转移矩阵后再次运行 sc.pp.neighbors 正确： 如果必须更改邻居图，请随后重新计算 velocity_graph 和 CellRank kernels 原因： Velocity graph 和转移矩阵依赖邻居图；更改邻居图会使它们失效，但不会发出提示

• 错误： 在 CellRank v2 的 compute_macrostates() 中使用 n_cells 参数 正确： 使用 n_states——v2 API 已重命名此参数 原因： n_cells 是 v1 语法，在 v2 中会引发错误或被静默忽略

• 错误： 未检查 macrostates 就调用 predict_terminal_states() 正确： 始终运行 g.plot_macrostates()，并在预测终末状态前检查粗粒化转移 原因： 终末状态预测取决于宏状态的质量；糟糕的宏状态会生成糟糕的命运图

• 错误： Velocity 箭头看起来随机或不连贯时，仍强行进行基于 velocity 的分析 正确： 如果没有 velocity 信号，请切换到 CytoTRACEKernel 或 PseudotimeKernel 原因： 随机的 velocity 箭头说明数据缺少动力学信号；强行使用 velocity 会生成毫无意义的轨迹

参考资料

• 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