企业官网建设流程全解析

这个问题几乎是每个用Seurat做完单细胞分析、想接着用monocle做轨迹的人都会卡一下的地方。我在各种生信交流群里见过太多次类似的提问："我用Seurat做了整合和聚类，能不能直接把integrated assay或者SCT assay喂给monocle？"底下的回答五花八门，有说可以的，有说不行的，还有人建议把SCT的残差转回整数再喂进去。作为一个把单细胞轨迹分析跑过上百个数据集的人，我可以直接告诉你结论：monocle的轨迹分析和拟时序分析，老老实实用RNA assay里的counts数据。这个结论背后不是个人习惯问题，而是算法模型的硬性要求。这篇就把道理讲清楚，顺便把我常用的完整流程和踩过的坑都整理出来，给正在为"到底该喂哪个数据"发愁的朋友一个明确参考。

1. monocle的模型假设决定了为什么不能用归一化后的数据

很多人不理解为什么一个简单的"数据输入选择"会引发这么多争论。根源在于monocle不是那种"随便给个表达矩阵就能跑"的工具，它内部的统计模型对数据结构有严格假设。你用错误的数据喂进去，轻则结果诡异，重则直接报错让你怀疑人生。

1.1 monocle2和monocle3背后的统计模型

先看monocle2。它的核心是使用VGAM包拟合基因表达随时间变化的广义线性模型，本质上假设基因表达量服从负二项分布（Negative Binomial）。负二项分布是针对非负整数计数数据的分布，原始UMI counts正好符合。在整个拟时序分析流程中，estimateSizeFactors和estimateDispersions这两个步骤都在为负二项模型的参数做估计，如果你塞进去的是连续值、有负数、有小数的归一化数据，这些估计就会出问题，后面所有的轨迹推断都是无源之水。

monocle3虽然改成了"先降维聚类，再学习轨迹图，最后排拟时序"的框架，看起来和monocle2差别很大，但它在估计基因表达随轨迹变化时（graph_test），底层使用的依然是基于负二项分布或准泊松模型的统计检验。所以不管你是用monocle2还是monocle3，表达矩阵层面的要求是一致的：必须是counts数据。

SCT（SCTransform）是基于正则化负二项回归计算的皮尔逊残差。皮尔逊残差是连续值，有正有负，均值为0，它的大小表示"实际表达量偏离模型预期的程度"。这个残差在找高变基因、做聚类时非常好用，但和"表达量本身"完全是两回事。把残差喂给monocle，等于让一个负二项模型去拟合一堆负数和小数，统计假设直接崩掉。

1.2 SCT和integration数据到底做了什么变换

SCTransform的本质是：对每个基因拟合一个正则化负二项回归模型，然后取残差。这个残差矩阵就是SCT@data，它代表的是"扣除技术噪声和测序深度影响后的表达偏差"，而不是真实的分子计数。它甚至不是表达量，而是一种标准化后的统计量。

integration（整合）后的数据就更复杂了。Seurat的integrated assay，是通过CCA、Harmony、MNN等方法做批次校正后，再进行缩放和中心化处理的结果。比如Harmony是直接在低维空间里做的校正，Seurat的整合流程还会把结果投影回基因空间。这些处理之后的矩阵值域已经彻底偏离原始counts，通常包含负值，分布形态也完全不同。用这种数据去套monocle的负二项模型，基本等于让鱼爬树。

如果扩展一下，还有人不死心，问log1p(counts)行不行。log转换后的数据虽然保留了表达的趋势，但同样不是计数数据。monocle官方文档对expressionFamily有明确说明：对于UMI counts，通常用negbinomial.size()；对于Smart-seq2这类非UMI全长转录组数据，可以选Tobit()或negbinomial()。无论哪种，前提都是底层矩阵是counts。你非要用log数据，除非手动指定expressionFamily = gaussianff()，但这就相当于抛弃了monocle最核心的统计推断基础，结果很难有说服力，不建议这么做。

2. 数据搬运前的关键准备：从Seurat对象到monocle对象

