单细胞轨迹分析Monocle¶

一句话概述：Monocle3通过拟时序（pseudotime）分析将单细胞按发育或分化的先后顺序排列，揭示细胞状态转变的动态轨迹和分支决策点。

核心知识点速览¶

概念	白话解释
拟时序（Pseudotime）	把真实时间替换为"发育进度"，给细胞排队
轨迹推断	把散落的细胞点连成一条"发育路线图"
分支点	轨迹分叉的地方，代表细胞"命运选择"的节点
根节点（Root）	轨迹的起点，通常是干细胞或祖细胞
UMAP	降维可视化方法，把高维数据画到2D平面
Partition	Monocle3把不相连的细胞群分成独立区域
Graph	Monocle3构建的轨迹骨架（主图）
CDS	Cell Data Set，Monocle3的核心数据结构
差异表达基因	沿轨迹变化的基因（pseudotime DEGs）
Seurat转换	从Seurat对象转为Monocle3的CDS对象

一、Monocle3安装¶

# 安装Monocle3（需要从GitHub安装）
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("BiocGenerics", "DelayedArray", "DelayedMatrixStats",
                        "limma", "lme4", "S4Vectors", "SingleCellExperiment",
                        "SummarizedExperiment", "batchelor", "HDF5Array",
                        "terra", "ggrastr"))  # 安装依赖

install.packages("devtools")  # 安装devtools
devtools::install_github("cole-trapnell-lab/monocle3")  # 安装monocle3

library(monocle3)  # 加载

二、完整分析流程¶

2.1 从Seurat对象转换¶

library(Seurat)
library(SeuratWrappers)  # Seurat转换工具
library(monocle3)

# 假设你已有处理好的Seurat对象
# seurat_obj 已完成归一化、聚类和UMAP

# 方法1：用SeuratWrappers转换
cds <- as.cell_data_set(seurat_obj)  # Seurat → CDS

# 方法2：手动构建CDS
expression_matrix <- GetAssayData(seurat_obj, layer = "counts")  # 获取计数矩阵
cell_metadata <- seurat_obj@meta.data  # 获取细胞元数据
gene_metadata <- data.frame(
  gene_short_name = rownames(expression_matrix),  # 基因名
  row.names = rownames(expression_matrix)
)

# 创建CDS对象
cds <- new_cell_data_set(
  expression_data = expression_matrix,  # 表达矩阵
  cell_metadata = cell_metadata,         # 细胞信息
  gene_metadata = gene_metadata          # 基因信息
)

2.2 预处理和降维¶

# 步骤1：预处理（归一化和选择高变基因）
cds <- preprocess_cds(
  cds,
  num_dim = 50,              # PCA维度数
  method = "PCA"              # 降维方法
)

# 查看PCA解释方差
plot_pc_variance_explained(cds)  # 画碎石图，确定维度数

# 步骤2：去批次效应（可选）
cds <- align_cds(
  cds,
  alignment_group = "batch"  # 按哪个变量校正批次
)

# 步骤3：UMAP降维
cds <- reduce_dimension(
  cds,
  reduction_method = "UMAP",  # 使用UMAP
  preprocess_method = "PCA"   # 基于PCA结果
)

# 可视化
plot_cells(cds,
           color_cells_by = "cell_type",  # 按细胞类型着色
           label_cell_groups = TRUE,       # 显示标签
           cell_size = 0.5)                # 点大小

2.3 聚类和分区¶

# 聚类（Leiden算法）
cds <- cluster_cells(
  cds,
  resolution = 1e-3  # 分辨率参数（越小越少cluster）
)

# 查看聚类结果
plot_cells(cds, color_cells_by = "cluster")     # 按cluster着色
plot_cells(cds, color_cells_by = "partition")    # 按partition着色

# Partition很重要！
# 轨迹只在同一个partition内学习
# 如果你想分析的细胞在不同partition，需要设 use_partition=FALSE

