企业官网建设流程全解析

这阵子把单细胞monocle3的分析流程重新过了一遍，起因很实际：手头项目需要从细胞聚类做到发育轨迹推断，旧笔记里那套流程在Monocle3几次版本更新之后已经有好几处跑不通了。这次整理比上次顺畅不少，但也确实踩了几个新坑，尤其是数据入口和根节点选择这两块，网上教程往往一两句话带过，真上手全是问题。我把从构建cds对象到聚类、轨迹、差异表达、可视化的完整流程重新梳了一遍，每步都标注了参数选择和背后的逻辑，如果你已经会用Seurat做基础单细胞分析，或者正打算给GEO下载的数据补一条轨迹分析，这篇基本能让你把Monocle3全流程走通。

1. 为什么Monocle3的流程需要“再整理”

1.1 版本迭代带来的API变化

很多新接触Monocle3的朋友会去翻早期教程，结果翻到一大半就发现代码跑不动。原因很简单：Monocle3和Monocle2是完全两套设计。Monocle2的核心是DDRTree降维加反向图嵌入，伪时间计算也依赖这一套；Monocle3改成了先UMAP降维，再在UMAP结构上学习主图，最后基于主图计算伪时间。思路变了，接口自然大改。早期教程里常见的orderCells、setOrderingFilter、detectGenes这些函数，在Monocle3里已经被order_cells、preprocess_cds、cluster_cells等一批新函数替代。新旧混用的结果就是一连串“could not find function”之类的报错，光排查这些问题就能耗掉半天。

1.2 工具定位认知要更新

Monocle3不是只能做轨迹分析。它从new_cell_data_set构建对象开始，preprocess_cds做预处理，reduce_dimension降维，cluster_cells聚类，learn_graph学图，graph_test做差异表达，plot_cells可视化，每个环节都有对应的函数，完全能作为一条独立的单细胞分析流程使用。只不过在实际项目中，大家还是习惯先用Seurat做前期的质控、聚类、注释，再用Monocle3专门做轨迹推断和轨迹相关的差异基因分析。

还有一个常见的认知偏差：以为只有整个数据集都要跑轨迹。真实场景中，轨迹分析的对象往往是“某一个感兴趣的亚群”，比如从全部细胞中先注释出T细胞亚群，再在这个子集上重新构建cds做轨迹分析。全量细胞上跑不是不行，但partition一多，learn_graph学出来的图会很乱，而且计算时间大幅上升，后续解释也会很费力。

2. 数据入口：如何正确构建cds对象

2.1 new_cell_data_set的三个核心输入

Monocle3的数据对象叫cell_data_set，简称cds，构建入口是new_cell_data_set函数。它有三个关键输入：expression_data、cell_metadata、gene_metadata。

expression_data是表达矩阵，行是基因，列是细胞。这里有一个非常关键的要求：传入的一定要是原始counts，而不是normalized后的数据，更不能是scaled数据。preprocess_cds在内部会自己做normalization和log转换，如果你在外面已经把数据标准化了一遍，再传进去，相当于做了两轮处理，聚出来的结果基本不能用。我在重跑流程时就犯过这个错：从Seurat对象里直接取了@assays$RNA@data（也就是normalized data）去构建cds，结果UMAP图怎么看怎么不对劲，后来换成@assays$RNA@counts才恢复正常。

cell_metadata是细胞信息表，行名必须和表达矩阵的列名一致，里面可以放样本编号、批次、细胞类型注释这些信息。gene_metadata是基因信息表，行名必须和表达矩阵的行名一致，特别注意里面一定要有一列叫gene_short_name，存的是基因名。很多教程在构建cds时没有单独准备这一列，结果后面plot_cells画基因表达图时直接报错，提示找不到gene_short_name。

2.2 从Seurat对象转换的实操写法

假设你手里的数据已经在Seurat对象里，常规做法是先从Seurat中提取counts矩阵，再带入new_cell_data_set。推荐写成这样：

library(monocle3) library(Seurat) obj <- readRDS("seurat_annotated.rds") expr_matrix <- GetAssayData(obj, assay = "RNA", slot = "counts") cell_meta <- obj@meta.data gene_meta <- data.frame(gene_short_name = rownames(expr_matrix)) rownames(gene_meta) <- rownames(expr_matrix) cds <- new_cell_data_set(expr_matrix, cell_metadata = cell_meta, gene_metadata = gene_meta)

这里有个细节容易被忽略：cell_meta的行名必须是细胞barcode，而且和expr_matrix的列名顺序可以不一致，new_cell_data_set会自动按行名对齐。但如果在构造cell_meta时不小心让行名变了，后面就会出现维度对不上的问题。建议转换后马上验证一下：