做单细胞分析，大部分人已经习惯用Seurat完成前期的质控、归一化、整合和聚类注释。到了轨迹分析这一步，很多人的困惑不是"不知道要用counts"，而是"不知道去哪里拿counts"以及"怎么拿才不破坏后面的分析流程"。这节把Seurat对象结构掰开来讲清楚。

2.1 先把Seurat对象的几个assay搞清楚

一个标准流程跑完之后，你的Seurat对象里通常躺着这几个东西：

• RNA@counts：原始UMI counts，整数矩阵，这是monocle该吃的东西。
• RNA@data：log1p归一化后的表达值，常用于绘图和找marker基因，但不要喂给monocle。
• SCT@counts、SCT@data：SCTransform之后的计数和残差矩阵。SCT@counts虽然也是整数，但它不是原始计数，已经经过了模型修正，当作原始counts用也不合适。
• integrated@data：整合缩放后的数据，通常已经中心化，存在负值。
• 可能还有integrated@scale.data，这个是标准化后的矩阵，和counts差别最大。

很多教程只会告诉你怎么画图，不会讲底层的对象结构，导致新手根本不了解该取哪个slot。判断标准其实就一条：看@counts里存的是不是整数计数。凡是@data或@scale.data里的数据，默认都不要直接喂给monocle。另外需要留意的是，如果你手动给Seurat对象做过JoinLayers或者修改过assay，取数前最好用dim()和class()确认一下矩阵类型，别盲目复制网上的代码。

2.2 正确提取counts矩阵并完成基因过滤

这里我直接给出一段我自己流程里的标准代码。假设你已经有一个跑完标准Seurat流程的对象seurat_obj：

library(Seurat) library(monocle) library(monocle3) # 提取RNA assay中的原始counts expr_matrix <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts") expr_matrix <- as.matrix(expr_matrix) # 基础过滤：至少在5个细胞中检测到表达 cell_counts <- Matrix::rowSums(expr_matrix > 0) genes_keep <- names(cell_counts[cell_counts >= 5]) expr_matrix <- expr_matrix[genes_keep, ] # 如果你有SCT结果，优先从SCT里取高变基因 hvgs <- VariableFeatures(seurat_obj, assay = "SCT") if (length(hvgs) == 0) { hvgs <- VariableFeatures(seurat_obj, assay = "RNA") }

注意，as.matrix这一步在大数据集上很消耗内存，如果你的数据特别大，可以在构建monocle对象之前先过滤基因数，再转成稠密矩阵。我遇到过不少人在这一步被内存卡死，明明100G内存的机器也吃不住几万个细胞的全量counts转矩阵。所以务必要先做基因过滤，常见的做法是保留高变基因或者找到的排序基因的并集，再传给newCellDataSet。

2.3 SCT高变基因才是你真正需要从SCT结果里借的东西

这里要强调一个很多人没意识到的关键点：虽然SCT数据本身不能直接喂给monocle，但SCT识别出来的高变基因列表，却是轨迹分析里最值得利用的资源。

SCTransform在拟合残差时，会对基因的均值-方差关系做正则化建模，因此它筛出来的高变基因比单纯用FindVariableFeatures(selection.method = "vst")更稳健，受测序深度干扰更小。我在实际项目中，通常会把SCT的高变基因作为排序基因的候选池，然后结合differentialGeneTest或者monocle3的graph_test进一步筛选。

一个经验是：直接用5000个高变基因喂给monocle往往得不到清晰轨迹，因为高变基因里混杂了很多和细胞状态转变无关的基因（比如细胞周期基因、应激基因）。更好的做法是先跑一遍differentialGeneTest，用细胞聚类标签作为fullModel公式项，选出随细胞类型变化显著的基因，再从中取top 500~1000个作为排序基因。这样做轨迹会更干净，分支也更容易解释。

3. 一套可以直接复现的Seurat到monocle拟时序流程

讲完理论，给出一套可复现的操作流程。我以monocle2为主要例子，因为目前很多做传统拟时序分析的人还在用它；monocle3的关键差异我会在流程末尾单独说明。

3.1 构建monocle2的CDS对象

从counts矩阵构建monocle2的CellDataSet（CDS）对象，关键是要同时准备细胞meta信息和基因注释信息：

