拟时序分析：Monocle / Slingshot

跟完这篇，您能把一份已经聚好类的单细胞数据，用 Slingshot 或 Monocle3 排出一条「细胞状态渐变」的轨迹，读懂它的起点、分支和沿轨迹变化的基因，并且清楚地知道：拟时序推的是相对顺序，不是日历时间——什么时候能信、什么时候只是工具替您画了一条线。

拟时序（pseudotime）解决的是一个具体问题：一群细胞其实处在同一个连续过程的不同阶段（比如从祖细胞逐步分化成熟），但它们是在同一时刻一次性测下来的。拟时序就是把这些「快照里被打散的时间」重新排回一条顺序上。下面先讲清它在推什么，再分别用两个主流工具跑一遍。

拟时序在推断什么

先建立一个准确的直觉，避免后面过度解读。

它推的是顺序，不是时间。 拟时序假设：转录组越相似的细胞，在过程中的位置越接近。于是它把细胞沿「表达渐变」的方向串起来，给每个细胞一个数值——这个数值只表示「谁在前、谁在后」，单位是任意的，不对应小时、天。
它只在「有连续过程」时才成立。 分化、激活、应激反应这类渐变过程，细胞之间存在中间态，拟时序才有意义。如果您手上是几群彼此独立、稳定的成熟细胞类型，它们之间根本没有过渡，硬跑轨迹工具照样会画出一条线，但那条线没有生物学含义。
方向要您来定。 算法给出的轨迹本身是无向的——它知道细胞排成一条链，但不知道哪头是开始。哪端是起点，得靠生物学先验指定（见后文「怎么定起点」）。

一句话：拟时序是用一个时间点上的细胞异质性，反推一个过程的先后。它是「排序」，不是「计时」。

开始之前：适不适合做拟时序

动手前先确认三件事，省得在不该做轨迹的数据上白忙一场。

数据已经聚类、最好也注释过。 拟时序是标准单细胞流程的下游步骤，前面的质控、降维、聚类、注释要先跑通。这条主线见单细胞分析全流程，注释见细胞注释。
确实存在一个连续过程。 问自己：这些细胞在生物学上是不是一个会「你变成我、我变成他」的渐变体系？UMAP 上是不是一条相连的「带状」结构，而不是彼此孤立的几团？如果是几团离散的稳定类型，别做拟时序。
起点心里有数。 您得能说出这个过程从哪种状态开始（哪群是祖细胞、或哪个采样时间点最早）。定不出起点，轨迹的方向就是猜的。

下面两个工具，Slingshot 和 Monocle3，思路一致、接口不同，选一个即可。

	Slingshot	Monocle3
输入	降维坐标 + 聚类标签	`cell_data_set` 对象
轨迹模型	簇质心最小生成树 + 主曲线	主图（principal graph）
分支	多条 lineage（谱系）	图上的多个分叉
定起点	`start.clus` 指定起始簇	`root_cells` / 交互选点
找变化基因	配 tradeSeq	自带 `graph_test`

用 Slingshot 推一条轨迹

Slingshot 的做法很直观：先在聚类质心之间连一棵最小生成树，定出「谁接谁」，再穿过这些簇拟合平滑的主曲线，把每个细胞投影到曲线上得到拟时序。它天然支持分叉——一棵树有多个末端，就有多条谱系。

library(slingshot)
library(SingleCellExperiment)

# sce 已带 UMAP 坐标(reducedDims 里名为 "UMAP")与聚类标签(colData 的 cluster 列)
sce <- slingshot(
  sce,
  clusterLabels = "cluster",   # 用哪一列聚类标签来连轨迹
  reducedDim    = "UMAP",       # 在哪个降维空间上拟合曲线
  start.clus    = "3"           # 起点簇：由生物学先验指定,见下文
)

# 每个细胞的拟时序(每条谱系一列;不在该谱系上的细胞为 NA)
pt <- slingPseudotime(sce)

# 各条谱系依次经过哪些簇
slingLineages(sce)

几点说明：

start.clus 是可选的。不给，Slingshot 也能连出无向的树，但方向随机——要得到有意义的先后，就得指定起点。
输出的 slingPseudotime(sce) 是一个矩阵：有几条谱系就有几列。同一个细胞若同时属于两条谱系（分叉前的共享段），这两列的值在分叉前是一样的。
想找沿轨迹显著变化的基因，Slingshot 官方推荐搭配 tradeSeq：对每个基因沿拟时序拟合广义可加模型（GAM），再检验它是否随拟时序变化。

library(tradeSeq)

# 用 slingshot 的谱系结构拟合每个基因的 GAM
sce <- fitGAM(counts = counts, sds = SlingshotDataSet(sce))

# 检验:表达是否随拟时序显著变化
assoc <- associationTest(sce)

用 Monocle3 推一条轨迹

Monocle3 把整条流程收进一个 cell_data_set 对象，从建对象、降维、聚类，到「学出一张主图」，再到「指定起点、算拟时序」，接口连贯。

library(monocle3)

# 1) 建对象:gene_metadata 须含 gene_short_name 列
cds <- new_cell_data_set(
  expression_data = counts,      # 基因 × 细胞 计数矩阵
  cell_metadata   = cell_meta,   # 每个细胞一行
  gene_metadata   = gene_meta    # 每个基因一行
)

# 2) 标准前处理:PCA -> UMAP -> 分群
cds <- preprocess_cds(cds, num_dim = 50)   # 主成分个数,按数据调
cds <- reduce_dimension(cds)               # 默认 UMAP
cds <- cluster_cells(cds)                  # 分群并划分 partition

# 3) 学出主图:这是轨迹的骨架
cds <- learn_graph(cds)