all(rownames(cell_meta) == colnames(expr_matrix)) all(rownames(gene_meta) == rownames(expr_matrix))

返回TRUE才继续往下走。这个习惯能帮你省掉很多诡异的报错。

2.3 从10X标准输出构建cds

有些项目数据是从GEO下载来的10X标准格式，也就是三个文件：matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz。这时通常先用Seurat的Read10X函数读出矩阵，再走上面的转换流程。如果是.h5文件，则用Seurat的Read10X_h5读取。Monocle3本身也提供了读取10X产物的小工具，但实际用下来，先转成Seurat对象再转cds是最省事的路径，因为可以先把Seurat里的质控和注释信息一起带过来，省得在cds阶段重新整理metadata。

如果数据只有表达矩阵，没有现成的metadata，就先构建一个以细胞名为行名的minimal data.frame：

cell_meta <- data.frame(row.names = colnames(expr_matrix))

后面需要补充样本、分组信息时，直接往cds的colData里加列就行。Monocle3的colData和Seurat的meta.data类似，是一个DataFrame对象，可以直接用$符号加新列。

3. 预处理、降维与聚类：把细胞分群做扎实

3.1 preprocess_cds里被低估的参数

构建完cds，下一步是preprocess_cds。这个函数默认用PCA降维，然后为后面的UMAP做准备。最常见的问题是两个：维度取多少，以及要不要做批次校正。

num_dim控制PCA保留的主成分数，默认是100。对大多数数据集，50左右是一个比较稳的起点。太小会丢失太多信息，太大容易引入噪声。一个简单的方法是看PC的方差贡献曲线，找到拐点，也可以用下游聚类结果的稳定性来调。

更重要的一个参数是residual_model_formula_str。这个参数的作用是把某些协变量（比如测序深度、线粒体比例、批次信息）从表达矩阵中回归掉，再做降维。比如你的数据来自多个样本或者多个测序批次，不处理的话，UMAP上首先分出来的可能就是批次，而不是生物学差异。一个常见的写法是：

cds <- preprocess_cds(cds, num_dim = 50, residual_model_formula_str = "~batch + n.umi")

这里面的n.umi是Monocle3默认记录的每细胞UMI数。如果你的metadata里没有这个字段，就要换成自己实际有的列名，比如“percent.mt”或者“Sample”。说实话，这个参数是Monocle3里最容易被忽略但影响最大的参数之一。很多人聚出来的群和样本信息高度重合，大概率就是没做这一步。

3.2 reduce_dimension与cluster_cells的分群参数

降维和聚类在Monocle3里是分开的两步。reduce_dimension默认使用UMAP，主要可调参数包括umap.metric和umap.n_neighbors。umap.n_neighbors默认30，值越小局部结构越突出，值越大全局结构越明显。如果样本量很大，可以考虑设成50或者更高。做轨迹分析的话，我一般倾向于让UMAP图保持局部结构清晰，这样learn_graph学出来的主图更容易贴合真实的细胞状态变化。

聚类函数是cluster_cells，内部用的是Leiden社区发现算法，有一个关键参数resolution，默认值是1e-3。这个参数控制聚类的粗细，数值越大分出的群越多。具体调多少要根据自己的数据看，没有绝对标准。我的习惯是先跑一个中间值，然后在UMAP图上染色看分群是否合理，再结合marker基因的表达验证每个群的生物学意义。一个简单可用的做法：

cds <- cluster_cells(cds, resolution = 1e-3)

分群结束后，cluster信息可以通过clusters(cds)取出来，会返回一个以细胞名为名的向量。如果你想调整resolution重新聚类，不需要从头跑，直接再调用一次cluster_cells覆盖结果就行。要注意的是，重新聚类之后，后面learn_graph这一步的partition也会改变，得重新跑。

3.3 聚类后马上要做的事：手动注释验证

聚类本质上只是把细胞按表达谱分成群，下一步必须做细胞类型注释。这也是现在单细胞分析流程里最耗精力的环节。Monocle3没有内置的自动注释功能，通常做法是结合已知marker基因，在UMAP图上逐个查看。

在Monocle3里看marker基因非常方便，plot_cells可以直接在UMAP上同时展示多个基因的表达，不用来回切工具：

plot_cells(cds, genes = c("CD3D", "CD14", "CD79A", "LYZ"), label_cell_groups = FALSE)

这里给到一组免疫细胞经典marker：CD3D是T细胞，CD14和LYZ是单核/巨噬细胞，CD79A是B细胞。手动注释的要点是不要只看单个基因阳/阴就下结论，要看组合。比如T细胞常见CD3D阳性、CD14阴性；单核细胞是CD14阳性、CD3D阴性。此外，注释结果最好和Seurat那边的注释结果比对一下，两边一致的情况下，后面轨迹解释会踏实得多。如果你发现Monocle3聚出来的某个群marker表达特征模糊，不要硬命名，宁可标记成“unknown”，也比错标强。