# 准备细胞meta信息 metadata <- seurat_obj@meta.data pd <- new("AnnotatedDataFrame", data = metadata) # 准备基因注释信息 gene_info <- data.frame( gene_short_name = rownames(expr_matrix), row.names = rownames(expr_matrix) ) fd <- new("AnnotatedDataFrame", data = gene_info) # 构建CDS对象 cds <- newCellDataSet( expr_matrix, phenoData = pd, featureData = fd, lowerDetectionLimit = 0.1, expressionFamily = negbinomial.size() ) # 估计size factor和离散度 cds <- estimateSizeFactors(cds) cds <- estimateDispersions(cds)

关于expressionFamily的选择，多说一句：对于UMI数据，推荐negbinomial.size()，因为它假设size factor已经通过UMI总数校正过，计算更快更稳定。如果你的数据是Smart-seq2这类全长转录组，不具备UMI特性，那用negbinomial()或Tobit()更合适。别小看这一步，选错family会导致后续estimateDispersions报错或者结果严重偏差。

另外，lowerDetectionLimit这个参数很多人不理解。它表示"被判定为检测到表达的最低值"，通常设置为0.1，主要是为了过滤掉那些在所有细胞中几乎不表达的基因。对UMI数据来说，0和1的区别很关键，设置一个略高于0的阈值有助于稳定后续的模型拟合。

3.2 挑选用于轨迹的排序基因

轨迹分析的核心是"哪些基因的变化最能定义细胞状态转变的顺序"。这一步直接影响轨迹的形状和生物学解释，不能马虎。

# 使用差异基因检测来选择排序基因 diff_test <- differentialGeneTest( cds, fullModelFormulaStr = "~Cluster", reducedModelFormulaStr = "~1" ) # 取top 800个显著差异基因 ordering_genes <- diff_test[order(diff_test$qval), ]$gene_short_name[1:800] cds <- setOrderingFilter(cds, ordering_genes) # 这一步完成后可以可视化看这些基因的表现 plot_ordering_genes(cds)

关于fullModelFormulaStr，这里有几个常见变体。~Cluster表示"基因表达随聚类变化"；如果你有明确的连续型协变量，比如时间点，可以写成~TimePoint；如果你想同时考虑多个因素，可以写成~Cluster + Batch。关键在于reducedModelFormulaStr，它代表零模型，通常写成~1即可。差异基因检测比较的是全模型和零模型的拟合差异，q值越小说明基因表达越受关注因素影响。

如果你用的是monocle3，排序基因的挑选过程不太一样，monocle3用graph_test来做基因模块分析，它可以找到在轨迹图上具有空间自相关性的基因。流程是在learn_graph之后跑graph_test(cds, neighbor_graph = "knn")，然后取q值最小的基因构建gene_modules。两种工具的思路不同，但都能达到"找到驱动轨迹的基因"这个目的。

3.3 降维、排序与可视化

monocle2默认使用DDRTree算法降维，将高维基因表达空间压缩到低维流形上，然后根据细胞在该流形上的位置推断拟时序。核心代码：

# 降维 cds <- reduceDimension(cds, max_components = 2, method = "DDRTree") # 推断拟时序 cds <- orderCells(cds) # 可视化 plot_cell_trajectory(cds, color_by = "seurat_clusters") plot_cell_trajectory(cds, color_by = "Pseudotime")

orderCells默认会选择一个细胞数量最多的节点作为根节点。但这个默认根节点未必是你生物学上想要的起点，比如你在研究T细胞分化，默认根点可能是中央记忆T细胞，但你更希望从naive T细胞开始算。这时候就要手动指定根节点：

# 先看轨迹图，确定根节点的state table(pData(cds)$State) # 指定根节点state cds <- orderCells(cds, root_state = 3)

这个操作在实操中非常常见，因为拟时序的"零点"决定了你如何解读分化方向。务必在出图之后检查根节点选择是否合理，必要时结合marker基因的表达来辅助判断。

如果你用的是monocle3，流程更现代一些：

cds <- preprocess_cds(cds, num_dim = 50) cds <- reduce_dimension(cds) cds <- cluster_cells(cds) cds <- learn_graph(cds) cds <- order_cells(cds) plot_cells(cds, color_cells_by = "pseudotime")

