Monocle拟时序分析：为何必须用原始counts而非SCT或整合数据？

上个月处理一批神经元分化的10x数据时，同组的师妹跑过来问我：Monocle做拟时序分析，到底该喂SCT整合后的数据，还是老老实实用RNA assay里的counts？她说自己用Seurat做了SCT整合，又用Harmony跑了integration，结果把整合后的表达矩阵丢给Monocle3，轨迹画出来乱成一团，跟我上一篇教程里的结果完全对不上。

这个问题其实问到了单细胞轨迹分析最关键的坑上。Monocle系列的输入数据选择，决定了你是拿到一条干净的发育轨迹，还是一堆看不出方向的细胞云。从Monocle2到Monocle3，从DDRTree到UMAP降维，很多教程都只说“把表达矩阵传进去”，却没人讲清楚这个矩阵到底应该是原始UMI counts、SCT变换后的counts，还是integration之后的矫正值。我当年也是踩了两次坑才弄明白：结论先说在前头——轨迹分析和拟时序推断，表达矩阵一定要用原始counts，而整合数据只能用来提供降维坐标或聚类标签。这篇就把原理、实操和排查全部拆开讲清楚，帮你少走弯路。

1. 为什么这个问题值得单独写一篇

1.1 一个卡了四天的问题

我的那个“四天”不是夸张。当时我在做一批肠类器官的分化数据，细胞从干细胞状态向吸收系和分泌系分化，理论上应该拉出一条Y字形轨迹。第一版我用Seurat里integration过后的assay跑Monocle3，具体操作是直接把seurat_obj@assays$integrated@data当作表达矩阵传入new_cell_data_set()。结果呢？learn_graph画出来之后，主轨迹直接从一群细胞里横穿过去，分支节点全堆在一端，秩序感的连续过渡完全看不出来。后来我又换成SCT的counts跑Monocle2，结果更惨，细胞被压成几个重叠的团，轨迹跟分群标签完全无关。

当时我以为是自己过滤参数没调好，试了不同num_dim、min_dispersion，折腾了三天毫无进展。第四天我决定回到最原始的RNA counts重跑一遍，轨迹一下子就恢复成了结构清晰的Y字形，分支点和Marker基因都严丝合缝。全局从头到尾，我只换了一个输入矩阵，其他参数完全没动。这才让我下决心把Monocle的数据输入底细彻底扒干净。

1.2 Monocle输入数据到底有什么硬性要求

先看硬性规定。Monocle3::new_cell_data_set()要求传入三个核心对象：expression_matrix、cell_metadata、gene_metadata。其中expression_matrix必须是基因在行、细胞在列的矩阵，而且官方推荐传入sparseMatrix，整数counts。这里的counts一般指UMI counts，也就是每个基因在每个细胞里被捕获到的转录本分子数。它不是标准化后的CPM、TPM，更不是log1p变换值，也不是批次校正后的残差。

为什么非要整数counts？因为Monocle3后续的很多步骤默认数据服从负二项分布或准泊松分布，整数counts恰好落在这些计数分布模型的支持域里。如果你丢进去一堆带小数的标准化值甚至负数，模型拟合时会出现各种怪异行为，最常见的表现就是graph_test()输出荒谬的Moran's I值，或者fit_models()卡住不收敛。

Monocle2的要求同样如此。它内部用DDRTree做降维时，表达矩阵会被用来计算细胞间距离和推断细胞状态转换。虽然它不是严格概率模型，但counts和log1p值产生的距离结构完全不同，直接用标准化值会导致DDRTree的树形拓扑被压缩，分支结构失真。所以无论你用Monocle2还是Monocle3，输入矩阵必须是原始counts，这是第一条铁律。

1.3 轨迹推断的算法逻辑决定了它吃不了“校正过”的值

光知道要求还不够，得理解为什么。轨迹分析的本质是从单细胞表达谱中重构细胞状态连续变化的路径。Monocle3的大致流程是：先对表达矩阵做PCA降维，再用UMAP把细胞投射到二维平面，随后通过聚类和learn_graph()在细胞云中学习主图结构，最后order_cells()找一个“根”细胞并沿着图计算拟时序。

关键在这几个环节全都在依赖表达量差异。PCA要最大化细胞间方差，UMAP保留局部与全局距离结构，learn_graph()寻找密度与距离定义的连通骨架。如果你喂进去的是SCT变换后的残差值，它已经把“测序深度”这个技术因素强行抹平，同时把高表达基因的离散度做了正则化。这个操作对聚类有好处，但对轨迹分析反而会掩盖真实表达差异——尤其当某些关键基因只是从低到高平缓变化时，SCT的压缩可能让这种变化变得不可见。

