05 轨迹推断与拟时序分析

聚类把细胞分成若干离散的 cluster，但发育、分化、应激响应这类现象更自然的描述是"细胞沿着一条连续路径逐渐变化"。轨迹推断（trajectory inference）的目标就是把单细胞数据里的连续结构显式地建模成图或曲线，然后给每个细胞一个在这条路径上的位置 —— 就是拟时序（pseudotime）。

这是一个"真正的时间轴"在单细胞里其实不可观测的近似：所有细胞都是同一刻被测的，我们只是用基因表达的变化推测出它们在分化进程里前后的相对位置。本章走三种主流方法的典型用法：Monocle3、Slingshot、以及基于 RNA velocity 的 scVelo。

这类分析能回答什么

• 分化路径上关键的基因开关在哪个拟时点打开
• 某个 cluster 在分化树上处在分支点还是端点
• 两个最终命运（比如 T vs B）的分叉什么时候发生
• 哪些基因是"驱动基因"、哪些只是被动跟随

如果你的项目里没有明显的连续过程（比如只是健康样本的静态图谱），拟时序分析往往得不出有用的结论。先想清楚要回答什么问题再选方法。

方法选择

实际选择建议：

• 新项目从 Monocle3 开始。它是现在最标准的选择。
• 需要 RNA velocity 的"方向性"（看细胞正在往哪变）时用 scVelo，前提是你有 BAM 文件或能跑 velocyto。
• 只是想给细胞一个"早-晚"排序、不关心轨迹形状的话 CytoTRACE 最简单。

Monocle3：完整流程

安装

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c(
  "BiocGenerics", "DelayedArray", "DelayedMatrixStats", "limma",
  "S4Vectors", "SingleCellExperiment", "SummarizedExperiment",
  "batchelor", "HDF5Array", "terra", "ggrastr"
))

install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle3")

从 Seurat 对象出发

Monocle3 的核心对象叫 cell_data_set（CDS）。如果前面用 Seurat 分析过，可以直接转：

library(Seurat)
library(monocle3)
library(SeuratWrappers)

cds <- as.cell_data_set(seurat_obj)

也可以手动从 counts 矩阵建：

cds <- new_cell_data_set(
  expression_data = seurat_obj@assays$RNA@counts,
  cell_metadata   = seurat_obj@meta.data,
  gene_metadata   = data.frame(
    gene_short_name = rownames(seurat_obj),
    row.names       = rownames(seurat_obj)
  )
)

预处理、降维和聚类

cds <- preprocess_cds(cds, num_dim = 50)
cds <- reduce_dimension(cds, reduction_method = "UMAP")
cds <- cluster_cells(cds, resolution = 1e-3)

plot_cells(cds, color_cells_by = "cell_type")
plot_cells(cds, color_cells_by = "cluster")

Monocle3 自己有一套 UMAP 和聚类。如果你已经在 Seurat 里做过一套，结果通常很相似，但 Monocle3 的 UMAP 会是后续 learn_graph 的输入，不要跳过这一步。

学习轨迹图

cds <- learn_graph(cds)

plot_cells(cds,
  color_cells_by       = "cell_type",
  label_groups_by_cluster = FALSE,
  label_leaves         = FALSE,
  label_branch_points  = FALSE
)

learn_graph 会在 UMAP 上拟合一条主干曲线（principal graph），分支点和叶子结点由算法自动判断。

选择起点并排序

拟时序需要一个起点（"分化程度最低的那个 cluster"）。两种方式：

# 方式 1：点图里手动点
cds <- order_cells(cds)

# 方式 2：根据已有注释自动选
cds <- order_cells(cds, root_cells = colnames(cds)[cds$cell_type == "Stem"])

plot_cells(cds,
  color_cells_by = "pseudotime",
  label_cell_groups = FALSE,
  label_leaves      = FALSE,
  label_branch_points = FALSE
)

找随拟时序变化的基因

track_genes <- graph_test(cds, neighbor_graph = "principal_graph")

track_genes_sig <- track_genes |>
  dplyr::filter(q_value < 0.05) |>
  dplyr::arrange(q_value)

head(track_genes_sig, 20)

plot_cells(
  cds,
  genes                    = head(track_genes_sig$gene_short_name, 4),
  show_trajectory_graph    = FALSE,
  label_cell_groups        = FALSE
)

基因模块

把这些随拟时序变化的基因聚成共表达模块，能看清分化不同阶段各自"开启"了哪一组程序：