monocle3的优势在于它直接支持UMAP降维，轨迹更符合直觉；劣势在于它对数据的过滤要求更严格，而且preprocess_cds内部也是用PCA，如果你喂给它一个没有经过合理基因过滤的counts矩阵，后面的轨迹同样会乱。

4. 整合数据的正确打开方式：批次校正另辟蹊径

既然SCT和integration数据不能直接喂给monocle，那单细胞数据常见的批次效应问题怎么解决？总不能完全不管吧。这节讲清楚批次校正的替代方案，以及什么时候真的有必要用SCT/integration。

4.1 monocle3的align_cds批次对齐

monocle3提供了一个专门应对批次效应的函数align_cds。它的工作原理是在轨迹构建过程中，对指定的扰动因素（比如样本批次、患者ID）做线性回归，把该因素带来的表达变异从数据中扣除，同时保留counts的整体结构。这种方法不会破坏负二项分布的假设，效果也相当不错。

# 在preprocess之后、learn_graph之前调用 cds <- preprocess_cds(cds, num_dim = 50) cds <- align_cds(cds, alignment_group = "batch") cds <- reduce_dimension(cds)

注意align_cds的alignment_group参数只能接受一个因子型变量，如果你的数据存在多个批次因素（比如平台+患者），建议把它们合并成一列再传入。这个函数使用的对齐逻辑和Seurat的整合有本质区别，它是直接对表达矩阵的线性模型残差做处理，而不是重新构建一个嵌入空间，所以更适合作为monocle内部的批次校正手段。

4.2 借用Seurat整合嵌入但保留counts的混合方案

还有一个非常实用的混合方案，我实际项目中经常使用。思路很简单：你可以把Seurat基于integration算出来的UMAP坐标直接覆盖到monocle3的CDS对象上，让monocle3使用这个整合后的低维嵌入来构建轨迹。这样既利用了Seurat整合消除批次效应的优势，又保证了表达矩阵层面仍然是counts数据。

cds <- new_cell_data_set( expr_matrix, cell_metadata = metadata, gene_metadata = gene_info ) # 先用counts做常规预处理 cds <- preprocess_cds(cds, num_dim = 50) # 用Seurat整合之后的UMAP替换monocle3默认的UMAP reducedDims(cds)[["UMAP"]] <- seurat_obj@reductions[["umap"]]@cell.embeddings # 后续正常走monocle3流程 cds <- cluster_cells(cds) cds <- learn_graph(cds) cds <- order_cells(cds)

这个方案的好处是：UMAP嵌入里包含了跨样本/跨批次的整合信息，learn_graph在UMAP上构建的轨迹自然会把这些信息带进来。但你必须理解一点——你只是借用了坐标，表达矩阵本身没有变，所以后续的graph_test、基因模块分析都还是基于counts模型完成的，统计性质是安全的。

这个技巧有个前提：Seurat的UMAP细胞顺序必须和CDS的细胞顺序一一对应。实际操作中要先确认colnames(expr_matrix)和rownames(seurat_obj@reductions[["umap"]]@cell.embeddings)完全一致，顺序不一致的用match函数重排一下，否则会出大问题，而且这种问题往往非常隐蔽，报错都不一定看得明白。

4.3 什么时候真的需要SCT/integration

说到底，SCT和integration的核心应用场景是聚类和差异表达分析。它们能有效消除测序深度、批次效应等噪声，让细胞亚群的分辨更准确。当你面对多个样本、多批次数据时，先用Seurat做整合得到正确的细胞分群，再在注释结果的基础上继续轨迹分析，这是最合理的整体流程。

但这不意味着你要把SCT或integrated数据直接丢给monocle。更合理的做法是：把整合结果作为"分群依据"和"嵌入坐标"，把counts数据作为"统计建模依据"。分析链路可以概括为：原始表达矩阵 → Seurat/SCT整合 → 获得细胞类型注释和高变基因 → 从raw counts重新构建monocle对象 → 用SCT高变基因或差异基因作为排序基因 → 做轨迹和拟时序。

