上个月处理一批神经元分化的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_objassays$integrateddata当作表达矩阵传入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()输出荒谬的Morans 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$SCTcounts和assays$SCTdata里其中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_objmeta.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整合定clusterMonocle3用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_objmeta.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) cdsclusters$UMAP - seurat_objmeta.data$seurat_clusters names(cdsclusters$UMAP) - rownames(cell_meta) # 4. 手动设定降维坐标直接用Seurat/Harmony的UMAP坐标 cdsreduce_dim_aux$UMAP - list() cdsreduce_dim_aux$UMAP$model - list(umap_coords seurat_objreductions$umapcell.embeddings) # 5. 继续走的流程 cds - learn_graph(cds) cds - order_cells(cds)这里有两个细节要特别注意。第一cdsclusters$UMAP是一个命名向量名字必须是细胞barcode值就是Seurat的cluster标签。这样cluster_cells(cds)会被我们用clusters覆盖learn_graph()会以这个分区结构去学习图。第二cdsreduce_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 各方案对比表方案表达矩阵降维坐标聚类来源批次效应处理适用场景风险方案一RNA countsMonocle自己计算Monocle聚类不处理依赖细胞自然汇聚单一/同批次样本谱系清晰多批次时轨迹易被批次割裂方案二RNA countsSeuratHarmony UMAPSeurat clusterHarmony在坐标层面校正多批次/多样本合并分析需手动传坐标操作略繁琐方案三SCT countsMonocle自己计算或外部Monocle/Seurat均可SCT层面部分校正快速出图或所有下游已完成SCT轨迹结构可能失真需验证方案二最符合“轨迹分析用原始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_matx round(expr_matx)) cat(is integer counts:, is_integer, \n) # 查看范围 range(expr_matx)如果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_objreductions$umapcell.embeddings cdsreduce_dim_aux[[UMAP]] - list(model list(umap_coords umap_coords)) # 手动注入聚类信息 cdsclusters[[UMAP]] - seurat_objmeta.data$seurat_clusters names(cdsclusters[[UMAP]]) - rownames(seurat_objmeta.data) cds - learn_graph(cds)这里需要注意umap_coords的行名顺序。Monocle内部会通过细胞名匹配坐标因此矩阵行名必须和colnames(cds)完全一致。如果不一致learn_graph()会报错或者画出错乱图。建议传之前先做一次umap_coords - umap_coords[colnames(cds), ]强制对齐顺序。另外cdsclusters$UMAP这个名字不是随便起的。它对应reduce_dim_aux$UMAP意思是在UMAP这个降维结果上做的聚类。Monocle里cluster_cells()运行后也会生成类似结构。如果你跳过cluster_cells()直接手动塞标签记得cdsclusters$UMAP的名和cdsreduce_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)来固定根。拟时序结果里cdsprincipal_graph_aux[[UMAP]]$pseudotime存储了每个细胞的拟时序值。后续plot_genes_in_pseudotime()或graph_test()都会用到它。注意如果你在自己绘图时发现拟时序最大最小值只有零和一很可能是传递根细胞的方式出错了或者你传的数据源不对不要让这种结果进入下游分析。4.4 拟时序下游分析graph_test与差异基因拟时序算出来不等于结束。真正有生物学意义的结论来自“哪些基因随拟时序变化”。Monocle3提供了graph_test()它基于空间自相关统计量Morans I找到在轨迹上表达模式非随机的基因。gene_fits - graph_test(cds, neighbor_graph principal_graph, cores 4)拿到结果后重点看morans_test_stat和q_value过小的q值乘以Moran统计量后基因往往就是关键状态转换因子。这个步骤同样强依赖表达矩阵的可靠性。如果用SCT或integration矩阵跑Morans I会被严重高估或低估你的候选基因列表会失控。我见过有人跑出上千个显著基因结果一半是线粒体基因和核糖体蛋白基因这就是典型的输入矩阵有问题。5. 常见问题与排查技巧实录5.1 错误1直接传SCTcounts症状轨迹分支乱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_matx)如果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矩阵匹配然后看cdsclusters$UMAP是否和colnames(cds)一致最后检查根细胞是不是选错了换个根细胞或根节点重跑。很多情况下问题不是算法参数而是前面数据传递时埋下的雷。6. 最后想说的我在实际项目里来回试过五六种组合最终固定下来的干活模板就是方案二Seurat负责整合和注释Monocle3只拿RNA counts做轨迹构建外部降维坐标负责给轨迹锚定细胞空间。这套组合省了我大量调参时间也不用担心批次效应把轨迹撕得四分五裂。另外一个特别想提醒的细节是不要迷信“整合后数据更准”这句话。integration解决的是跨样本比较问题轨迹分析解决的是细胞状态过渡问题两者的数学目标并不一样。先想清楚你的生物学问题到底需要哪个工具去回答再决定把什么数据喂进什么算法。拟时序分析本质上是在基因表达梯度里寻找方向任何一步的过度校正都可能把这个梯度抹平。如果你现在正被Monocle的轨迹问题折磨第一件事去检查你传进去的表达矩阵到底是不是原始counts。不是的话直接停机改数据别在参数和图形美化上浪费时间。改完之后你会回来感谢这条建议的。