4. 轨迹推断：learn_graph和order_cells的实操细节

4.1 learn_graph到底在学什么

轨迹推断是Monocle3的重头戏，核心函数是learn_graph。这一步是在UMAP坐标上学习一个主图，你可以把它理解成在细胞分布中找到一条或几条“主干道”，细胞沿着这些路径从一种状态过渡到另一种状态。

learn_graph有一个非常重要的参数use_partition，默认是TRUE。这个参数的意思是：聚类得到的每个大partition之间不强制连线。也就是说，如果两个细胞群在UMAP上离得很远，而且被聚类算法分到了不同partition，learn_graph不会强行把它们连起来。这个设计很合理，因为生物学上很多细胞类型之间本来就不存在连续分化关系，强行连线会产生假轨迹。

如果确认你的数据是一个连续分化过程，比如从干细胞到各系祖细胞的发育，可以尝试把use_partition设为FALSE。但这一步要非常谨慎，因为一旦关闭，learn_graph会把所有细胞都连到一个图里，非常容易产生看起来很“漂亮”但生物学解释不了的假轨迹。

另一个值得注意的参数是close_loop。默认情况下，learn_graph会尝试识别轨迹中的环状结构。如果你的数据里有类似细胞周期的过程，环状轨迹是有意义的；但如果只是普通的发育过程，出现了环状结构往往是降维或者聚类的假象。我的做法是：先保持默认跑一遍，看结果；如果UMAP上明显没有环状结构，但learn_graph画出了环，就设置close_loop = FALSE重新学一次。

4.2 根节点选择决定伪时间方向

轨迹学出来了，接下来要用order_cells定义伪时间的起点。根节点选哪里，直接决定下游所有细胞伪时间的排序，也就决定了你后面做差异基因分析看到的“趋势”方向。根选反了，原本向终末分化的基因看起来会像反向表达，后面解读就会出大问题。

Monocle3的order_cells支持几种指定根节点的方式。最简单的交互式操作是直接运行：

cds <- order_cells(cds)

这时会弹出一个UMAP图窗口，手动点击图中的某个位置作为根节点。点击的位置最好选在UMAP上处于“起始状态”的那群细胞附近，而不是任意一个角落。

另一种更可控的方式是指定细胞barcode。比如你手头有手动注释好的细胞类型，已知这群细胞是发育起点，就可以这样：

root_cell_barcodes <- colData(cds)$cell_type == "Stem" cds <- order_cells(cds, root_cells = rownames(colData(cds))[root_cell_barcodes])

这种方式的优点是明确、可重复，适合写进分析脚本里。但前提是你的注释结果足够可靠。如果对哪个群是起点没有把握，可以在UMAP上先看几个候选marker的表达，比如干性相关的基因（如Kit、Prom1、Cd34等）是否集中在某一群。哪群干性marker阳性，哪群大概率是起点。

伪时间计算完成后，可以用pseudotime(cds)取出每个细胞的伪时间值。此时在UMAP上染色检查一下：

plot_cells(cds, color_cells_by = "pseudotime")

如果伪时间从起点向外平滑扩散，说明根节点选得没问题。如果颜色分布很乱，或者起点不在预期位置，就需要重新选根。

4.3 分支与伪时间的生物学解读

learn_graph学出来的主图往往不是简单的直线，而是有分支的结构。分支点意味着在这个位置，细胞命运发生了分化，一部分走向A类细胞，另一部分走向B类细胞。

对于分支结构的解读，建议先用choose_graph_segments把感兴趣的分支选出来，把目标分支的细胞子集单独拿出来做后续分析。这样能避免全图的复杂结构干扰判断。选中分支后，可以对这个子集重新做graph_test，找到分支特异性的驱动基因。我以前处理过一个造血发育的数据，一开始在全图上跑graph_test，出来的基因列表覆盖了太多过程，很难聚焦。后来换成只对目标分支做分析，找到的分支特异基因明显更干净，也更符合文献报道。

5. 差异表达与可视化：把结果真正用起来

5.1 graph_test：找沿轨迹变化的基因

trajectory分析最终要落脚到基因层面，不然只有一条图，讲不了生物学故事。Monocle3的graph_test是专门用来找“随轨迹位置变化”的基因的工具。它的统计量是Moran's I，一种空间自相关指数。简单理解就是：一个基因的表达值如果在相邻的细胞之间都相似，而在远离的部位有差异，那它在轨迹上就有明显的空间自相关性，说明它的表达变化和伪时间或者分支位置有关。