2.4 学习轨迹图¶

# 构建轨迹（核心步骤）
cds <- learn_graph(
  cds,
  use_partition = TRUE,    # 每个partition独立学习轨迹
  close_loop = FALSE       # 是否允许环形轨迹
)

# 可视化轨迹
plot_cells(cds,
           color_cells_by = "cell_type",   # 按细胞类型着色
           label_groups_by_cluster = FALSE, # 不按cluster标注
           label_leaves = FALSE,            # 不标注叶节点
           label_branch_points = TRUE,      # 标注分支点
           graph_label_size = 3)            # 图标签大小

2.5 排序细胞（计算Pseudotime）¶

# 选择根节点（交互式或程序化）

# 方法1：交互式选择（弹出图形界面）
cds <- order_cells(cds)  # 在弹出的图上点击起始点

# 方法2：程序化选择（推荐，可重复）
# 找到特定细胞类型的节点作为根
get_earliest_principal_node <- function(cds, cell_type = "HSC") {
  cell_ids <- which(colData(cds)$cell_type == cell_type)  # 找到干细胞
  closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[cell_ids, ])  # 最近的图节点
  root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name  # 所有图节点
  root_pr_nodes[which.min(as.numeric(table(closest_vertex)[root_pr_nodes]))]  # 选择根
}

cds <- order_cells(cds, root_pr_nodes = get_earliest_principal_node(cds, "HSC"))

# 可视化Pseudotime
plot_cells(cds,
           color_cells_by = "pseudotime",   # 按pseudotime着色
           label_cell_groups = FALSE,
           label_leaves = FALSE,
           label_branch_points = FALSE,
           cell_size = 0.5) +
  scale_color_viridis_c()  # 用viridis配色

2.6 沿轨迹的差异基因¶

# 找沿pseudotime变化的基因
track_genes_result <- graph_test(
  cds,
  neighbor_graph = "principal_graph",  # 使用轨迹图
  cores = 4                             # 并行线程数
)

# 筛选显著基因
sig_genes <- track_genes_result %>%
  filter(q_value < 0.05) %>%           # FDR < 0.05
  arrange(q_value) %>%                  # 按q值排序
  head(100)                             # 取前100个

# 查看特定基因的表达随pseudotime变化
plot_cells(cds,
           genes = c("GATA1", "SPI1", "KLF1"),  # 指定基因
           show_trajectory_graph = TRUE,
           label_cell_groups = FALSE,
           cell_size = 0.5)

# 基因表达沿pseudotime的折线图
plot_genes_in_pseudotime(
  cds[rownames(cds) %in% c("GATA1", "SPI1", "KLF1"), ],  # 选择基因
  color_cells_by = "cell_type",  # 按细胞类型着色
  min_expr = 0.5                  # 最小表达阈值
)

2.7 分支分析¶

# 在分支点找差异表达基因
# 识别分支决策相关基因

# 获取分支信息
branch_genes <- graph_test(cds, neighbor_graph = "principal_graph")

# 用模块分析（基因分组）
gene_modules <- find_gene_modules(
  cds[rowData(cds)$gene_short_name %in% sig_genes$gene_short_name, ],
  resolution = 1e-2    # 模块分辨率
)

# 可视化基因模块
plot_cells(cds,
           genes = gene_modules,
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)

三、将Pseudotime整合回Seurat¶

# 提取pseudotime值
pseudotime_values <- pseudotime(cds)  # 获取每个细胞的pseudotime

# 添加到Seurat对象
seurat_obj$pseudotime <- pseudotime_values[colnames(seurat_obj)]

# 在Seurat中可视化
FeaturePlot(seurat_obj, features = "pseudotime") +
  scale_color_viridis_c()

# 按pseudotime排序做差异分析
VlnPlot(seurat_obj, features = c("GATA1", "SPI1"),
        group.by = "cell_type", pt.size = 0)

常见报错与解决¶