gene_modules <- find_gene_modules(
  cds[track_genes_sig$gene_short_name, ],
  resolution = 1e-2
)

plot_cells(
  cds,
  genes                 = gene_modules,
  label_cell_groups     = FALSE,
  show_trajectory_graph = FALSE
)

Slingshot：轻量替代

分支结构简单、已经聚好类的情况下，Slingshot 更直接：

BiocManager::install("slingshot")

library(slingshot)
library(SingleCellExperiment)

sce <- as.SingleCellExperiment(seurat_obj)

sce <- slingshot(
  sce,
  clusterLabels = "seurat_clusters",
  reducedDim    = "UMAP",
  start.clus    = "0"
)

pseudotime <- slingPseudotime(sce)
head(pseudotime)

输出的 pseudotime 是一个 细胞 × 分支数 的矩阵，每列对应一条分支的拟时序值。可视化：

library(RColorBrewer)
colors <- colorRampPalette(brewer.pal(11, "Spectral")[-6])(100)

plot(
  reducedDims(sce)$UMAP,
  col  = colors[cut(pseudotime[, 1], breaks = 100)],
  pch  = 16,
  asp  = 1
)
lines(SlingshotDataSet(sce), lwd = 2, col = "black")

找随拟时序变化的基因走 tradeSeq：

library(tradeSeq)

sce       <- fitGAM(sce)
assocRes  <- associationTest(sce)
sig_genes <- rownames(assocRes)[assocRes$pvalue < 0.05]

scVelo：用 RNA velocity 给出方向

标准拟时序只回答"这些细胞的相对位置"，但不回答"他们正在往哪变"。RNA velocity 利用未剪接 mRNA 和已剪接 mRNA 的比例估计每个细胞的"即时变化方向"，在分化终点还没进入数据的项目里特别有用。

前提：有包含 spliced / unspliced counts 的 loom 文件。通常走 velocyto 生成：

velocyto run10x -m repeat_msk.gtf /path/to/sample01 genes.gtf

Python 分析：

import scanpy as sc
import scvelo as scv

adata = scv.read("sample01.loom", cache=True)
adata_seurat = sc.read_h5ad("analyzed.h5ad")
adata = scv.utils.merge(adata, adata_seurat)

scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)

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

需要更精细的估计时用动态模型（dynamical mode），慢但准：

scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)

scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color="velocity_pseudotime", cmap="gnuplot")

PAGA 和 CytoTRACE：快速概览

PAGA 把聚类结构抽成一张图，每个 cluster 是节点，节点之间边的粗细代表"有多少细胞邻居跨 cluster"：

sc.tl.paga(adata, groups="leiden")
sc.pl.paga(adata, color=["leiden", "CST3"])

sc.tl.draw_graph(adata, init_pos="paga")
sc.pl.draw_graph(adata, color="leiden", legend_loc="on data")

import numpy as np
adata.uns["iroot"] = np.flatnonzero(adata.obs["leiden"] == "0")[0]
sc.tl.dpt(adata)
sc.pl.umap(adata, color=["dpt_pseudotime"], color_map="viridis")

CytoTRACE 更简单，不建轨迹、只给每个细胞一个"分化程度估计值"：

devtools::install_github("digitalcytometry/cytotrace")

library(CytoTRACE)
results <- CytoTRACE(as.matrix(seurat_obj@assays$RNA@counts))
seurat_obj$CytoTRACE <- results$CytoTRACE

FeaturePlot(seurat_obj, features = "CytoTRACE")

CytoTRACE 的假设是"分化程度高的细胞表达基因更少"，对造血、上皮等几种经典系统比较可靠，放到任何系统上都用要谨慎。

结果解读的几个注意事项

• 拟时序值是相对的，不是真实时间。两个细胞 pt 差 2 和差 20 不代表真实时间差 10 倍。
• 起点错了整条轨迹都错。选 root 时如果没有分化早期的 marker，可以用 CytoTRACE 或 velocity 方向辅助判断。
• 分支点要交叉验证：Monocle3 的分支点和 marker 基因切换、已发表分化时序一致时才有说服力。
• 不要把轨迹当作因果。拟时序只表明"这些细胞在基因表达空间里排成了这样"，转化成 "A 变成 B" 这个因果陈述需要额外证据（比如谱系追踪实验或 velocity 方向）。

下一步

• 06 细胞通讯分析
• 07 多模态数据分析

参考资源

• Monocle3 文档
• Slingshot Bioconductor 页面
• scVelo 教程
• Saelens et al. 2019, 轨迹推断方法评测