如果批次效应严重到影响轨迹推断，优先考虑align_cds或混合嵌入方案，而不是硬把归一化数据塞进monocle。这个原则你掌握了，大部分轨迹分析项目都不会翻车。

5. 实测踩坑记录与常见报错处理

最后把自己实际测试过的几个"错误输入"结果分享出来，给那些非想试试的朋友省点时间。以下都是我在真实数据集上跑出来的现象，不是凭空猜测。

5.1 直接喂SCT数据会发生什么

我试过把SCT@data（皮尔逊残差矩阵）直接作为表达矩阵喂给monocle2。第一反应是estimateDispersions直接报错，报错信息大意是"variance cannot be modeled"之类的，指向离散度拟合失败。这是因为残差矩阵有大量负值，负二项模型的对数链接函数无法处理负的均值。

即便你强行跳过了离散度估计，后面differentialGeneTest的结果也会非常离谱——大量的基因q值趋近于0，排序基因集中在一堆低表达基因上，轨迹完全没有生物学意义。这是因为残差矩阵中基因的方差结构和原始counts完全不同，差异基因检测的意义被扭曲了。

SCT@counts这个夹在中间的选项我也试过。它虽然是整数，但已经被SCTransform模型修正过，不再代表真实的分子计数。用它构建的CDS，轨迹整体方向和用原始counts的结果可能大体一致，但分支长度和拟时序数值会偏。更重要的是，如果你和别人共享数据或者发表文章，评审问起来"你喂给monocle的到底是什么"，你需要解释清楚SCT@counts的来源，反而麻烦。

5.2 直接喂integrated数据会发生什么

integrated数据的问题更明显。Seurat的integrated assay在整合后通常经过ScaleData处理，矩阵里大量负值，且值域分布完全不是计数形态。负二项模型直接罢工，我遇到过的报错是：

Error in log(x) : NaNs produced Error in calcVarianceModel(...) : variance model failed

如果你用的不是monocle2而是monocle3，预处理阶段可能没有立刻报错，因为preprocess_cds内部做PCA时并不严格要求非负。但等到learn_graph和graph_test阶段，结果就会开始飘。我在一个测试集上看到，用integrated数据跑出来的轨迹，分支方向和已知的细胞分化顺序完全反了。这就是统计假设被破坏后累积的误差，不是调参能救回来的。

5.3 我的一些经验性结论

用一张表格总结不同数据输入的实测表现，方便你对照参考：

这里补充几个我自己总结的避坑经验：

第一，多数报错不是"monocle不好用"，而是输入数据格式不对。遇到"Error in log(x)"、"variance model failed"、"size factor estimate failed"这类信息，要养成条件反射：先查自己的表达矩阵是不是counts。

第二，estimateDispersions非常耗时，尤其当基因数很多时。我通常会把表达矩阵先过滤到5000个基因以内再构建CDS。如果过滤后仍然很慢，可以考虑用cores参数做并行计算，也能显著提速。

第三，monocle3的preprocess_cds内部会跑PCA，如果你的counts矩阵里有基因在所有细胞中都是0，PCA会给出警告。事前过滤低表达基因不仅为了速度，也为了数值稳定性。

第四，拟时序分析完成后，一定要用已知的marker基因验证轨迹方向。我习惯把几个关键marker基因的表达叠加到轨迹图上，比如plot_cell_trajectory(cds, color_by = "MyMarkerGene")。如果marker的表达变化和已知生物学方向矛盾，首先怀疑根节点选择，其次怀疑排序基因列表是否混入了无关基因，而不是马上怀疑monocle本身。

就我个人习惯来说，现在跑单细胞轨迹分析基本固定两种路线：要么直接用原始counts喂monocle3，再用align_cds处理批次；要么先做Seurat整合，但只借用整合后的UMAP坐标，表达矩阵始终用counts。这个原则我实践了很久，很少翻车。如果非要对标题里的问题给一个最直接的答案：monocle的拟时序分析，用RNA assay里的counts数据。SCT和integration的结果可以用来挑基因、用来定坐标，但它们不是monocle的"主食"。希望这篇能把你的纠结一次性解决，让轨迹分析少走弯路。