integration就更麻烦了。Harmony或Seurat整合的目标是让不同样本/不同批次中同类型细胞在降维空间里相互靠近，它本质上是在用“批次来源”这个信息去主动扭曲表达空间。轨迹重建需要的是保留真实的生物学变化梯度，而整合后的空间里，真实的发育差异可能被当成“批次差异”被抹除，或者反过来把批次差异残留成伪轨迹。所以直接用整合后的表达矩阵跑Monocle，相当于拿一张被人为PS过的地图去爬山，路线准不准全凭运气。

2. SCT与integration：确实很香，但不是给轨迹分析用的

2.1 SCT到底对数据做了什么

SCTransform是一个广受欢迎的标准化方法。它用正则化负二项回归对每个基因建模，把UMI counts中混入的测序深度因素回归掉，同时保留生物学变异。相比传统的LogNormalize方法，SCT能更好处理高深度样本的过度离散问题，也更适合多个样本合并后的下游分析。

但有一个事实很多人忽略：SCT之后的数据存储在assays$SCT@counts和assays$SCT@data里，其中@counts已经不再是原始UMI整数，而是经过变换后取整或重新归一化的表达值（实际是修正后的counts，通常带有小数；在部分Seurat版本中是整数化处理的近似），@data则是对应的log1p变换。也就是说，SCT counts本身就不等于原始分子计数，它已经混合了模型修正信息。这样的矩阵用来做聚类和差异表达没问题，但用来做Monocle的拟时序推断，会让算法在拟合计数分布时产生偏差。

我做过一个对照实验：同一批数据，一个输入原始RNA counts，一个输入SCT counts（注意这里还特意挑了Seurat内部round成整数的SCT counts），然后全部参数保持一致跑Monocle3。结果原始counts跑出的轨迹在已知分化节点上Marker基因呈现平滑波浪式变化，而SCT counts跑出的轨迹上，同一个Marker的表达模式变成了锯齿状，且分支点的位置偏移了至少两个cluster宽度。原因并不复杂：SCT的修正逻辑是等化技术噪音，但等化过程中也会顺带削弱基因表达的内在动态范围，而轨迹分析恰恰需要这份动态。

2.2 integration又是如何改变表达矩阵的

Seurat的整合流程，尤其是IntegrateData()，会基于共享的“锚点”对细胞表达谱进行相互修正，最终产出一个新的integrated assay。这个assay的@data是经过批次矫正的标准化表达值，它是为了跨样本聚类服务的，不是为了还原真实分子状态。Harmony稍有不同，它不直接生成新的表达矩阵，而是修正PCA空间坐标，但它同样会通过迭代混合的方式改变细胞在低维空间中的相对距离。

问题就出在这里。Monocle3的preprocess_cds()如果用了你传入的矩阵PCA，它的低维空间会和Seurat/Harmony产生的低维空间完全不同。如果你用integration后的数据作为表达矩阵，PCA轴被批次矫正主导，UMAP距离结构被强行拉平，轨迹图自然就失真了。哪怕你把Harmony的降维坐标硬塞给Monocle3（后面会讲方法），只要表达矩阵本身是整合矫正值，learn_graph()和拟时序计算时依然基于这套被扭曲的表达值，结果照样不可靠。

2.3 错误案例：拿SCT counts或integrated data跑Monocle会怎样

我见过几种典型症状，基本都是这么作出来的：

• 轨迹主干不明显，细胞变成一个大团，所有分支纠缠在中心，看不清连续路径；
• 轨迹的分支点特别多，且大量分支只连接两三个细胞，看起来像个毛发团；
• 拟时序值和已知分化方向相反，或者连续状态被压缩成三级跳；
• 用plot_genes_in_pseudotime()画Marker基因，表达曲线不是平滑平滑的，而是剧烈抖动或直接平台。

这些症状的核心成因都一样：Monocle在它认为的“细胞状态梯度”里没有找到足够的连续信号。整合或SCT后的表达矩阵已经部分抹平了梯度，算法只能乱抓局部噪音来构图。所以如果你现在跑出的Monocle结果怎么看怎么别扭，第一件事不是调参，而是先检查自己到底喂了什么矩阵进去。

3. 实操：到底怎么组合才是最佳实践