plot_cells(cds, color_cells_by = "cluster")

learn_graph 只画出「骨架」，此时还没有方向。指定起点后才有拟时序：

# 指定起点:用已知的祖细胞/起始状态的细胞作 root
cds <- order_cells(cds, root_cells = root_ids)

# 现在每个细胞有了拟时序,画出来看
plot_cells(cds, color_cells_by = "pseudotime")

# 取出拟时序数值
pt <- pseudotime(cds)

找沿轨迹变化的基因，Monocle3 自带 graph_test，基于主图上的 Moran's I（空间自相关）打分：

# neighbor_graph = "principal_graph" 表示沿轨迹找变化基因
deg <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)

# q_value 是校正后的显著性;morans_I 越接近 1,越是沿轨迹强烈、连续地变化
sig <- subset(deg, q_value < 0.05)

q_value < 0.05 是常用惯例、非铁律，可按需收紧；结合 morans_I 的大小一起看更稳，只卡 q 值容易把变化很弱的基因也放进来。
num_dim（主成分个数）没有标准答案，按数据的复杂度调，思路和普通单细胞降维一致。

怎么定起点

这是整个分析里最容易被忽略、又最影响结论的一步。请记住一句话：起点是生物学决定，不是算法算出来的。

算法给出的轨迹是无向的骨架，把哪一端设成「开始」，会整条翻转拟时序——原本「早期」的细胞变成「晚期」，「上调」的基因读成「下调」。定错起点，后面所有「先后」「早晚」的结论全反。所以起点必须有据可依，而不是随手指一个。常用的三种依据：

依据	怎么用	适用场景
已知起始态标志基因	找出该体系公认的祖细胞/干性标志基因表达最高的那群，设为起点	分化体系有明确的祖细胞 marker
采样时间点	若实验本身带时间标签（如第 0 天样本），把最早时间点的细胞设为起点	有时间序列设计的实验
过程方向的领域知识	依据这个生物学过程公认的方向（比如从未成熟到成熟）定端点	方向在领域内是共识

落到代码上，就是把上表选出的那群细胞的编号（或那个簇）交给工具：Slingshot 用 start.clus = "..."，Monocle3 用 order_cells(cds, root_cells = ...)（Monocle3 也支持在图上交互点选起点节点）。

定不出可靠起点时，诚实的做法是明说方向未定，而不是随便挑一端然后把结论讲得斩钉截铁。

轨迹与分支怎么解读

拿到带方向的轨迹后，这样读它：

一条路径 = 一个命运。 从起点走到某个末端的每一条路径，代表一条分化/演变谱系。有几个末端，就有几条谱系。
分支点 = 命运分岔。 轨迹分叉的地方，通常对应一个祖细胞群「往两个方向分头走」的决策点。贴着分支点的细胞是过渡态，介于两种命运之间，别硬把它归成某一类。
沿轨迹变化的基因讲出了故事。 用前面 graph_test（Monocle3）或 associationTest（tradeSeq）挑出的基因，看它们在拟时序上何时开、何时关：起点高、随拟时序下降的，多是干性/祖细胞程序；末端才升起来的，多是成熟/功能基因。分支两侧各自特异开启的基因，则勾勒出两种命运的分野。
末端 = 相对稳定的终态。 谱系末端一般对应分化完成、较稳定的成熟细胞状态。

务必警惕：分支可能是假的。 聚类分辨率、降维参数、离群细胞都可能凭空造出一个分叉。判断一个分支真不真，回到生物学——分叉两侧是不是真有各自的 marker、是不是符合已知的命运图谱。一条画得出来的分支，不等于一条真实存在的分支。marker 怎么核对见细胞注释。

拟时序 ≠ 真实时间

这一节请务必看完，它决定了您的结论能讲多硬。拟时序和真实时间的差别，不是细节，是原则：

它是相对排序，单位任意。 拟时序值只表达先后，不对应真实时长。「拟时序 5」和「拟时序 10」之间，不代表隔了固定的时间，更不代表隔了 5 个单位的真实时间。
跨谱系不可直接比。 不同谱系的拟时序各自成立、量纲不同。谱系 A 的拟时序 10 和谱系 B 的拟时序 10 不是同一件事，别把它们摆在一起比早晚。
距离不等比于耗时。 细胞可能在某个状态里停留很久，也可能很快掠过另一段。拟时序上走得远，不代表真实世界里花的时间多。
它默认「过程是连续且被采到的」。 拟时序靠一个时间点上的细胞异质性反推先后，前提是中间态确实存在并被测到。若过程不连续、或关键中间态没采到，排出的顺序就不可信。
方向是您赋予的。 前面反复强调：方向来自您指定的起点，不是数据自己「说」出来的。想要数据驱动的方向线索，可以了解 RNA velocity（用剪接/未剪接读段推断细胞「正朝哪走」）——它和拟时序互补，但也有各自的假设与局限，不是绝对真理。

一句话收尾：拟时序是一把排序的尺，不是一只计时的表。把它当排序用，结论稳；当计时用，容易过度解读。

在平台上跑一遍

上面这两条路线，在「百沐一下」不写代码也能走完：上传聚好类（最好已注释）的单细胞对象，用一句话说清目的——例如「对这份分化数据做拟时序，起点是 XX 祖细胞群，用 Slingshot 排轨迹并找出沿拟时序变化的基因」，组学智能体会自动完成建轨迹、按您指定的起点定方向、算拟时序、挑变化基因，并连同轨迹图、基因表达随拟时序的曲线和可复现的 R 代码一起交付。用法见组学智能体使用指南，或直接到 app.baimuyixia.com 试一遍。