这个问题几乎是每个用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数据。SCTSCTransform是基于正则化负二项回归计算的皮尔逊残差。皮尔逊残差是连续值有正有负均值为0它的大小表示实际表达量偏离模型预期的程度。这个残差在找高变基因、做聚类时非常好用但和表达量本身完全是两回事。把残差喂给monocle等于让一个负二项模型去拟合一堆负数和小数统计假设直接崩掉。1.2 SCT和integration数据到底做了什么变换SCTransform的本质是对每个基因拟合一个正则化负二项回归模型然后取残差。这个残差矩阵就是SCTdata它代表的是扣除技术噪声和测序深度影响后的表达偏差而不是真实的分子计数。它甚至不是表达量而是一种标准化后的统计量。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对象里通常躺着这几个东西RNAcounts原始UMI counts整数矩阵这是monocle该吃的东西。RNAdatalog1p归一化后的表达值常用于绘图和找marker基因但不要喂给monocle。SCTcounts、SCTdataSCTransform之后的计数和残差矩阵。SCTcounts虽然也是整数但它不是原始计数已经经过了模型修正当作原始counts用也不合适。integrateddata整合缩放后的数据通常已经中心化存在负值。可能还有integratedscale.data这个是标准化后的矩阵和counts差别最大。很多教程只会告诉你怎么画图不会讲底层的对象结构导致新手根本不了解该取哪个slot。判断标准其实就一条看counts里存的是不是整数计数。凡是data或scale.data里的数据默认都不要直接喂给monocle。另外需要留意的是如果你手动给Seurat对象做过JoinLayers或者修改过assay取数前最好用dim()和class()确认一下矩阵类型别盲目复制网上的代码。2.2 正确提取counts矩阵并完成基因过滤这里我直接给出一段我自己流程里的标准代码。假设你已经有一个跑完标准Seurat流程的对象seurat_objlibrary(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的CellDataSetCDS对象关键是要同时准备细胞meta信息和基因注释信息# 准备细胞meta信息 metadata - seurat_objmeta.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_objreductions[[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_objreductions[[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数据会发生什么我试过把SCTdata皮尔逊残差矩阵直接作为表达矩阵喂给monocle2。第一反应是estimateDispersions直接报错报错信息大意是variance cannot be modeled之类的指向离散度拟合失败。这是因为残差矩阵有大量负值负二项模型的对数链接函数无法处理负的均值。即便你强行跳过了离散度估计后面differentialGeneTest的结果也会非常离谱——大量的基因q值趋近于0排序基因集中在一堆低表达基因上轨迹完全没有生物学意义。这是因为残差矩阵中基因的方差结构和原始counts完全不同差异基因检测的意义被扭曲了。SCTcounts这个夹在中间的选项我也试过。它虽然是整数但已经被SCTransform模型修正过不再代表真实的分子计数。用它构建的CDS轨迹整体方向和用原始counts的结果可能大体一致但分支长度和拟时序数值会偏。更重要的是如果你和别人共享数据或者发表文章评审问起来你喂给monocle的到底是什么你需要解释清楚SCTcounts的来源反而麻烦。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 我的一些经验性结论用一张表格总结不同数据输入的实测表现方便你对照参考输入数据是否推荐实测表现RNAcounts强烈推荐符合负二项假设轨迹稳定差异基因检测可靠RNAdata (log1p)不推荐非整数size factor估计困难结果无统计依据SCTdata残差不推荐含负值直接报错或结果异常SCTcounts不推荐非原始计数轨迹偏差且难解释integrateddata强烈不推荐负值多报错频繁轨迹方向可能反转这里补充几个我自己总结的避坑经验第一多数报错不是monocle不好用而是输入数据格式不对。遇到Error in log(x)、variance model failed、size factor estimate failed这类信息要养成条件反射先查自己的表达矩阵是不是counts。第二estimateDispersions非常耗时尤其当基因数很多时。我通常会把表达矩阵先过滤到5000个基因以内再构建CDS。如果过滤后仍然很慢可以考虑用cores参数做并行计算也能显著提速。第三monocle3的preprocess_cds内部会跑PCA如果你的counts矩阵里有基因在所有细胞中都是0PCA会给出警告。事前过滤低表达基因不仅为了速度也为了数值稳定性。第四拟时序分析完成后一定要用已知的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的主食。希望这篇能把你的纠结一次性解决让轨迹分析少走弯路。