3.1 方案一：Monocle3全流程用原始counts（最稳妥）

最简单可靠的方式就是在Monocle3里全程使用原始RNA counts，不掺入任何整合信息。流程如下：

library(Seurat) library(monocle3) # 假设seurat_obj已经做过常规QC、LogNormalize和聚类 counts_matrix <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts") cell_meta <- seurat_obj@meta.data gene_meta <- data.frame(gene_short_name = rownames(counts_matrix)) rownames(gene_meta) <- rownames(counts_matrix) cds <- new_cell_data_set(counts_matrix, cell_metadata = cell_meta, gene_metadata = gene_meta) cds <- preprocess_cds(cds, num_dim = 50) cds <- reduce_dimension(cds, reduction_method = "UMAP") cds <- cluster_cells(cds) cds <- learn_graph(cds) cds <- order_cells(cds)

这样跑出来的轨迹只依赖原始的转录组成像关系，所有生物学梯度都来自真实的分子计数。缺点也明显：如果你分析的是多批次、多样本的数据，原始counts里混入了不容忽略的批次效应，细胞可能先按样本分成几个大群，轨迹被批次分隔成几个碎片。这时就需要方案二或方案三。

3.2 方案二：Seurat整合定cluster，Monocle3用counts重建（推荐）

这是我在多批次项目里最常用的方案。思路是：用Seurat做SCT整合和Harmony，因为它擅长把同类型细胞从批次效应中拉到一起，聚类分群很稳定；但进入Monocle3时，只取整合得到的聚类标签和Harmony降维坐标，表达矩阵仍然用原始RNA counts。

大概流程是：

# 1. Seurat标准流程 seurat_obj <- SCTransform(seurat_obj, vars.to.regress = "percent.mt") seurat_obj <- RunPCA(seurat_obj) seurat_obj <- RunHarmony(seurat_obj, group.by.vars = "sample") seurat_obj <- RunUMAP(seurat_obj, reduction = "harmony", dims = 1:30) seurat_obj <- FindNeighbors(seurat_obj, reduction = "harmony", dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = 0.8) # 2. 构建Monocle3 cds counts_matrix <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts") cell_meta <- seurat_obj@meta.data gene_meta <- data.frame(gene_short_name = rownames(counts_matrix)) rownames(gene_meta) <- rownames(counts_matrix) cds <- new_cell_data_set(counts_matrix, cell_metadata = cell_meta, gene_metadata = gene_meta) # 3. 用Seurat聚类标签初始化分区 cds <- preprocess_cds(cds, num_dim = 50) cds@clusters$UMAP <- seurat_obj@meta.data$seurat_clusters names(cds@clusters$UMAP) <- rownames(cell_meta) # 4. 手动设定降维坐标：直接用Seurat/Harmony的UMAP坐标 cds@reduce_dim_aux$UMAP <- list() cds@reduce_dim_aux$UMAP$model <- list(umap_coords = seurat_obj@reductions$umap@cell.embeddings) # 5. 继续走的流程 cds <- learn_graph(cds) cds <- order_cells(cds)

这里有两个细节要特别注意。第一，cds@clusters$UMAP是一个命名向量，名字必须是细胞barcode，值就是Seurat的cluster标签。这样cluster_cells(cds)会被我们用@clusters覆盖，learn_graph()会以这个分区结构去学习图。第二，cds@reduce_dim_aux$UMAP$model$umap_coords需要是一个矩阵，行名是细胞barcode，列是UMAP1、UMAP2。设好之后，Monocle的后续图构建和拟时序计算会基于Seurat整合出来的细胞空间，但表达矩阵始终是RNA counts，既控制了批次效应，又保住了真实的表达梯度。

这个方法我实测下来最稳。细胞分群来自整合后的共识，轨迹结构又由真实counts驱动，Marker基因在拟时序上的变化流畅连贯，分支点也能对应上已知的谱系决定基因。

3.3 方案三：用SCT整合但显式传入counts（硬核方案）

如果你比较早就做了SCT，所有下游分析和注释都基于SCT assay，不想推翻重来，那么可以在传入Monocle时显式指定使用哪套counts。

counts_matrix <- Seurat::GetAssayData(seurat_obj, assay = "SCT", slot = "counts")

等等——这里的SCT counts能不能直接用？我的建议是不到万不得已不要用。因为SCT counts虽然也取名叫counts，但它是经过模型修正后的值，不等于原始分子计数。如果你想尽量保留SCT优势，又不想丢掉原始counts轨迹信号，可以把SCT用于聚类、注释和挑选根细胞，轨迹表达矩阵只用RNA counts，也就是顺着方案二走。

