Monocle2的拟时分析做完之后很多人停在了“画出一条漂亮的轨迹图”这一步——细胞按颜色从浅到深排开伪时间轴铺在树状结构上确实好看。但一张轨迹图只能证明“分化确实发生了”它回答不了“沿着这个轨迹细胞的哪些生物学程序在改变”这个问题。要把轨迹和功能联系起来就必须走到基因模块层面先让数千个随拟时变化的基因按表达趋势聚成几组再对每一组做GO富集解析才能把“这个模块在干嘛”翻译成生物学结论。这篇文章把我自己跑这套流程的完整思路、代码细节和踩过的坑整理出来从Monocle2的拟时值提取、基因模块划分到GO富集和结果解读一条线讲清楚。适合正在做单细胞分化轨迹分析、对Monocle2已经基本入门但卡在“轨迹之后做什么”的同学参考。1. 为什么拟时分析要走到“基因模块”这一步1.1 单基因轨迹的三个常见瓶颈最早我做拟时分析时习惯把几个明星marker基因拉出来用plot_genes_in_pseudotime画三张趋势曲线比如干性基因下降、分化基因上升然后结论就写完了。但做多了就会发现这条路有三个瓶颈。第一单基因只能回答“这一个基因是否随时间变化”回答不了“这一群基因是否协同变化”。拟时轴上真正重要的生物学过程比如糖酵解重编程、线粒体代谢转换、细胞周期退出动辄涉及几百个基因协同表达。你挑三五个marker只是管中窥豹。第二成千上万个基因都在随拟时波动靠肉眼挑根本不现实。differentialGeneTest通常会返回一两千个显著基因一个人手工去看这些基因的趋势既低效又主观。第三也是最关键的单基因没有“通路”语义。哪怕你看到基因A在后期显著升高你也很难立刻说清楚这代表什么过程但如果是一组基因同时富集到“氧化磷酸化”结论就非常明确。这就是为什么需要基因模块。所谓基因模块我理解为一群表达趋势高度相似的基因在拟时早期高表达或者在拟时后期骤然升高或者中间出现一个峰值。把这些基因聚成4到8个模块之后每个模块内部基因的行为一致模块和模块之间又互相区分于是每一段拟时上发生的“程序切换”就变成了几个模块之间的接力赛。1.2 基因模块和GO富集的配合逻辑基因模块划分出来以后下一步自然就是功能注释。模块本身是一串基因ID如果不做富集它就是无名无姓的几十几百个符号。GO富集做的事情很简单拿模块里的基因和一个背景基因集比较用超几何检验看哪些基因本体论条目Gene OntologyGO条目在这批基因里出现得异常多比如“细胞周期”“DNA复制”“突触传递”等等。GO富集是给模块“贴标签”的过程。有了标签你才能在论文里写“module 1在早期富集到细胞周期相关条目提示细胞在拟时早期仍处于增殖状态module 4在后期富集到神经递质传递相关条目与功能成熟阶段吻合”。所以完整流程是一条链Monocle2拟时轨迹 → 拟时相关基因 → 基因模块 → GO富集 → 生物学解读。第2节先把Monocle2这一段快速过一遍因为拟时值的质量决定后面一切。2. 从表达矩阵到拟时值Monocle2全流程回顾2.1 输入对象构建一步错步步错的CellDataSetMonocle2的所有操作都围绕一个CellDataSet对象cds展开。构建这个对象需要三个输入表达矩阵、细胞注释信息phenoData、基因注释信息featureData。我在实际项目里最常踩的坑都在这一层后面专门讲这里先给一个可用的模板。library(monocle) # 表达矩阵基因在行细胞在列行为矩阵或Matrix都行 expr_matrix - as.matrix(GetAssayData(seurat_obj, slot counts)) # 细胞注释数据框行名必须是细胞barcode cell_metadata - data.frame( cell colnames(expr_matrix), cell_type seurat_obj$cell_type, sample seurat_obj$sample, row.names colnames(expr_matrix) ) # 基因注释数据框行名是基因名至少要有一列gene_short_name gene_annotation - data.frame( gene_short_name rownames(expr_matrix), row.names rownames(expr_matrix) ) pd - new(AnnotatedDataFrame, data cell_metadata) fd - new(AnnotatedDataFrame, data gene_annotation) cds - newCellDataSet( expr_matrix, phenoData pd, featureData fd, expressionFamily negbinomial.size() )expressionFamily参数值得多说一句。UMI计数数据10X、Drop-seq、InDrop这些都是通常用negbinomial.size()因为UMI本身是离散计数负二项分布符合它的特征。如果你的数据是TPM/RPKM这类连续值才考虑tobit()或gaussianff()。选错分布族会导致后面的estimateSizeFactors和estimateDispersions结果很奇怪轨迹也可能被噪声带跑。构建之后必须依次做两步estimateSizeFactors和estimateDispersions。不要跳过我看过有人直接对cds调reduceDimension结果轨迹完全随机。cds - estimateSizeFactors(cds) cds - estimateDispersions(cds)2.2 轨迹推断与拟时方向判定接下来是降维和排序cds - reduceDimension( cds, max_components 2, reduction_method DDRTree, norm_method log, pseudo_expr 1 ) cds - orderCells(cds)DDRTree是Monocle2默认推荐的降维方法它本质上做的是反向图嵌入在高维空间里找一条树形主曲线然后把每个细胞投射到这条曲线上投射位置转换成拟时值。这里的max_components2并不代表只保留2个基因它是指主曲线的低维表示是二维。绝大多数情况下保持默认2就够设成3虽然偶尔能让分支更清楚但后续可视化、根状态指定都会麻烦。orderCells之后pData(cds)$Pseudotime里就是每个细胞的拟时值pData(cds)$State是每个细胞落在哪个分支状态。到这里Monocle2常规流程就算跑完了。但我想特别强调拟时值本身是有方向的而这个方向不一定符合你的生物学预期。orderCells默认会把细胞数最多的那个状态当作根状态把根状态附近的细胞设为拟时0。如果你的研究对象是“从祖细胞向终末细胞分化”祖细胞数量反而少默认根就可能落在终末分化细胞那头导致整个拟时轴方向是反的。怎么检查方向最直观的办法是把已知的祖细胞marker和终末marker分别投影到轨迹上或者直接看plot_cell_trajectory(cds, color_by Pseudotime)时颜色渐变是否从你预期的起点出发。2.3 拟时方向纠正的技巧如果发现方向反了不用重新跑Monocle2提供了显式指定根状态的接口。你先画轨迹图看每个State的位置找到应该作为起点的那个状态编号然后# 假设确定根状态是State 3 cds - orderCells(cds, root_state 3)一句话就能把拟时方向翻过来。之后的基因模块、GO富集全部基于方向纠正后的拟时值这一步千万别省。我在一个造血分化的项目里就吃过这个亏没纠正方向k-means分出来的模块趋势完全颠倒早期模块富集到的全是终末功能条目差点得出一个反的结论。另外还要提醒一句如果你用到differentialGeneTest建议用纠正方向后的cds重新跑一次。虽然理论上~sm.ns(Pseudotime)只关心基因表达和拟时之间的非线性关系方向反转对显著基因列表影响不大但对拟合曲线的形状有影响进而影响下游聚类的边界。3. 两种基因模块划分路线实操3.1 路线一BEAM找分支相关基因再用热图切模块如果你的轨迹有明确的分支结构比如从共同前体分化为两种细胞命运那BEAM就是最常用的起点。BEAM的全称是Branched Expression Analysis Modeling它检验的是基因表达是否在某个分支点前后出现显著性差异也就是寻找“命运决定基因”。# branch_point数值从轨迹图上看一般取1有多个分支时逐个尝试 BEAM_res - BEAM(cds, branch_point 1, cores 4) BEAM_res - BEAM_res[order(BEAM_res$qval), ] sig_genes - row.names(subset(BEAM_res, qval 0.01)) cat(BEAM显著基因数量, length(sig_genes), \n) cds_sig - cds[sig_genes, ] plot_pseudotime_heatmap( cds_sig, num_clusters 6, return_heatmap TRUE )plot_pseudotime_heatmap是Monocle2自带的函数它做的事情是把显著基因的表达量按细胞拟时排序后做成热图并在此基础上对基因做k-means聚类。num_clusters就是你想分的模块数。热图每一列是一个细胞从左到右拟时递增每一行是一个基因基因按聚类结果排序并用颜色条标出归属模块。这个函数的优点是快、参数少适合跑通流程、先看个大概。但缺点也很明显它内部对基因做的是z-score标准化热图上只能看出每组基因相对自身均值的波动模式看不出模块间的绝对表达量差异而且num_clusters选了以后具体每个聚类包含哪些基因你需要从热图对象里手动把kmeans结果取出来。3.2 路线二自取表达矩阵做k-means自由度更高如果轨迹没有明显分支或者你想完全掌控聚类细节我建议自己提取拟时排序后的表达矩阵来做k-means。这里给出我常用的代码模板# 先准备拟时相关基因用differentialGeneTest找随拟时显著变化的基因 diff_test - differentialGeneTest( cds, fullModelFormulaStr ~sm.ns(Pseudotime) ) sig_genes - row.names(subset(diff_test, qval 0.01)) # 提取表达矩阵并按拟时排序 plot_df - exprs(cds)[sig_genes, ] pseudotime_vec - pData(cds)$Pseudotime plot_df - plot_df[, order(pseudotime_vec)] # 每行做z-score标准化让不同表达量水平的基因可比较 scaled_mat - t(scale(t(plot_df))) # k-means聚类k是模块数 set.seed(42) km - kmeans(scaled_mat, centers 6, nstart 25, iter.max 50) # 按聚类编号整理基因 module_list - split(rownames(scaled_mat), km$cluster)这里的nstart 25比较重要。k-means对初始中心敏感nstart是随机初始中心重复次数最后取组内平方和最小的结果。如果设成1结果可能每次跑都不一样。设种子set.seed是为了结果可复现这点投稿时尤其重要审稿人如果让你补实验记录你至少要能跑出同一个模块划分。另外我强烈建议在聚类前把表达量过滤一遍。differentialGeneTest返回的显著基因可能有两三千个但里面混着大量低表达、低方差的“垃圾基因”。k-means对这类基因很敏感它们随机抖动会把聚类边界带偏。建议至少用rowMeans和rowSds做一轮过滤keep - rowSds(scaled_mat) 0.5 scaled_mat - scaled_mat[keep, ]3.3 模块数量怎么定以及两条路线的取舍num_clusters或centers到底选几没有标准答案。我个人的经验是4到8之间先试最终结合富集结果定。4个模块通常是“早期高、早期低、后期高、后期低”6个模块会多出“中间峰值型”“后期骤升型”等更细的模式超过8个模块每个模块的基因数变少GO富集检验效力下降容易富出一堆不显著的结果。判断标准我认为有两个一是热图里每个模块的轮廓是否清晰模块内部基因的趋势是否一致二是每个模块做完GO之后能不能解读出明确的生物学故事。如果两个相邻模块富集到几乎一模一样的GO条目说明模块切得过细了合并它们更合理如果某个模块富集不出任何条目先别急着加阈值看看是不是模块里基因数太少可能要把相邻模块合并。路线一BEAM热图适合有分支轨迹、关注细胞命运选择的场景路线二差异基因k-means适合单一连续分化过程、你想把整个拟时轴划成几个连续阶段的场景。实际项目中我经常两条路线都跑先用路线一快速确认分支是否与已知细胞类型对应再用路线二把主干的基因模块彻底拆清楚。4. GO富集从基因ID转换到结果过滤的完整操作4.1 基因ID转换是常见的翻车点模块拿到手最兴奋的时候最容易翻车的环节是基因ID转换。clusterProfiler的enrichGO要求基因ID是ENTREZID至少也要求是能识别的格式。但你的模块基因大概率是SYMBOL因为Seurat和Monocle的featureData经常用gene_short_name所以必须转一次。library(clusterProfiler) library(org.Hs.eg.db) # 以人类的模块为例模块1 module_gene_symbol - module_list[[1]] gene_entrez - bitr( module_gene_symbol, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db )这里最常见的坑是bitr返回的基因数比输入少。因为有些SYMBOL在数据库中匹配不到ENTREZID比如线粒体tRNA基因、新注释的lncRNA、版本更新后改了名称的基因。如果丢失比例很高整个模块的富集结果就会失真。所以bitr之后我会加一行输出cat(输入, length(module_gene_symbol), 个基因成功转换, nrow(gene_entrez), 个\n)如果转换率低于80%我一般会检查是不是OrgDb选错了物种。人的数据必须用org.Hs.eg.db小鼠用org.Mm.eg.db斑马鱼用org.Dr.eg.db这个不多说。还要注意基因名格式有的上游流程给你的是Ensembl ID或NCBI RefSeq不能直接当SYMBOL传进去fromType要跟着改成ENSEMBL或REFSEQ。4.2 enrichGO参数设置转换完之后正式跑富集ego - enrichGO( gene gene_entrez$ENTREZID, universe background_entrez, # 背景基因下面解释 OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, # BP / CC / MF 三选一或ALL pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.2, readable TRUE # 结果里把ENTREZID换回SYMBOL )readable TRUE这个参数强烈建议打开。不然你拿到富集结果geneID列全是ENTREZID数字还得再转一次才能看懂基因是谁很多新手在这个细节上浪费时间。pvalueCutoff和qvalueCutoff是两道关卡。pvalueCutoff过滤原始p值qvalueCutoff过滤BH校正后的q值。实际项目中我倾向于设定pvalueCutoff 0.05qvalueCutoff 0.2因为在几十个模块里同时做富集时q值会非常保守卡到0.05可能很多模块什么都富不出来。如果你只追求结果稳健可以都设0.05但要做好很多模块没词条的心理准备。4.3 背景基因最容易被忽略却影响最大的参数universe这个参数很多人不填或者随便填。不填的话clusterProfiler默认使用你的OrgDb里的全部基因作为背景。问题在于你的单细胞转录组只检测到了一万多个基因而人类基因组注释有接近两万个蛋白编码基因直接用全基因组做背景会显著低估富集显著性。可能某个模块里的基因在全基因组背景里占比不大但放在你实际检测到的基因里已经是高度富集了。我的建议是把Monocle2分析之前过滤后剩下的所有表达基因作为背景也就是rownames(cds)对应的全部基因。用SYMBOL转一次ENTREZID存成背景向量all_symbol - rownames(cds) all_entrez - bitr(all_symbol, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) background_entrez - unique(all_entrez$ENTREZID) ego - enrichGO( gene gene_entrez$ENTREZID, universe background_entrez, ... )如果你的研究思路是先筛出“所有随拟时显著变化的基因”再在这个范围内看哪些GO富集那也可以把背景设为这些显著基因全体。这种做法的逻辑是既然没有随拟时变化的基因不在讨论范围内那么背景就应该是所有变化基因。两种策略都可以但必须在方法部分写清楚。我自己的习惯是用全部表达基因做背景这样结论更普适。5. 富集结果怎么看为什么有些模块富不出来5.1 BP、CC、MF各看什么GO富集结果有BP生物过程、CC细胞组分、MF分子功能三个本体很多人一上来就选ont ALL然后结果表里三个混在一起看着乱也不好解读。我的经验是追踪分化或发育过程的动态变化优先看BP和CC。BP回答的是细胞在“做什么”比如“细胞周期G1/S转换”“糖酵解过程”“轴突引导”CC回答的是这些基因产物“在哪里工作”比如“线粒体内膜”“突触后膜”“核小体”。两个本体的优先级视项目而定。MF相对抽象一般是“结合”“催化活性”这类分子层面的描述对于“细胞从增殖到分化”这种宏观问题帮助不大。如果你是想找转录因子或者激酶这类有明确分子功能的模块再单独看MF。5.2 怎么从结果表里判断富集质量enrichGO返回的是一个enrichResult对象可以直接as.data.frame()转成表格。我拿到表后不会只看p值而是看三个数第一个是GeneRatio表示模块里落在该GO条目上的基因数占模块总基因数的比例。第二个是BgRatio表示背景基因落在该GO条目上的比例。GeneRatio / BgRatio就是富集因子enrichment factor大致代表这个条目在模块里的富集强度。第三个是Count即模块里实际注释到这个条目的基因数。这三个数要合起来看。一个GO条目即使p值极其显著如果Count只有3个基因那这个结论脆得像纸换个过滤参数可能就没了。我一般要求Count 5才保留。GeneRatio太小的条目比如0.02哪怕q值小于0.05说服力也不强。5.3 一个典型的解读路径我举一个例子。假设你的拟时轨迹是一条从干细胞到神经元的连续分化曲线k-means分了6个模块module 1在拟时早期高表达GO富集到“有丝分裂细胞周期”“DNA复制”“姐妹染色单体分离”这个结论在生物学上非常合理干细胞阶段增殖活跃。module 3在中期出现峰值富集到“神经嵴细胞迁移”“轴突导向”“细胞形态发生”提示细胞正在经历形态和位置变化。module 6在后期高表达富集到“化学突触传递”“离子跨膜运输”“神经元成熟”说明终末阶段功能基因全面开启。这样一条解读路径下来模块1→3→6就是“增殖→迁移改造→功能成熟”的完整故事。如果你的结果出现module 1富集到突触传递、module 6却富集到细胞周期赶紧回头检查拟时方向大概率是方向设反了。5.4 多模块富集结果的清晰可视化每个模块单独跑一次enrichGO然后挑每个模块top 5到top 10的GO条目拼在一起画点图。我用的比较多的是clusterProfiler自带的dotplot以及把多个模块结果合并后自定义ggplot2画图。# 单个模块的点图 dotplot(ego, showCategory 10)如果想在一张图里横向比较多个模块可以手动给每个模块的富集结果加一列module标签合并后按模块和p.adjust排序再用ggplot2画气泡图。气泡大小是Count颜色是p.adjust。这种图在论文里很常见信息密度大审稿人也容易看明白。6. 实操中遇到的坑和补救建议6.1 一个模块富集出几百个GO怎么收敛模块基因数量多时比如一个模块500个基因跑BP富集可能哗啦一下出来四五百条显著GO。这时候别急着全写上论文先用simplify去冗余。GO条目之间不是独立的很多条目语义高度重叠比如“细胞周期”和“细胞周期过程”这两个词基因集基本一样富集结果里会同时出现。simplify用语义相似度把冗余条目合并成一个代表性条目ego_simple - simplify(ego, cutoff 0.7, measure Wang)cutoff 0.7是相似度阈值越大合并越少越小合并越狠。我用0.7比较多。去掉冗余之后结果从四五百条能收敛到几十条接下来选top条目就轻松多了。6.2 拟时相关基因太多时的过滤策略differentialGeneTest的q值经常给到0.001以下一筛就有两三千个基因。如果把所有显著基因都丢进k-means模块会大而杂每类趋势里都混着一些和主流趋势无关的基因。我的做法是在显著基因基础上再加一道“表达量”过滤。具体来说要求基因至少在10%的细胞里检测到且平均表达量normalized之后大于某个阈值。过滤之后基因数量通常会从三千降到一千五左右k-means聚出来的模块轮廓明显更干净。6.3 多分支轨迹要小心“主干模块”和“分支模块”混在一起如果轨迹有两个以上的分支比如一个共同前体分别分出神经元和胶质细胞两条命运那么BEAM找到的分支基因里面有一部分是只在神经元分支特异高表达的另一部分只在胶质分支特异。把它们混在一起做k-means可能会聚出一个“两种命运基因的混合物模块”这个模块的GO富集会非常混乱既有神经元发育又有胶质细胞分化。遇到这种情况我建议先按细胞状态State把分支拆开。只保留一个分支的细胞重新提取表达矩阵并排序再分别做基因模块和GO富集。Monocle2本身支持用subset按State过滤细胞branch_cells - row.names(pData(cds))[pData(cds)$State %in% c(2, 3)] cds_branch - cds[, branch_cells]然后对这个分支子集重新跑聚类。这样虽然工作量翻倍但得到的结论清晰很多每个分支都能讲出一个独立的故事。6.4 关于Monocle2的版本问题最后说一个每次写教程都得提的事。Monocle2这个包从Bioconductor 3.14之后就不再随新版本发布如果你在最新版R里直接BiocManager::install(monocle)大概率会报错。我目前常用的安装方式是从GitHub装install.packages(devtools) devtools::install_github(cole-trapnell-lab/monocle-releasemonocle2)装完之后加载包名仍然是library(monocle)不影响使用。如果你的集群上R版本太老也可以去Bioconductor的归档版本里找对应R版本的旧二进制包。另外Monocle2和Monocle3不要混着装两者的函数名和对象结构有冲突我见过同事同时加载两个gem导致newCellDataSet都调不对的情况。这套流程跑下来你手里的产物不只是一张轨迹图而是一组“模块GO条目”的对应关系。写论文时表里可以放每个模块的top GO讨论里可以按“早期模块→中期模块→后期模块”的顺序把分化过程的程序切换讲清楚。我个人在最近的几个项目里都靠这套流程在审稿时少挨了不少关于“你的标记基因选择为什么主观”的质疑——因为基因模块和GO富集至少在聚类和统计意义上是有可复现依据的。希望这篇整理的实操笔记能帮你少走几步弯路。