Monocle3 轨迹分析 — 单细胞发育轨迹与伪时间推断工具¶

一句话说明¶

Monocle3 v1.4.26 是 R 语言中最权威的单细胞轨迹分析工具，通过"伪时间"（pseudotime）将细胞排列在发育路径上，揭示分化、发育、激活等动态过程。

安装与配置¶

# 安装 Bioconductor 依赖（R 4.4+ 环境）
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

# 安装 Monocle3 依赖包
BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                        'limma', 'lme4', 'S4Vectors', 'SingleCellExperiment',
                        'SummarizedExperiment', 'batchelor', 'HDF5Array',
                        'terra', 'ggrastr'))

# 从 GitHub 安装 Monocle3（当前最新版 v1.4.26）
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/monocle3')

# 验证安装
library(monocle3)
packageVersion("monocle3")   # 应显示 "1.4.26" 或以上

核心用法¶

library(monocle3)              # 加载 Monocle3
library(ggplot2)               # 可视化
library(dplyr)                 # 数据操作

# ── 创建 CellDataSet 对象 ─────────────────────────────
# Monocle3 使用 CellDataSet (cds) 对象，相当于 Scanpy 的 AnnData

# 从表达矩阵创建（行为基因，列为细胞）
expression_matrix <- readRDS('expression_matrix.rds')   # 读取表达矩阵
cell_metadata <- readRDS('cell_metadata.rds')           # 细胞元数据（含聚类等信息）
gene_metadata <- readRDS('gene_metadata.rds')           # 基因元数据

# 创建 CellDataSet 对象
cds <- new_cell_data_set(
  expression_matrix,          # 原始计数矩阵
  cell_metadata = cell_metadata,
  gene_metadata = gene_metadata
)

# ── 预处理 ────────────────────────────────────────────
cds <- preprocess_cds(
  cds,
  num_dim = 50,               # PCA 维数（类似 Scanpy 的 n_comps）
  norm_method = "log"         # 标准化方法
)

# ── 降维 ──────────────────────────────────────────────
cds <- reduce_dimension(
  cds,
  reduction_method = 'UMAP', # 降维方法（UMAP 或 tSNE）
  preprocess_method = 'PCA'  # 基于 PCA 做 UMAP
)

# ── 聚类 ──────────────────────────────────────────────
cds <- cluster_cells(
  cds,
  resolution = 1e-5           # 聚类分辨率
)

# 可视化聚类结果
plot_cells(cds, color_cells_by = "cluster")

参数详解¶

参数	函数	说明
`num_dim`	`preprocess_cds`	PCA 主成分数（建议 30-100）
`resolution`	`cluster_cells`	聚类分辨率（越大 cluster 越多）
`min_branch_len`	`learn_graph`	轨迹最短分支长度
`use_partition`	`learn_graph`	是否按聚类学习局部轨迹
`close_loop`	`learn_graph`	是否允许轨迹形成环路

实战案例¶

library(monocle3)
library(ggplot2)

# ── 完整轨迹分析流程 ──────────────────────────────────

# 1. 从 Seurat 对象转换（最常见的工作流）
library(Seurat)
seurat_obj <- readRDS('seurat_annotated.rds')  # 读取已注释的 Seurat 对象

# 提取数据构建 CDS 对象
gene_annotation <- data.frame(
  gene_short_name = rownames(seurat_obj),       # 基因名
  row.names = rownames(seurat_obj)
)
cds <- new_cell_data_set(
  seurat_obj@assays$RNA@counts,                # 原始计数矩阵
  cell_metadata = seurat_obj@meta.data,         # 细胞元数据
  gene_metadata = gene_annotation
)

# 2. 使用 Seurat 的 UMAP 坐标（保持一致性）
# 先在 Monocle3 内部预处理
cds <- preprocess_cds(cds, num_dim = 30)
cds <- reduce_dimension(cds)

# 将 Seurat UMAP 坐标嵌入 Monocle3
reducedDims(cds)$UMAP <- seurat_obj@reductions$umap@cell.embeddings
# 将 Seurat 聚类信息传入
cds@clusters$UMAP$clusters <- factor(seurat_obj@meta.data$seurat_clusters)

# 3. 学习轨迹图（PAGA 图主干）
cds <- learn_graph(
  cds,
  use_partition = TRUE,        # 在各 partition 内学习局部轨迹
  verbose = TRUE
)
plot_cells(
  cds,
  color_cells_by = "cell_type",  # 按细胞类型着色
  label_groups_by_cluster = FALSE,
  label_leaves = FALSE,
  label_branch_points = FALSE
)

# 4. 设置轨迹起点（需要选择祖先细胞所在的节点）
# 交互式选择起点（在 RStudio 中使用）
# cds <- order_cells(cds)

# 脚本中自动指定起点（根据已知祖先细胞类型）
root_cell <- rownames(colData(cds))[colData(cds)$cell_type == "Progenitor"][1]
cds <- order_cells(cds, root_cells = root_cell)   # 指定根细胞

# 5. 可视化伪时间
plot_cells(
  cds,
  color_cells_by = "pseudotime",    # 按伪时间着色（从蓝到黄代表从早到晚）
  label_cell_groups = FALSE,
  label_leaves = FALSE,
  label_branch_points = FALSE,
  graph_label_size = 1.5
)

# 6. 寻找与伪时间相关的基因（轨迹差异表达）
pr_test_res <- graph_test(
  cds,
  neighbor_graph = "principal_graph",  # 基于轨迹主图检验
  cores = 4                            # 使用 4 个 CPU 核心
)
# 筛选 p 值显著的基因
pr_deg_ids <- row.names(
  subset(pr_test_res, q_value < 0.05)  # FDR 校正后 p 值 < 0.05
)

# 7. 可视化基因表达沿伪时间的变化
plot_cells(
  cds,
  genes = head(pr_deg_ids, 4),       # 展示前 4 个显著基因
  show_trajectory_graph = FALSE,
  label_cell_groups = FALSE
)

# 8. 模块分析（寻找共表达基因模块）
gene_module_df <- find_gene_modules(
  cds[pr_deg_ids,],                  # 只分析显著的差异基因
  resolution = 0.1
)

常见报错与解决¶

报错	原因	解决方法
`Error: BiocManager not found`	缺少 Bioc 管理器	`install.packages("BiocManager")`
UMAP 坐标替换后轨迹不对	维度顺序不匹配	检查 UMAP 列名为 `UMAP_1, UMAP_2`
`learn_graph` 运行很慢	细胞数太多	先 subsample 到 5 万细胞
伪时间全是 Inf	起点细胞选错	重新用 `order_cells` 选择合适起点
`graph_test` 报内存错误	基因数太多	只分析高变基因

速查表¶

# 完整轨迹分析 5 步走
cds <- preprocess_cds(cds, num_dim=30)   # 步骤1：PCA预处理
cds <- reduce_dimension(cds)              # 步骤2：UMAP降维
cds <- cluster_cells(cds)                # 步骤3：聚类
cds <- learn_graph(cds)                  # 步骤4：学习轨迹
cds <- order_cells(cds)                  # 步骤5：排序/赋伪时间

# 可视化命令
plot_cells(cds, color_cells_by="pseudotime")   # 伪时间图
plot_cells(cds, color_cells_by="cell_type")    # 细胞类型图