用法很直接：

deg_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)

结果是一个data.frame，每一行是一个基因，关键列是morans_I和q_value。morans_I越接近1，说明基因表达在轨迹图上的空间聚集性越强；q_value是校正后的p值，通常用q_value < 0.05作为筛选阈值，可以结合morans_I排序，取靠前的基因往下游分析。

graph_test和常规的cluster差异表达不一样，它不关心某个基因在A群高还是在B群高，而是看表达量是否随着轨迹位置连续变化。这个设计特别适合用来找发育过程中渐变的关键调控因子，而不是仅仅在不同群之间跳跃表达的marker。

5.2 plot_cells与plot_genes_in_pseudotime的使用技巧

可视化部分，plot_cells是最常用的函数。它可以在UMAP上按不同方式着色：按cluster、按pseudotime、按基因表达量。

plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE)

如果想叠加查看某个基因的表达，把genes参数加上：

plot_cells(cds, genes = c("Kit", "Gata3"), color_cells_by = "pseudotime", label_cell_groups = FALSE)

此时图上会出现一个小面板展示这些基因在UMAP上的表达分布，能直观看到哪些基因和伪时间走向一致。

伪时间趋势图用plot_genes_in_pseudotime，但它的输入格式有点绕。必须先把基因放成一个data.frame，其中第一列是gene_short_name，再传入函数。一个示例：

gene_df <- data.frame(gene_short_name = c("Kit", "Gata3")) gene_df$gene_short_name <- as.character(gene_df$gene_short_name) plot_genes_in_pseudotime(cds, gene_group_df = gene_df)

出来的图每一行是一个基因，横轴是伪时间，纵轴是表达量，还会把拟合曲线画出来。这个图对判断基因在发育过程是先升后降还是持续上升很有帮助，写文章配图也常用。

6. 实战中的报错与问题排查记录

6.1 常见报错和对应的解决方法

重新整理流程的过程中，有几个报错反复出现，这里集中记录一下。

第一个是plot_cells报找不到gene_short_name。这个报错最常见的原因就是前面说的，构建cds时gene_metadata没有包含gene_short_name列。解决办法是在构建前就补好，或者对已有cds补一列：

rowData(cds)$gene_short_name <- rownames(cds)

第二个是learn_graph跑起来非常慢。如果数据量大，比如几万个细胞，learn_graph确实要花不少时间。可以先在子集上调试参数，确定没问题再跑全量。另外可以适当降低preprocess_cds的num_dim，或者增加机器的内存配额。跑learn_graph之前也可以先用choose_cells选择一部分感兴趣细胞来降载。

第三个是order_cells交互式选根时点的位置怎么也选不上。这通常是点击位置离主图太远。把图形窗口放大一些，尽量点在图上的节点附近，或者在UMAP图上先用plot_cells找到目标细胞群的中心位置，再点那里。

第四个是把Seurat对象转cds时报维度对不上。大多数情况是cell_meta的行名和表达矩阵的列名没对齐。处理办法很简单，构造metadata时就保证行名是barcode，换数据格式后先跑一遍前面写的验证代码。

6.2 根节点选错后的排查思路

伪时间方向和预期不一致，是轨迹分析里最隐蔽的问题。有时候不是报错，而是结果看起来“能分析”，但结论是反的。我的排查思路是这样的：先看UMAP上伪时间颜色的扩散方向，如果从终点细胞群开始向起点扩散，那多半是根选反了。这时候不要急着改代码，先用marker基因验证一下哪群是真正的起点细胞，再基于这群细胞的barcode重新order_cells。

还有一种情况比较麻烦：不同partition之间的伪时间没有可比性。因为learn_graph默认use_partition=TRUE，每个partition是分开学图的，伪时间在不同partition之间可能不连续。如果后续分析需要跨partition比较伪时间，就得跑一次use_partition=FALSE，或直接在一个选定的partition子集上分析。

改完根节点之后，记得重新跑graph_test，因为伪时间方向变了，空间自相关检验的结果也会相应变化。这个顺序搞错的话，报告里的基因列表可能完全是反的。

这次整理流程，我最大的体会是：Monocle3的代码本身不复杂，真正的门槛在参数怎么选、数据怎么入口、根节点怎么定。这些细节决定了分析结果的可靠性，也是网上教程最容易跳过的地方。如果你手头正好在跑单细胞数据，建议先用一个小规模子集把全流程跑通，再把参数固定下来跑全量，会顺手很多。后面有机会，我还会把分支特异性基因的筛选和模块分析再单独整理一篇，那部分也是实际操作中另一个容易卡壳的环节。