企业官网建设流程全解析

Monocle2的拟时分析做完之后，很多人停在了“画出一条漂亮的轨迹图”这一步——细胞按颜色从浅到深排开，伪时间轴铺在树状结构上，确实好看。但一张轨迹图只能证明“分化确实发生了”，它回答不了“沿着这个轨迹，细胞的哪些生物学程序在改变”这个问题。要把轨迹和功能联系起来，就必须走到基因模块层面：先让数千个随拟时变化的基因按表达趋势聚成几组，再对每一组做GO富集解析，才能把“这个模块在干嘛”翻译成生物学结论。这篇文章把我自己跑这套流程的完整思路、代码细节和踩过的坑整理出来，从Monocle2的拟时值提取、基因模块划分，到GO富集和结果解读，一条线讲清楚。适合正在做单细胞分化轨迹分析、对Monocle2已经基本入门但卡在“轨迹之后做什么”的同学参考。

1. 为什么拟时分析要走到“基因模块”这一步

1.1 单基因轨迹的三个常见瓶颈

最早我做拟时分析时，习惯把几个明星marker基因拉出来，用plot_genes_in_pseudotime画三张趋势曲线，比如干性基因下降、分化基因上升，然后结论就写完了。但做多了就会发现，这条路有三个瓶颈。

第一，单基因只能回答“这一个基因是否随时间变化”，回答不了“这一群基因是否协同变化”。拟时轴上真正重要的生物学过程，比如糖酵解重编程、线粒体代谢转换、细胞周期退出，动辄涉及几百个基因协同表达。你挑三五个marker，只是管中窥豹。

第二，成千上万个基因都在随拟时波动，靠肉眼挑根本不现实。differentialGeneTest通常会返回一两千个显著基因，一个人手工去看这些基因的趋势，既低效又主观。

第三，也是最关键的：单基因没有“通路”语义。哪怕你看到基因A在后期显著升高，你也很难立刻说清楚这代表什么过程；但如果是一组基因同时富集到“氧化磷酸化”，结论就非常明确。

这就是为什么需要基因模块。所谓基因模块，我理解为一群表达趋势高度相似的基因：在拟时早期高表达，或者在拟时后期骤然升高，或者中间出现一个峰值。把这些基因聚成4到8个模块之后，每个模块内部基因的行为一致，模块和模块之间又互相区分，于是每一段拟时上发生的“程序切换”就变成了几个模块之间的接力赛。

1.2 基因模块和GO富集的配合逻辑

基因模块划分出来以后，下一步自然就是功能注释。模块本身是一串基因ID，如果不做富集，它就是无名无姓的几十几百个符号。GO富集做的事情很简单：拿模块里的基因和一个背景基因集比较，用超几何检验看哪些基因本体论条目（Gene Ontology，GO条目）在这批基因里出现得异常多，比如“细胞周期”“DNA复制”“突触传递”等等。

GO富集是给模块“贴标签”的过程。有了标签，你才能在论文里写“module 1在早期富集到细胞周期相关条目，提示细胞在拟时早期仍处于增殖状态；module 4在后期富集到神经递质传递相关条目，与功能成熟阶段吻合”。所以完整流程是一条链：Monocle2拟时轨迹 → 拟时相关基因 → 基因模块 → GO富集 → 生物学解读。第2节先把Monocle2这一段快速过一遍，因为拟时值的质量决定后面一切。

2. 从表达矩阵到拟时值：Monocle2全流程回顾

2.1 输入对象构建：一步错步步错的CellDataSet

Monocle2的所有操作都围绕一个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_components=2并不代表只保留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.05，qvalueCutoff = 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-release@monocle2")

装完之后加载包名仍然是library(monocle)，不影响使用。如果你的集群上R版本太老，也可以去Bioconductor的归档版本里找对应R版本的旧二进制包。另外，Monocle2和Monocle3不要混着装，两者的函数名和对象结构有冲突，我见过同事同时加载两个gem导致newCellDataSet都调不对的情况。

这套流程跑下来，你手里的产物不只是一张轨迹图，而是一组“模块+GO条目”的对应关系。写论文时，表里可以放每个模块的top GO；讨论里可以按“早期模块→中期模块→后期模块”的顺序把分化过程的程序切换讲清楚。我个人在最近的几个项目里，都靠这套流程在审稿时少挨了不少关于“你的标记基因选择为什么主观”的质疑——因为基因模块和GO富集至少在聚类和统计意义上是有可复现依据的。希望这篇整理的实操笔记能帮你少走几步弯路。