搞单细胞转录组的人绕不开拟时轨迹分析而只要做到拟时轨迹就迟早会面对Monocle2这个老牌工具。很多刚上手的朋友跑完orderCells出了一条漂亮的轨迹就以为万事大吉但下一步往往卡住那些沿伪时间变化的基因模块到底对应了什么生物学功能怎么才能从一堆候选基因里提炼出细胞命运决定的线索这篇文章就来把Monocle2拟时轨迹基因模块的GO富集解析串成一条可复现的完整流程从软件原理、差异基因筛选、模块聚类到GO富集实操把每一步背后的逻辑和容易踩的坑都讲清楚。整个过程适合刚跑通Monocle2但不知道怎么继续深入的人也适合那些已经拿到基因列表、却不知道如何往下游功能分析的人。我会按自己的实际操作习惯来写代码直接给参数直接讲希望能帮你少走一些弯路。1. 拟时轨迹分析的核心思路与Monocle2选型逻辑1.1 Monocle2如何把细胞“捋”成一条线反向图嵌入原理Monocle2的核心逻辑不复杂但你必须先明白它算出来的不是“时间”而是细胞在表达状态空间里的距离。我们可以把每个细胞的数以万计基因表达量看成高维空间里的一个坐标点拟时分析要做的就是在这个高维点云里找出一条能串联起不同状态的路径。Monocle2采用反向图嵌入Reverse Graph EmbeddingDDRTree作为降维工具先在高维空间里学习一个树状结构的低维流形再把所有细胞按表达相似度投影到这棵树上。这棵树的根到叶子的距离就是细胞之间的“拟时距离”本质上是状态转换的进度而不是真实的昼夜时间。所以你会发现在同一个时间点取样不同细胞也可以被放在伪时间轴的不同位置因为它们在分化过程中所处阶段不一样。我最初刚接触Monocle2时总把轨迹上的顺序想成严格的先后序列这其实是个认识误区。DDRTree生成的路径可以是分叉的一个母节点可以分出两个子节点这就对应了细胞命运决定的分叉点。Monocle2的拟时轨迹能识别这种分叉所以在分析免疫细胞分化、肿瘤细胞干性演进、发育过程这类有“分叉”故事的课题里特别适用。1.2 为什么选Monocle2而不是Monocle3或Slingshot单说拟时分析市面上现在有不少新工具。Monocle3已经发布多年Slingshot、Scanpy里的diffmap也都有各自的优势。但Monocle2目前仍然有不可替代的使用场景尤其是做基因模块下游富集解析时它的差异基因检验框架和伪时间热图配合得最好。我用一张表把这几个工具的使用特点大致说清楚工具原理优点适合场景Monocle2DDRTree反向图嵌入轨迹连续稳定BEAM分支检验成熟热图可视化方便中小规模样本数万细胞以内需要提取伪时间相关基因模块Monocle3UMAP学习图大规模数据支持好分区模块更细超大样本量需要精细分区轨迹Slingshot主曲线拟合方法透明速度较快已有明确起点和聚类标签如果要追问为什么现在很多文献还在用Monocle2一个很现实的原因是它的伪时间基因模块提取流程非常顺。Monocle3虽然也能排序但下游做BEAM和伪时间热图时需要额外改造很多朋友的踩坑成本反而更高。Slingshot则强依赖你预设的起点和聚类注释对没有先验知识的流程不太友好。Monocle2的局限我也直说它对细胞数量有瓶颈超过十万个细胞后DDRTree降维会比较吃力而且它强制所有细胞被安排在一棵树上如果数据里存在相互独立的多条谱系分化Monocle2的处理效果会不那么理想。所以选型前需要先看清自己的数据形态。1.3 从表达矩阵到CellDataSet的数据准备Monocle2一切分析都建立在CellDataSet对象上。你手里常见的表达矩阵通常是基因×细胞的数值矩阵首先要把行名设为基因ID可以是symbol也可以是Ensembl ID列名设为细胞条形码。接下来需要给细胞和基因分别准备元数据。用代码来演示一下从矩阵构建CellDataSet的基本流程library(monocle) expr_matrix - readRDS(expr_matrix.rds) # 基因 x 细胞 cell_metadata - readRDS(cell_metadata.rds) # 细胞注释 gene_annotation - data.frame(gene_short_name rownames(expr_matrix)) rownames(gene_annotation) - rownames(expr_matrix) pd - new(AnnotatedDataFrame, data cell_metadata) fd - new(AnnotatedDataFrame, data gene_annotation) cds - newCellDataSet( expr_matrix, phenoData pd, featureData fd, lowerDetectionLimit 0.5, expressionFamily negbinomial.size() )其中expressionFamily参数的选择很重要。如果你的表达量是UMI count或者read count推荐用negbinomial.size()如果是经过大小归一化后的TPM/FPKM这类连续值则需要用gaussianff()。我见过不少朋友在构建时直接用默认的negbinomial导致后面estimateSizeFactors和estimateDispersions报错或者结果异常这一步务必检查清楚。数据准备阶段的核心目的是让后续分析输入干净把低表达的基因留下太多会增加降维噪音把细胞过滤得不干净轨迹会出假分支。按我的经验至少要先过滤掉在极少数细胞中表达的基因。1.4 必须做的事过滤、归一化与离散度估计在我实际跑过的流程里正式降维之前有四步是雷打不动的detectGenes、estimateSizeFactors、estimateDispersions以及选高变基因。cds - detectGenes(cds, min_expr 0.1) expressed_genes - row.names(subset(fData(cds), num_cells_expressed 10)) cds - cds[expressed_genes, ] cds - estimateSizeFactors(cds) cds - estimateDispersions(cds)这里num_cells_expressed 10是我常用的经验值但具体的阈值要根据你的细胞总数和测序深度来调整。细胞数越多阈值可以适当提高测序深度越低阈值要放宽避免丢失真实表达模式。之后的降维有些教程会先选高变基因有些则直接对全部过滤后的基因跑reduceDimension。我的建议是先做一步差异表达基因筛选再由这些基因参与DDRTree降维。这样既能保留有用的生物学信号又能让运算速度加快很多尤其是细胞数量在几万级别时全基因跑DDRTree是真的会等到怀疑人生。2. 基因模块提取从拟时轨迹到候选基因2.1 判断轨迹形态有分支原理没分支走另一套逻辑拿到一条轨迹后第一件事不是立刻提取基因而是看轨迹有没有分叉点。你可以用plot_cell_trajectory(cds, color_by Pseudotime)观察细胞在树上的分布也可以用plot_complex_cell_trajectory查看分支结构。Monocle2的分支节点在内部叫Branch Point是轨迹中某个细胞状态可以朝向两个不同方向发展的位置。如果你的轨迹只有一个分支比如细胞从A状态分化为B和C两种终末状态就可以用BEAMBranch Expression Analysis Modeling来检验哪些基因沿着从分支点到终点的表达模式有显著变化。这类基因往往正好对应了细胞命运决定的关键调控子。但如果你的数据没有明显分支轨迹基本是一根连续的线那么就不要强行用BEAM改用differentialGeneTest以伪时间作为协变量来筛选随时间显著变化的基因。2.2 BEAM与differentialGeneTest的参数细节先说说有分支的情况BEAM的代码非常简单但参数值得仔细讲一下cds - orderCells(cds) BEAM_res - BEAM(cds, branch_point 1, cores 8) BEAM_res - BEAM_res[order(BEAM_res$qval), ] BEAM_genes - subset(BEAM_res, qval 1e-4)其中branch_point 1需要你确认在使用orderCells后细胞被正确分配到树枝上。如果树枝顺序不对BEAM跑出来的基因分会乱掉导致后续模块完全失去生物学意义。我习惯在跑BEAM前先用plot_cell_trajectory(cds, color_by State)检查State分配和分支结构是否合理。cores参数影响计算速度如果BEAM运行时间过长多半是这里没设置好。加上cores 8之后实际运行速度会有质的提升。对于没有分支的轨迹用differentialGeneTestdiff_test - differentialGeneTest(cds, fullModelFormulaStr ~sm.ns(Pseudotime)) diff_test - diff_test[order(diff_test$qval), ] diff_genes - subset(diff_test, qval 0.01)~sm.ns(Pseudotime)是Monocle2里默认使用的自然样条基函数用来拟合伪时间的非线性表达趋势。如果不加sm.ns模型拟合的是线性关系很多表达先升后降或者波动型基因就会漏掉。筛选阈值项目中常见用的是qval 0.05但如果你后续还要做模块聚类我会建议适当收紧到qval 1e-4甚至更严格。因为模块聚类对噪音基因很敏感宽松阈值会把大量不相关基因拉进模块最后GO富集结果非常分散难以解释。2.3 把筛选出的基因按表达趋势聚成模块拿到差异基因列表后Querier通常会用plot_pseudotime_heatmap直接可视化这个函数内部会按伪时间把基因表达拟合成曲线然后用k-means聚类。代码长这样plot_pseudotime_heatmap( cds[candidate_genes, ], num_clusters 4, cores 8, show_rownames TRUE )这个函数输出的热图非常直观每一行是一个基因列是按伪时间排序的细胞行方向按照聚类模块分组能立刻看出来哪些基因在分化早期高表达、哪些基因在终末阶段才被激活。不过我建议不要只用这一个函数如果直接在R里跑plot_pseudotime_heatmap返回的聚类结果比较“一次性”不方便直接拿到每个模块的基因列表做后续GO富集。我常用的方法是手动提取表达趋势矩阵并自行聚类library(tidyverse) pseudo_expr - data.frame(t(exprs(cds[candidate_genes, ]))) pseudo_expr$Pseudotime - cds$Pseudotime tet - pseudo_expr %% pivot_longer(-Pseudotime, names_to gene, values_to expr) # 这里更推荐先用smooth.spline拟合每个基因的表达趋势 # 再把拟合值按伪时间排序做scale归一化 gene_matrix - vapply( candidate_genes, function(g) { fit - smooth.spline(pseudo_expr$Pseudotime, pseudo_expr[[g]]) predict(fit, x order(pseudo_expr$Pseudotime))$y }, numeric(nrow(pseudo_expr)) ) gene_matrix_scaled - t(scale(t(gene_matrix))) km - kmeans(gene_matrix_scaled, centers 4, nstart 20) module_list - split(names(km$cluster), km$cluster)这么做的好处是你能把每个模块的基因列表直接存成文本文件方便批量导入clusterProfiler之类的工具。另外先做趋势拟合再做聚类可以有效去除单细胞表达矩阵里的采样噪声如果直接对原始表达量做聚类会发现同一模块里的基因表达波形很碎模块划分极不稳定。2.4 模块数量怎么选宁可少不要贪多关于num_clusters的取值很多新手会在4还是在6之间反复横跳。我自己的经验是模块数量不宜过多每个模块里的基因数至少要保持在20到50个以上。过少基因的模块做GO富集时统计功效会严重不足经常会得到一堆语义宽泛但毫无区分度的通路对解释细胞状态转变帮助不大。一开始可以把模块数设成3到6然后观察每个模块基因的平均表达曲线是否具备不同的时序特征。如果两个模块的表达趋势几乎重合说明聚类数设大了如果一个模块里既有早高表达基因又有晚高表达基因说明聚类数偏小。热图配合平均曲线检查比任何统计指标都直观。模块本身也不是越少越好。基因太少则富集无意义基因太多则模块内功能混杂富集出来一堆本体论大词难以提炼核心结论。在实际项目里我通常会在num_clusters 5的基础上根据富集结果微调最终选定4到6个模块这样既有清晰生物学故事下游分析也足够稳。3. GO富集分析把模块变成功能注释3.1 选对工具与基因ID模块基因列表准备好之后就到了整个流程的收尾环节GO富集。这一环节的目标是搞清楚每个模块里的基因到底参与了哪些生物过程、分子功能和细胞组分。工具层面我强烈推荐R包clusterProfiler它支持批量输入多个模块也能直接可视化省去切换网页工具的麻烦。但要用它首先得把基因ID统一成注释包支持的格式。如果是人类数据需要org.Hs.eg.db小鼠数据用org.Mm.eg.db斑马鱼、拟南芥等也有对应的OrgDb包只是支持范围不一定全。如果你的物种注释偏冷门也可以考虑用DAVID和g:Profiler做在线富集但那样不方便对多个模块做统一管理和批量可视化。基因ID转换有一个细节clusterProfiler用了不同的keyType而很多单细胞注释文件默认给的是gene symbol。你要先用bitr转成ENTREZID再执行enrichGO。千万别把gene symbol直接丢进enrichGO的gene参数然后指望着readable TRUE帮你解决一切那样通常会报错或结果为空。3.2 一个可以直接跑的clusterProfiler流程下面这段代码是我在项目里反复使用的标准流程适用于人类数据其它物种只需替换OrgDb和相应的keyTypelibrary(clusterProfiler) library(org.Hs.eg.db) module_genes - readLines(module_1_genes.txt) gene_entrez - bitr( module_genes, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db ) # 如果担心背景基因问题可以计算所有检测到基因的Entrez列表 all_genes - readLines(all_detected_genes.txt) all_entrez - bitr( all_genes, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db ) ego - enrichGO( gene gene_entrez$ENTREZID, universe all_entrez$ENTREZID, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.2, readable TRUE ) head(egoresult)这里universe参数值得单独说明。GO富集的本质是超几何检验它需要知道“你从多少个基因里挑出当前这些候选基因”。如果universe不指定clusterProfiler默认用注释包里所有基因这会让你的显著性偏向于那些在基因组里本身就和某通路紧密相关的基因导致很多模块富集到相同的“基础代谢”或“拼接”通路丢失特异性。在实际分析时我认为更合理的universe是在你的单细胞数据里经过detectGenes过滤后仍然被检测到的所有基因。这样富集比较的基线才是“你能检测到什么”而不是“基因组上有什么”结果才更贴近真实生物学。3.3 多模块批量富集与可视化既然做了多个基因模块逐个跑enrichGO效率太低了。clusterProfiler提供了一个非常方便的函数compareCluster直接把模块列表打包处理module_list - list( module1 readLines(module_1_genes.txt), module2 readLines(module_2_genes.txt), module3 readLines(module_3_genes.txt) ) ck - compareCluster( module_list, fun enrichGO, OrgDb org.Hs.eg.db, keyType SYMBOL, ont BP, pvalueCutoff 0.05, qvalueCutoff 0.2, readable TRUE ) dotplot(ck, showCategory 5)compareCluster会自动把每个模块当作一个输入组返回一个合并的富集结果再用dotplot绘制点图每个模块一行点的大小表示富集到的基因数目颜色深浅表示显著性。这么做的好处是可以快速比较不同模块在功能上的差异比如模块1富集到细胞周期模块2富集到免疫应答模块3富集到代谢那整个细胞状态演进的逻辑骨架就出来了。如果你在dotplot输出里发现所有模块都富集到同一批通路不要急着往下写结论。先检查模块划分是否有问题再检查universe是否设置合理最后检查是不是差异基因列表有大量重叠。这三点经常是结果看起来“千篇一律”的根源。3.4 富集结果表格怎么看GO富集的结果表格里我平时重点看四列Description、GeneRatio、p.adjust和geneID。Description是功能通路名字比如“regulation of cell differentiation”“response to cytokine”GeneRatio表示你的模块基因里有多少个落在这个通路分母是模块总基因数分子是富集到该通路的基因数p.adjust是多重检验校正后的显著性一般小于0.05算有意义geneID列是具体的基因列表末尾如果显示/分隔的多个基因可以用来核对该通路是否真的和你想讲的故事有关。很多新手只看p值不看GeneRatio容易被一个只富集到2个基因的通路吸引住然后写出一大段解释。我的习惯是如果富集到某个通路的基因数少于5即使p值再小也只作为辅助参考不放到核心结论里。少于5个基因的富集往往经不起抽样的随机性考验换个样品或换个聚类参数可能就消失了。4. 结果解读、避坑与效率优化4.1 从GO术语反推生物学意义GO富集结果的解读本质上是要把轨迹上的基因模块和细胞命运决定事件联系起来。拟时轨迹和模块分析本身只负责“排序”和“找基因”真正告诉你怎么回事的是功能富集。举个例子假设你的轨迹是从naive T细胞到Th1/Th2分化分支点处BEAM筛出的差异基因会被聚成两个模块分支点前模块主要富集到T细胞受体信号通路、细胞激活相关过程分支点后某个模块富集到细胞因子介导的信号通路等。这时你就不能说“这些基因只在分化后期表达”而要结合模块表达曲线和分支顺序去讲哪一批基因先启动作出决定哪一批基因在终末阶段执行功能。我一直提醒自己和身边同事GO富集得到的生物学过程只揭示了“这一类基因倾向参与什么”而不等于“细胞一定发生了什么”。比如富集到“cell cycle”不代表细胞真的在增殖可能只是模块里包含了若干周期相关基因实际表达趋势未必和细胞周期状态同步。这个分寸感很重要尤其是写文章时千万别把相关性说成因果性。4.2 实操中最容易翻车的5个细节第一branch_point参数没确认。BEAM输出的p值是针对某个指定分支的如果轨迹有多个分支而你把branch_point设成1结果可能只反映了第一个分支的变化。最好先用plot_complex_cell_trajectory看节点排列再逐一跑BEAM。第二聚类前没有做标准化。不同基因的表达量本身就差异巨大高表达基因会主导欧氏距离导致聚类结果里高表达基因自成一团低表达但趋势重要的基因被混在一起。用scale按行标准化再聚类几乎是必须的。第三GO分析的universe没设置。如果没有指定背景显著性会被所有注释基因稀释模块差异性的信号会被冲淡。一定要把universe设为你实际检测到的基因集。第四基因ID转换不完整。bitr有时会返回NAs这是因为有新旧symbol不一致、线粒体基因、核糖体RNA等情况。建议转换后过滤掉NA再富集不要强行把NA丢进去跑。第五富集结果为空就盲目放宽阈值。我理解没有结果时想放宽一点的心理但盲目把qvalueCutoff调到0.9只会得到一堆无意义的宽泛通路。更理性的做法是检查模块基因数量是不是太少、universe是否合理、ont是否选得合适实在不行就换用MSigDB或者Reactome富集作为补充。4.3 常见报错速查表整理一个我在答疑中经常碰到的报错速查表方便你排查报错信息可能原因解决办法Error in BEAM: cds must have a branch轨迹没有分叉改用differentialGeneTestError in estimateDispersionsexpressionFamily使用错误UMI/count用negbinomial.size()TPM/FPKM用gaussianff()Error in bitr: some genes cannot be mapped基因名规范或物种匹配问题检查基因标识类型调整fromType过滤映射不上的基因No enrichment found in any gene set模块基因过少、universe过大或者阈值太严增加模块基因、缩小universe、合理放宽pvalueCutoffKilled或内存不足BEAM计算量过大用cores并行减少参与基因数或细胞数还有一个容易被忽略的地方单细胞矩阵里的基因名如果是Ensembl ID而你的差异基因列表用的是symbol两者没对齐时GO富集大概率会全军覆没。处理这种问题我通常在样品注释文件生成阶段就把基因名统一成同一种格式宁可在上游多花十分钟也不要在下游排查半天。4.4 效率优化怎样让流程跑得更顺手Monocle2最大的槽点是跑得慢尤其是BEAM。如果你的细胞数量在几万以上基因数量又在几千BEAM可能跑数小时。为了兼顾效率和结果稳定性我一般会采用“先粗筛再精算”的策略先用differentialGeneTest以比较宽松的阈值筛出一批候选基因再做BEAM因为BEAM的模型参数估计远比普通差异检验复杂基因数量少一半时间能省一大半。GO富集阶段clusterProfiler本身速度不算慢但对多个模块跑多次时尽量用compareCluster合并执行这样只用加载一次注释包内存占用和耗时都更友好。另外如果你的OrgDb包版本很旧会导致部分基因ID映射失败建议定期更新org.*.db包别让环境里的旧注释拖后腿。写在最后的一点心里话这整套流程我自己跑过很多遍最大的体会是Monocle2拟时轨迹只是给你一个细胞状态变化的抽象骨架真正让这个骨架长出肉、让生物学故事立住的是后续的基因模块和GO富集解析。很多朋友把70%的时间花在跑轨迹上留给功能分析的时间很少这其实有点本末倒置。如果让我回头重做一次我会先把数据质控和聚类前期工作做得更扎实再跑拟时拿到模块列表后优先用GO富集和文献相互印证而不是一味追求更小的p值。最后再分享一个小习惯把每个模块的基因列表单独存成文本文件旁边附带一列模块编号和一列简单注释。这样不管是后续换用Reactome、MSigDB还是做体外实验验证都能随时取用。看似不起眼但在项目周期拉长到几个月之后你会感谢当初这个微不足道的整理动作。
Monocle2拟时轨迹到GO富集：基因模块解析完整指南
搞单细胞转录组的人绕不开拟时轨迹分析而只要做到拟时轨迹就迟早会面对Monocle2这个老牌工具。很多刚上手的朋友跑完orderCells出了一条漂亮的轨迹就以为万事大吉但下一步往往卡住那些沿伪时间变化的基因模块到底对应了什么生物学功能怎么才能从一堆候选基因里提炼出细胞命运决定的线索这篇文章就来把Monocle2拟时轨迹基因模块的GO富集解析串成一条可复现的完整流程从软件原理、差异基因筛选、模块聚类到GO富集实操把每一步背后的逻辑和容易踩的坑都讲清楚。整个过程适合刚跑通Monocle2但不知道怎么继续深入的人也适合那些已经拿到基因列表、却不知道如何往下游功能分析的人。我会按自己的实际操作习惯来写代码直接给参数直接讲希望能帮你少走一些弯路。1. 拟时轨迹分析的核心思路与Monocle2选型逻辑1.1 Monocle2如何把细胞“捋”成一条线反向图嵌入原理Monocle2的核心逻辑不复杂但你必须先明白它算出来的不是“时间”而是细胞在表达状态空间里的距离。我们可以把每个细胞的数以万计基因表达量看成高维空间里的一个坐标点拟时分析要做的就是在这个高维点云里找出一条能串联起不同状态的路径。Monocle2采用反向图嵌入Reverse Graph EmbeddingDDRTree作为降维工具先在高维空间里学习一个树状结构的低维流形再把所有细胞按表达相似度投影到这棵树上。这棵树的根到叶子的距离就是细胞之间的“拟时距离”本质上是状态转换的进度而不是真实的昼夜时间。所以你会发现在同一个时间点取样不同细胞也可以被放在伪时间轴的不同位置因为它们在分化过程中所处阶段不一样。我最初刚接触Monocle2时总把轨迹上的顺序想成严格的先后序列这其实是个认识误区。DDRTree生成的路径可以是分叉的一个母节点可以分出两个子节点这就对应了细胞命运决定的分叉点。Monocle2的拟时轨迹能识别这种分叉所以在分析免疫细胞分化、肿瘤细胞干性演进、发育过程这类有“分叉”故事的课题里特别适用。1.2 为什么选Monocle2而不是Monocle3或Slingshot单说拟时分析市面上现在有不少新工具。Monocle3已经发布多年Slingshot、Scanpy里的diffmap也都有各自的优势。但Monocle2目前仍然有不可替代的使用场景尤其是做基因模块下游富集解析时它的差异基因检验框架和伪时间热图配合得最好。我用一张表把这几个工具的使用特点大致说清楚工具原理优点适合场景Monocle2DDRTree反向图嵌入轨迹连续稳定BEAM分支检验成熟热图可视化方便中小规模样本数万细胞以内需要提取伪时间相关基因模块Monocle3UMAP学习图大规模数据支持好分区模块更细超大样本量需要精细分区轨迹Slingshot主曲线拟合方法透明速度较快已有明确起点和聚类标签如果要追问为什么现在很多文献还在用Monocle2一个很现实的原因是它的伪时间基因模块提取流程非常顺。Monocle3虽然也能排序但下游做BEAM和伪时间热图时需要额外改造很多朋友的踩坑成本反而更高。Slingshot则强依赖你预设的起点和聚类注释对没有先验知识的流程不太友好。Monocle2的局限我也直说它对细胞数量有瓶颈超过十万个细胞后DDRTree降维会比较吃力而且它强制所有细胞被安排在一棵树上如果数据里存在相互独立的多条谱系分化Monocle2的处理效果会不那么理想。所以选型前需要先看清自己的数据形态。1.3 从表达矩阵到CellDataSet的数据准备Monocle2一切分析都建立在CellDataSet对象上。你手里常见的表达矩阵通常是基因×细胞的数值矩阵首先要把行名设为基因ID可以是symbol也可以是Ensembl ID列名设为细胞条形码。接下来需要给细胞和基因分别准备元数据。用代码来演示一下从矩阵构建CellDataSet的基本流程library(monocle) expr_matrix - readRDS(expr_matrix.rds) # 基因 x 细胞 cell_metadata - readRDS(cell_metadata.rds) # 细胞注释 gene_annotation - data.frame(gene_short_name rownames(expr_matrix)) rownames(gene_annotation) - rownames(expr_matrix) pd - new(AnnotatedDataFrame, data cell_metadata) fd - new(AnnotatedDataFrame, data gene_annotation) cds - newCellDataSet( expr_matrix, phenoData pd, featureData fd, lowerDetectionLimit 0.5, expressionFamily negbinomial.size() )其中expressionFamily参数的选择很重要。如果你的表达量是UMI count或者read count推荐用negbinomial.size()如果是经过大小归一化后的TPM/FPKM这类连续值则需要用gaussianff()。我见过不少朋友在构建时直接用默认的negbinomial导致后面estimateSizeFactors和estimateDispersions报错或者结果异常这一步务必检查清楚。数据准备阶段的核心目的是让后续分析输入干净把低表达的基因留下太多会增加降维噪音把细胞过滤得不干净轨迹会出假分支。按我的经验至少要先过滤掉在极少数细胞中表达的基因。1.4 必须做的事过滤、归一化与离散度估计在我实际跑过的流程里正式降维之前有四步是雷打不动的detectGenes、estimateSizeFactors、estimateDispersions以及选高变基因。cds - detectGenes(cds, min_expr 0.1) expressed_genes - row.names(subset(fData(cds), num_cells_expressed 10)) cds - cds[expressed_genes, ] cds - estimateSizeFactors(cds) cds - estimateDispersions(cds)这里num_cells_expressed 10是我常用的经验值但具体的阈值要根据你的细胞总数和测序深度来调整。细胞数越多阈值可以适当提高测序深度越低阈值要放宽避免丢失真实表达模式。之后的降维有些教程会先选高变基因有些则直接对全部过滤后的基因跑reduceDimension。我的建议是先做一步差异表达基因筛选再由这些基因参与DDRTree降维。这样既能保留有用的生物学信号又能让运算速度加快很多尤其是细胞数量在几万级别时全基因跑DDRTree是真的会等到怀疑人生。2. 基因模块提取从拟时轨迹到候选基因2.1 判断轨迹形态有分支原理没分支走另一套逻辑拿到一条轨迹后第一件事不是立刻提取基因而是看轨迹有没有分叉点。你可以用plot_cell_trajectory(cds, color_by Pseudotime)观察细胞在树上的分布也可以用plot_complex_cell_trajectory查看分支结构。Monocle2的分支节点在内部叫Branch Point是轨迹中某个细胞状态可以朝向两个不同方向发展的位置。如果你的轨迹只有一个分支比如细胞从A状态分化为B和C两种终末状态就可以用BEAMBranch Expression Analysis Modeling来检验哪些基因沿着从分支点到终点的表达模式有显著变化。这类基因往往正好对应了细胞命运决定的关键调控子。但如果你的数据没有明显分支轨迹基本是一根连续的线那么就不要强行用BEAM改用differentialGeneTest以伪时间作为协变量来筛选随时间显著变化的基因。2.2 BEAM与differentialGeneTest的参数细节先说说有分支的情况BEAM的代码非常简单但参数值得仔细讲一下cds - orderCells(cds) BEAM_res - BEAM(cds, branch_point 1, cores 8) BEAM_res - BEAM_res[order(BEAM_res$qval), ] BEAM_genes - subset(BEAM_res, qval 1e-4)其中branch_point 1需要你确认在使用orderCells后细胞被正确分配到树枝上。如果树枝顺序不对BEAM跑出来的基因分会乱掉导致后续模块完全失去生物学意义。我习惯在跑BEAM前先用plot_cell_trajectory(cds, color_by State)检查State分配和分支结构是否合理。cores参数影响计算速度如果BEAM运行时间过长多半是这里没设置好。加上cores 8之后实际运行速度会有质的提升。对于没有分支的轨迹用differentialGeneTestdiff_test - differentialGeneTest(cds, fullModelFormulaStr ~sm.ns(Pseudotime)) diff_test - diff_test[order(diff_test$qval), ] diff_genes - subset(diff_test, qval 0.01)~sm.ns(Pseudotime)是Monocle2里默认使用的自然样条基函数用来拟合伪时间的非线性表达趋势。如果不加sm.ns模型拟合的是线性关系很多表达先升后降或者波动型基因就会漏掉。筛选阈值项目中常见用的是qval 0.05但如果你后续还要做模块聚类我会建议适当收紧到qval 1e-4甚至更严格。因为模块聚类对噪音基因很敏感宽松阈值会把大量不相关基因拉进模块最后GO富集结果非常分散难以解释。2.3 把筛选出的基因按表达趋势聚成模块拿到差异基因列表后Querier通常会用plot_pseudotime_heatmap直接可视化这个函数内部会按伪时间把基因表达拟合成曲线然后用k-means聚类。代码长这样plot_pseudotime_heatmap( cds[candidate_genes, ], num_clusters 4, cores 8, show_rownames TRUE )这个函数输出的热图非常直观每一行是一个基因列是按伪时间排序的细胞行方向按照聚类模块分组能立刻看出来哪些基因在分化早期高表达、哪些基因在终末阶段才被激活。不过我建议不要只用这一个函数如果直接在R里跑plot_pseudotime_heatmap返回的聚类结果比较“一次性”不方便直接拿到每个模块的基因列表做后续GO富集。我常用的方法是手动提取表达趋势矩阵并自行聚类library(tidyverse) pseudo_expr - data.frame(t(exprs(cds[candidate_genes, ]))) pseudo_expr$Pseudotime - cds$Pseudotime tet - pseudo_expr %% pivot_longer(-Pseudotime, names_to gene, values_to expr) # 这里更推荐先用smooth.spline拟合每个基因的表达趋势 # 再把拟合值按伪时间排序做scale归一化 gene_matrix - vapply( candidate_genes, function(g) { fit - smooth.spline(pseudo_expr$Pseudotime, pseudo_expr[[g]]) predict(fit, x order(pseudo_expr$Pseudotime))$y }, numeric(nrow(pseudo_expr)) ) gene_matrix_scaled - t(scale(t(gene_matrix))) km - kmeans(gene_matrix_scaled, centers 4, nstart 20) module_list - split(names(km$cluster), km$cluster)这么做的好处是你能把每个模块的基因列表直接存成文本文件方便批量导入clusterProfiler之类的工具。另外先做趋势拟合再做聚类可以有效去除单细胞表达矩阵里的采样噪声如果直接对原始表达量做聚类会发现同一模块里的基因表达波形很碎模块划分极不稳定。2.4 模块数量怎么选宁可少不要贪多关于num_clusters的取值很多新手会在4还是在6之间反复横跳。我自己的经验是模块数量不宜过多每个模块里的基因数至少要保持在20到50个以上。过少基因的模块做GO富集时统计功效会严重不足经常会得到一堆语义宽泛但毫无区分度的通路对解释细胞状态转变帮助不大。一开始可以把模块数设成3到6然后观察每个模块基因的平均表达曲线是否具备不同的时序特征。如果两个模块的表达趋势几乎重合说明聚类数设大了如果一个模块里既有早高表达基因又有晚高表达基因说明聚类数偏小。热图配合平均曲线检查比任何统计指标都直观。模块本身也不是越少越好。基因太少则富集无意义基因太多则模块内功能混杂富集出来一堆本体论大词难以提炼核心结论。在实际项目里我通常会在num_clusters 5的基础上根据富集结果微调最终选定4到6个模块这样既有清晰生物学故事下游分析也足够稳。3. GO富集分析把模块变成功能注释3.1 选对工具与基因ID模块基因列表准备好之后就到了整个流程的收尾环节GO富集。这一环节的目标是搞清楚每个模块里的基因到底参与了哪些生物过程、分子功能和细胞组分。工具层面我强烈推荐R包clusterProfiler它支持批量输入多个模块也能直接可视化省去切换网页工具的麻烦。但要用它首先得把基因ID统一成注释包支持的格式。如果是人类数据需要org.Hs.eg.db小鼠数据用org.Mm.eg.db斑马鱼、拟南芥等也有对应的OrgDb包只是支持范围不一定全。如果你的物种注释偏冷门也可以考虑用DAVID和g:Profiler做在线富集但那样不方便对多个模块做统一管理和批量可视化。基因ID转换有一个细节clusterProfiler用了不同的keyType而很多单细胞注释文件默认给的是gene symbol。你要先用bitr转成ENTREZID再执行enrichGO。千万别把gene symbol直接丢进enrichGO的gene参数然后指望着readable TRUE帮你解决一切那样通常会报错或结果为空。3.2 一个可以直接跑的clusterProfiler流程下面这段代码是我在项目里反复使用的标准流程适用于人类数据其它物种只需替换OrgDb和相应的keyTypelibrary(clusterProfiler) library(org.Hs.eg.db) module_genes - readLines(module_1_genes.txt) gene_entrez - bitr( module_genes, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db ) # 如果担心背景基因问题可以计算所有检测到基因的Entrez列表 all_genes - readLines(all_detected_genes.txt) all_entrez - bitr( all_genes, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db ) ego - enrichGO( gene gene_entrez$ENTREZID, universe all_entrez$ENTREZID, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.2, readable TRUE ) head(egoresult)这里universe参数值得单独说明。GO富集的本质是超几何检验它需要知道“你从多少个基因里挑出当前这些候选基因”。如果universe不指定clusterProfiler默认用注释包里所有基因这会让你的显著性偏向于那些在基因组里本身就和某通路紧密相关的基因导致很多模块富集到相同的“基础代谢”或“拼接”通路丢失特异性。在实际分析时我认为更合理的universe是在你的单细胞数据里经过detectGenes过滤后仍然被检测到的所有基因。这样富集比较的基线才是“你能检测到什么”而不是“基因组上有什么”结果才更贴近真实生物学。3.3 多模块批量富集与可视化既然做了多个基因模块逐个跑enrichGO效率太低了。clusterProfiler提供了一个非常方便的函数compareCluster直接把模块列表打包处理module_list - list( module1 readLines(module_1_genes.txt), module2 readLines(module_2_genes.txt), module3 readLines(module_3_genes.txt) ) ck - compareCluster( module_list, fun enrichGO, OrgDb org.Hs.eg.db, keyType SYMBOL, ont BP, pvalueCutoff 0.05, qvalueCutoff 0.2, readable TRUE ) dotplot(ck, showCategory 5)compareCluster会自动把每个模块当作一个输入组返回一个合并的富集结果再用dotplot绘制点图每个模块一行点的大小表示富集到的基因数目颜色深浅表示显著性。这么做的好处是可以快速比较不同模块在功能上的差异比如模块1富集到细胞周期模块2富集到免疫应答模块3富集到代谢那整个细胞状态演进的逻辑骨架就出来了。如果你在dotplot输出里发现所有模块都富集到同一批通路不要急着往下写结论。先检查模块划分是否有问题再检查universe是否设置合理最后检查是不是差异基因列表有大量重叠。这三点经常是结果看起来“千篇一律”的根源。3.4 富集结果表格怎么看GO富集的结果表格里我平时重点看四列Description、GeneRatio、p.adjust和geneID。Description是功能通路名字比如“regulation of cell differentiation”“response to cytokine”GeneRatio表示你的模块基因里有多少个落在这个通路分母是模块总基因数分子是富集到该通路的基因数p.adjust是多重检验校正后的显著性一般小于0.05算有意义geneID列是具体的基因列表末尾如果显示/分隔的多个基因可以用来核对该通路是否真的和你想讲的故事有关。很多新手只看p值不看GeneRatio容易被一个只富集到2个基因的通路吸引住然后写出一大段解释。我的习惯是如果富集到某个通路的基因数少于5即使p值再小也只作为辅助参考不放到核心结论里。少于5个基因的富集往往经不起抽样的随机性考验换个样品或换个聚类参数可能就消失了。4. 结果解读、避坑与效率优化4.1 从GO术语反推生物学意义GO富集结果的解读本质上是要把轨迹上的基因模块和细胞命运决定事件联系起来。拟时轨迹和模块分析本身只负责“排序”和“找基因”真正告诉你怎么回事的是功能富集。举个例子假设你的轨迹是从naive T细胞到Th1/Th2分化分支点处BEAM筛出的差异基因会被聚成两个模块分支点前模块主要富集到T细胞受体信号通路、细胞激活相关过程分支点后某个模块富集到细胞因子介导的信号通路等。这时你就不能说“这些基因只在分化后期表达”而要结合模块表达曲线和分支顺序去讲哪一批基因先启动作出决定哪一批基因在终末阶段执行功能。我一直提醒自己和身边同事GO富集得到的生物学过程只揭示了“这一类基因倾向参与什么”而不等于“细胞一定发生了什么”。比如富集到“cell cycle”不代表细胞真的在增殖可能只是模块里包含了若干周期相关基因实际表达趋势未必和细胞周期状态同步。这个分寸感很重要尤其是写文章时千万别把相关性说成因果性。4.2 实操中最容易翻车的5个细节第一branch_point参数没确认。BEAM输出的p值是针对某个指定分支的如果轨迹有多个分支而你把branch_point设成1结果可能只反映了第一个分支的变化。最好先用plot_complex_cell_trajectory看节点排列再逐一跑BEAM。第二聚类前没有做标准化。不同基因的表达量本身就差异巨大高表达基因会主导欧氏距离导致聚类结果里高表达基因自成一团低表达但趋势重要的基因被混在一起。用scale按行标准化再聚类几乎是必须的。第三GO分析的universe没设置。如果没有指定背景显著性会被所有注释基因稀释模块差异性的信号会被冲淡。一定要把universe设为你实际检测到的基因集。第四基因ID转换不完整。bitr有时会返回NAs这是因为有新旧symbol不一致、线粒体基因、核糖体RNA等情况。建议转换后过滤掉NA再富集不要强行把NA丢进去跑。第五富集结果为空就盲目放宽阈值。我理解没有结果时想放宽一点的心理但盲目把qvalueCutoff调到0.9只会得到一堆无意义的宽泛通路。更理性的做法是检查模块基因数量是不是太少、universe是否合理、ont是否选得合适实在不行就换用MSigDB或者Reactome富集作为补充。4.3 常见报错速查表整理一个我在答疑中经常碰到的报错速查表方便你排查报错信息可能原因解决办法Error in BEAM: cds must have a branch轨迹没有分叉改用differentialGeneTestError in estimateDispersionsexpressionFamily使用错误UMI/count用negbinomial.size()TPM/FPKM用gaussianff()Error in bitr: some genes cannot be mapped基因名规范或物种匹配问题检查基因标识类型调整fromType过滤映射不上的基因No enrichment found in any gene set模块基因过少、universe过大或者阈值太严增加模块基因、缩小universe、合理放宽pvalueCutoffKilled或内存不足BEAM计算量过大用cores并行减少参与基因数或细胞数还有一个容易被忽略的地方单细胞矩阵里的基因名如果是Ensembl ID而你的差异基因列表用的是symbol两者没对齐时GO富集大概率会全军覆没。处理这种问题我通常在样品注释文件生成阶段就把基因名统一成同一种格式宁可在上游多花十分钟也不要在下游排查半天。4.4 效率优化怎样让流程跑得更顺手Monocle2最大的槽点是跑得慢尤其是BEAM。如果你的细胞数量在几万以上基因数量又在几千BEAM可能跑数小时。为了兼顾效率和结果稳定性我一般会采用“先粗筛再精算”的策略先用differentialGeneTest以比较宽松的阈值筛出一批候选基因再做BEAM因为BEAM的模型参数估计远比普通差异检验复杂基因数量少一半时间能省一大半。GO富集阶段clusterProfiler本身速度不算慢但对多个模块跑多次时尽量用compareCluster合并执行这样只用加载一次注释包内存占用和耗时都更友好。另外如果你的OrgDb包版本很旧会导致部分基因ID映射失败建议定期更新org.*.db包别让环境里的旧注释拖后腿。写在最后的一点心里话这整套流程我自己跑过很多遍最大的体会是Monocle2拟时轨迹只是给你一个细胞状态变化的抽象骨架真正让这个骨架长出肉、让生物学故事立住的是后续的基因模块和GO富集解析。很多朋友把70%的时间花在跑轨迹上留给功能分析的时间很少这其实有点本末倒置。如果让我回头重做一次我会先把数据质控和聚类前期工作做得更扎实再跑拟时拿到模块列表后优先用GO富集和文献相互印证而不是一味追求更小的p值。最后再分享一个小习惯把每个模块的基因列表单独存成文本文件旁边附带一列模块编号和一列简单注释。这样不管是后续换用Reactome、MSigDB还是做体外实验验证都能随时取用。看似不起眼但在项目周期拉长到几个月之后你会感谢当初这个微不足道的整理动作。