报错信息	原因	解决方案
`No root nodes specified`	没有选择轨迹起点	用order_cells()选择根节点
`Cells in different partitions`	要分析的细胞在不同分区	设use_partition=FALSE
`learn_graph: insufficient cells`	细胞太少	至少需要几百个细胞
`as.cell_data_set: missing metadata`	Seurat对象缺少信息	检查meta.data是否完整
`pseudotime is Inf`	某些细胞没有被连接到轨迹	检查分区和轨迹图连通性
`graph_test: slow`	基因太多太慢	只分析高变基因子集

速查表¶

# Monocle3标准流程
Seurat → as.cell_data_set() → preprocess_cds()
→ reduce_dimension() → cluster_cells() → learn_graph()
→ order_cells() → plot_cells(color_cells_by="pseudotime")
→ graph_test() → plot_genes_in_pseudotime()

# 核心函数
preprocess_cds()     — 预处理（PCA）
reduce_dimension()   — UMAP降维
cluster_cells()      — 聚类
learn_graph()        — 学习轨迹
order_cells()        — 计算pseudotime
graph_test()         — 轨迹差异基因
plot_cells()         — 可视化

# 参数选择
num_dim: 通常30-50（看碎石图）
resolution: 越小cluster越少
use_partition: TRUE（默认）或FALSE
close_loop: FALSE（大多数情况）

# 根节点选择策略
已知起点: 选干细胞/祖细胞所在节点
未知起点: 选最未分化的cluster
多个起点: 可以运行多次比较

面试高频问题¶

Q1：什么是Pseudotime？和真实时间有什么关系？¶

答：Pseudotime是一个抽象的"发育进度"度量，表示细胞在发育/分化过程中走了多远。它不等于真实时间——一个pseudotime=5的细胞不一定比pseudotime=3的细胞晚5-3=2天。Pseudotime衡量的是从起点到该细胞的最短路径上，转录组变化的总量。它假设处于不同分化阶段的细胞在同一时间点被捕获，通过表达模式相似性推断先后顺序。

Q2：Monocle3的轨迹推断原理是什么？¶

答：Monocle3先用UMAP降维，然后在降维空间中用逆图嵌入（reversed graph embedding）学习一个主图（principal graph）。这个图是细胞分布的"骨架"——连接的是细胞群的中心而非单个细胞。然后把每个细胞映射到图上最近的点，从用户指定的根节点出发，沿图的最短路径计算每个细胞的pseudotime。分支对应细胞命运的分叉点。

Q3：如何选择轨迹的起点（根节点）？¶

答：需要先验生物学知识。常见策略：①已知干细胞/祖细胞标记（如造血分化中的CD34+HSC）——选这些细胞所在的图节点作为根；②如果不确定，可以看哪些cluster表达已知的"未分化"标记基因；③也可以用RNA velocity作为辅助——velocity指向分化方向，起点应该是velocity的源头。根节点的选择直接影响pseudotime方向和结果解读。

Q4：轨迹分析有什么局限性？¶

答：①需要中间状态存在——如果数据中只有起点和终点细胞没有过渡态，轨迹推断会失败；②时间分辨率有限——pseudotime不等于真实时间；③分支不一定是真实的——噪声和批次效应可能产生虚假分支；④只能分析连续过程——不适合分析不连续的细胞状态变化（如免疫细胞激活）；⑤对参数敏感——不同参数可能产生不同轨迹。

Q5：Monocle3和其他轨迹推断工具有什么区别？¶

答：Monocle3的特点是使用主图（principal graph）学习轨迹，能处理复杂的分支和环形结构。与之对比：①Slingshot——基于主曲线（principal curves），更简单但不支持环形轨迹；②RNA velocity（scVelo）——利用剪接/未剪接RNA比例推断方向性，可以和Monocle3互补；③PAGA——基于图抽象化，适合探索性分析和粗粒度轨迹。实际分析中建议多工具验证。