所谓方案三，更多是一种折中：你已经用SCT counts做了所有分析，临时想快速接入Monocle，那只有接受SCT counts带来的潜在轨迹扭曲。操作层面没问题，代码也能跑，但结果解读时一定要谨慎，最好跟RNA counts的结果做交叉验证。

3.4 各方案对比表

方案二最符合“轨迹分析用原始counts、细胞聚类用整合信息”这一经验法则，也是我目前给所有人的默认推荐。

4. 从Seurat到Monocle3的完整代码和参数细节

4.1 数据准备与检查

构建Monocle3对象之前，必须先确认counts矩阵的质量。肉眼检查至少包括三项：

• 矩阵里有没有负数？counts矩阵不允许负数，有负数说明你拿错slot了。
• 矩阵是不是整数？虽然Monocle3对整数要求不是绝杀死，但非整数counts会让后续负二项拟合偏差变大。
• 行名和列名对不对得上？行名是基因名，列名是细胞barcode，顺序无所谓，但名字必须和metadata一致。

检查代码可以这样写：

# 检查counts矩阵是否为整数 expr_mat <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts") is_integer <- all(expr_mat@x == round(expr_mat@x)) cat("is integer counts:", is_integer, "\n") # 查看范围 range(expr_mat@x)

如果is_integer为FALSE，多半是用了SCT assay的counts或某些工具输出的小数矩阵。这时要么回头找原始RNA counts，要么用round()强行取整——但这不是好习惯，取整只能骗过检查，不能骗过统计模型。

4.2 降维坐标与聚类标签的传递

如果你选方案二，这里有一个很多人会踩的坑：Monocle3里的reduce_dimension()不是必须调的。你完全可以直接用Seurat的UMAP坐标，甚至用Harmony的PC坐标继续Monocle的preprocess_cds()和learn_graph()。

最稳妥的传参方式是：

cds <- preprocess_cds(cds, num_dim = 50) # 手动注入降维坐标 umap_coords <- seurat_obj@reductions$umap@cell.embeddings cds@reduce_dim_aux[["UMAP"]] <- list(model = list(umap_coords = umap_coords)) # 手动注入聚类信息 cds@clusters[["UMAP"]] <- seurat_obj@meta.data$seurat_clusters names(cds@clusters[["UMAP"]]) <- rownames(seurat_obj@meta.data) cds <- learn_graph(cds)

这里需要注意umap_coords的行名顺序。Monocle内部会通过细胞名匹配坐标，因此矩阵行名必须和colnames(cds)完全一致。如果不一致，learn_graph()会报错或者画出错乱图。建议传之前先做一次umap_coords <- umap_coords[colnames(cds), ]，强制对齐顺序。

另外，cds@clusters$UMAP这个名字不是随便起的。它对应reduce_dim_aux$UMAP，意思是在UMAP这个降维结果上做的聚类。Monocle里cluster_cells()运行后也会生成类似结构。如果你跳过cluster_cells()直接手动塞标签，记得cds@clusters$UMAP的名和cds@reduce_dim_aux$UMAP$model$umap_coords的行名保持一致。

4.3 learn_graph和order_cells参数说明

learn_graph()是Monocle3里最容易“黑箱”的一步。它要用cluster_cells()得到的分区信息去学习细胞图结构，核心参数有use_partition和close_loop。默认use_partition = TRUE，意思是算法会把不同的cluster partition当作独立的图来学习，这样有利于呈现分支结构；但如果你的数据里存在环状过渡（比如细胞周期），可以尝试close_loop = TRUE让首尾连接。

order_cells()需要你指定根节点。可以用root_cells参数直接给一个或几个细胞barcode，也可以用root_pr_nodes指定根节点。我一般会在learn_graph()之后调用plot_cells()，目测选一个处于最早分化状态的细胞群，然后把这群细胞中的一个作为根。如果你有明确的Marker基因，也可以用order_cells(cds, root_cells = cells)来固定根。

拟时序结果里，cds@principal_graph_aux[["UMAP"]]$pseudotime存储了每个细胞的拟时序值。后续plot_genes_in_pseudotime()或graph_test()都会用到它。注意，如果你在自己绘图时发现拟时序最大最小值只有零和一，很可能是传递根细胞的方式出错了，或者你传的数据源不对，不要让这种结果进入下游分析。

