AI优化细胞轨迹分析

细胞轨迹分析的核心问题

单细胞轨迹分析旨在重建细胞从一种状态（如干细胞）转变为另一种状态（如分化终末细胞）的连续过程。传统方法如Monocle、PAGA、RNA velocity面临：

难以处理复杂的分支结构（如多谱系分化）
对起始点和终点的选择敏感
无法建模循环轨迹（如细胞周期）

AI赋能轨迹推断的三种范式

1. 基于图神经网络（GNN）的拓扑学习

代表工具：CellPath, DeepVelo
核心思想：将细胞看作图中的节点，使用GNN学习细胞-细胞转换的概率邻接矩阵

# 使用scVelo + GNN增强RNA velocity
import scvelo as scv
import torch_geometric

adata = scv.read('data.h5ad')
scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)

# 标准velocity分析
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)

# 用GNN优化velocity图
# （需要额外的GNN模块，此处为示意）
from deepvelo import GNNRefiner
refined_graph = GNNRefiner.refine(adata.uns['velocity_graph'])

2. 基于VAE的连续轨迹嵌入

代表工具：Palantir, STREAM, scVI dynamics
核心思想：使用变分自编码器学习低维流形，在流形上拟合平滑轨迹

import palantir

# 预处理
adata = sc.read_h5ad('hematopoiesis.h5ad')
sc.pp.neighbors(adata, n_neighbors=30)
sc.tl.diffmap(adata)

# 运行Palantir
pr_res = palantir.core.run_palantir(
    adata,
    early_cell='HSC_001',  # 指定起始细胞
    num_waypoints=500
)

# 可视化分化势能和伪时序
palantir.plot.plot_palantir_results(pr_res, adata)

3. 基于物理模型的轨迹预测

代表工具：Waddington-OT, CellRank
核心思想：结合最优传输理论和RNA velocity，建模细胞命运转换的能量景观

import cellrank as cr

# 计算细胞转换核（transition kernel）
cr.tl.terminal_states(adata, cluster_key='celltype')
cr.tl.initial_states(adata, cluster_key='celltype')

# 计算吸收概率（细胞最终分化到各终末状态的概率）
cr.tl.lineages(adata)

# 可视化Waddington景观
cr.pl.lineages(adata, same_plot=True)

推荐工作流：CellRank 2.0

CellRank整合了RNA velocity、细胞-细胞相似性和实验时间信息，使用马尔可夫链建模细胞状态转换，支持多种计算后端（velocity、伪时序、Waddington-OT）。

完整示例：造血分化轨迹

import scanpy as sc
import scvelo as scv
import cellrank as cr

# 1. 加载数据并预处理
adata = sc.datasets.pbmc3k()
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)

# 2. 降维和聚类
sc.pp.highly_variable_genes(adata)
sc.tl.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata)

# 3. RNA velocity分析
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

# 4. CellRank推断细胞命运
cr.tl.terminal_states(adata, cluster_key='leiden', method='leiden')
cr.tl.initial_states(adata, cluster_key='leiden')
cr.tl.lineages(adata)

# 5. 识别关键转换驱动基因
cr.tl.lineage_drivers(adata)

# 6. 可视化
cr.pl.lineages(adata, same_plot=False)
scv.pl.velocity_embedding_stream(adata, basis='umap', color='leiden')

轨迹分析的质量控制

生物学验证：检查轨迹是否符合已知的分化顺序（如干细胞→祖细胞→成熟细胞）
Marker基因表达梯度：沿伪时序绘制已知marker基因，验证表达模式
velocity一致性：检查velocity方向是否与伪时序方向一致
分支点基因：识别命运决定点的差异表达基因，验证是否为已知调控因子

特殊场景处理

循环轨迹（如细胞周期）

# 使用Cyclum检测周期性轨迹
from cyclum import Cyclum

cy = Cyclum(adata.X)
theta = cy.fit_predict()  # 返回每个细胞在周期中的相位
adata.obs['cell_cycle_phase'] = theta

多条件比较（如治疗前后）

# 使用CellRank比较两个条件的轨迹差异
cr.tl.terminal_states(adata, cluster_key='celltype', method='eigengap')
cr.tl.lineages(adata, cluster_key='celltype')

# 分别计算每个条件的吸收概率
for condition in ['control', 'treated']:
    adata_subset = adata[adata.obs['condition'] == condition]
    cr.tl.lineages(adata_subset)