4.4 拟时序下游分析：graph_test与差异基因

拟时序算出来不等于结束。真正有生物学意义的结论来自“哪些基因随拟时序变化”。Monocle3提供了graph_test()，它基于空间自相关统计量Moran's I，找到在轨迹上表达模式非随机的基因。

gene_fits <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)

拿到结果后，重点看morans_test_stat和q_value，过小的q值乘以Moran统计量后基因往往就是关键状态转换因子。这个步骤同样强依赖表达矩阵的可靠性。如果用SCT或integration矩阵跑，Moran's I会被严重高估或低估，你的候选基因列表会失控。我见过有人跑出上千个显著基因，结果一半是线粒体基因和核糖体蛋白基因，这就是典型的输入矩阵有问题。

5. 常见问题与排查技巧实录

5.1 错误1：直接传SCT@counts

症状：轨迹分支乱，Marker基因在拟时序上的曲线不光滑。不少人看到“SCT assay也有一个counts slots”就想当然用了，但实际上这个counts不等于原始分子数。排查方式很简单：把SCT counts和RNA counts相加后看总量，SCT counts的矩阵总量会和原始UMI总量有明显差异，尤其高深度样本差异更大。正确做法永远是回退到RNA counts。

5.2 错误2：用integrated assay的@data

症状：细胞全部挤在一团，learn_graph()结果没有明显分支。integrated data是经过批次校正和中心化的数据，包含负值，PCA距离被压缩。如果已经跑了这步，赶紧换成RNA counts重跑。如果你担心批次效应，用方案二的Harmony坐标即可，不要整个矩阵替换。

5.3 怎么快速判断表达矩阵是不是counts

除了上面提到的整数检查，还有一个实用技巧：看矩阵的稀疏比例。原始UMI counts是高稀疏矩阵，绝大多数基因在绝大多数细胞里是0，稀疏率通常在90%以上。而log1p标准化后的data矩阵虽然0也多，但它已经失去了整数的特性；integrated data则可能产生很多负值。直接在R里运行summary(expr_mat@x)，如果min是0且max是几百到几万之间，大概率是counts；如果min是负数，绝对是标准化data。

5.4 多个样本到底该不该批次整合

如果你的数据包含多个样本或批次，单用raw counts跑Monocle常常会看到细胞先按样本分成几团，然后每团内部再沿着轨迹排列。这不是数据有问题，而是批次效应在counts层面没有被校正。解决方案就是在方案二里利用Harmony或Seurat整合坐标，把样本间的偏移在坐标层面纠正，同时保住counts梯度。注意，这里不要用vars.to.regress = "sample"去回归掉样本信息，那会把生物学差异也抹掉。最好在Seurat里用SCT整合或Harmony，只把低维坐标传给Monocle。

5.5 拟时序结果不受控制的排查思路

如果跑完order_cells()后，拟时序值在某个cluster内部突然断层，或者拟时序与已知Marker表达完全冲突，我建议按这个顺序排查：

• 先确认表达矩阵是raw counts，不是SCT/integrated；
• 再确认降维坐标是否来自整合后的UMAP，是的话有没有和counts矩阵匹配；
• 然后看cds@clusters$UMAP是否和colnames(cds)一致；
• 最后检查根细胞是不是选错了，换个根细胞或根节点重跑。

很多情况下，问题不是算法参数，而是前面数据传递时埋下的雷。

6. 最后想说的

我在实际项目里来回试过五六种组合，最终固定下来的干活模板就是方案二：Seurat负责整合和注释，Monocle3只拿RNA counts做轨迹构建，外部降维坐标负责给轨迹锚定细胞空间。这套组合省了我大量调参时间，也不用担心批次效应把轨迹撕得四分五裂。

另外一个特别想提醒的细节是：不要迷信“整合后数据更准”这句话。integration解决的是跨样本比较问题，轨迹分析解决的是细胞状态过渡问题，两者的数学目标并不一样。先想清楚你的生物学问题到底需要哪个工具去回答，再决定把什么数据喂进什么算法。拟时序分析本质上是在基因表达梯度里寻找方向，任何一步的过度校正都可能把这个梯度抹平。

如果你现在正被Monocle的轨迹问题折磨，第一件事去检查你传进去的表达矩阵到底是不是原始counts。不是的话，直接停机改数据，别在参数和图形美化上浪费时间。改完之后你会回来感谢这条建议的。