ATAC + RNA 多组学：基于 ATAC 的 Monocle3 分析

作者: SeekGene

时长: 45 分钟

字数: 11.1k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于单细胞 ATAC-seq（scATAC-seq）数据，利用 Monocle3 轨迹推断算法，重建细胞从"起点"到"终点"的连续发育或分化过程。拟时序（Pseudotime）分析的核心思想是：将高维空间中离散的细胞，沿着一条隐含的生物学进程轨迹进行排序，从而揭示基因调控、染色质开放和转录因子活性在这一过程中的动态变化规律。

本教程围绕三个核心模态展开，系统性地重建"调控发起 → 染色质变化 → 基因表达"的完整逻辑链条：

Peak（染色质开放峰）：代表基因组中的顺式调控元件（如增强子、启动子），其开放程度随拟时序发生动态变化。
Motif（转录因子结合基序）：富集于开放 Peak 内的 DNA 序列模式，反映 TF 在染色质水平的潜在结合活性。
Gene（基因表达）：调控事件在转录层面的最终功能输出，沿拟时序呈现时空特异性的表达模式。

分析流程概览：

轨迹推断：基于降维嵌入（UMAP），利用 Monocle3 学习细胞状态之间的主图结构（Principal Graph），并以指定的根细胞群为起点，为每个细胞计算拟时序值。
差异开放峰检测：利用广义线性模型（GLM），识别沿拟时序显著变化的 Peak。
差异Motif 活性（chromVAR）检测：计算每个细胞中各 TF Motif 的染色质可及性偏差值，并通过线性模型评估其随拟时序的显著性变化。
差异表达基因检测（graph_test）：基于 Monocle3 主图拓扑结构，检测沿拟时序呈现空间自相关的基因。
多模态联合可视化：将 RNA、Motif、Peak 沿拟时序的动态变化整合到同一热图框架中，实现跨模态信号传递的全局展示。

💡 说明：文档中展示的图表仅为部分代表性的可视化结果。如需查看所有细胞群或对比组的完整分析结果（包含详细的数据表格与所有高清绘图），请前往 ./output/ 目录进行查阅。

suppressPackageStartupMessages({
    library(Signac)
    library(Seurat)
    library(SeuratWrappers)
    library(SummarizedExperiment)
    library(monocle3)
    library(cicero)
    library(Matrix)
    library(ggplot2)
    library(patchwork)
    library(parallel)
    library(Biostrings)
    library(TFBSTools)
    library(JASPAR2020)
    library(GenomeInfoDb)
    library(ComplexHeatmap)
    library(circlize)
})

2. 输入文件准备

2.1 输入文件准备与参数配置

在运行分析之前，请确保以下输入文件和参数已正确配置。脚本顶部提供了统一的参数配置区，用户可根据实际项目需求进行修改：

input.rds：预处理好的单细胞 ATAC-seq Seurat 对象（必须包含 ATAC Assay、聚类信息及 UMAP 降维结果）。
meta.tsv（可选）：附加的细胞元数据文件，包含 barcodes 列及额外的注释列。
genome.fa：对应物种的参考基因组 FASTA 文件，用于 chromVAR Motif 活性分析中的背景校正。

💡 说明此处参数需按照自身数据实际情况填写，尤其是文件路径、物种信息和目标细胞类型。所有路径建议使用绝对路径以避免工作目录混淆。

## outdir：输出目录路径，所有分析结果（图表、表格）将保存至此目录下的 output 子文件夹中
outdir = "./"

## ref_genome：参考基因组 FASTA 文件路径，用于 chromVAR Motif 活性分析中的 GC 含量背景校正
ref_genome = "/path/to/genome.fa"

## rds：输入的 Seurat 对象 RDS 文件路径，需包含 ATAC 和 RNA 两个 Assay 及预处理结果
rds = "/path/to/input.rds"

## meta：（可选）细胞元数据 TSV 文件路径，需包含 barcodes 列用于与 Seurat 对象匹配
meta = "/path/to/meta.tsv"

## Reduction：用于轨迹推断的降维嵌入方式，可选 "UMAP"、"ATACUMAP"、"WNNUMAP"（不区分大小写）
Reduction = "WNNUMAP"

## species：物种信息，"human" 或 "mouse" 等，用于自动选择对应的 JASPAR Motif 数据库
species = "human"

## root：指定作为拟时序起点的细胞类型名称，需存在于 clusters_col 对应的列中
root = "CMP"

## clusters_col：Seurat 对象 meta.data 中存储细胞类型/簇信息的列名
clusters_col = "Celltype"

## celltypes：待分析的细胞类型列表，用逗号分隔，仅保留这些类型的细胞用于分析
celltypes = "Monocytes,Pro B cells,CMP"

## sample_col：Seurat 对象 meta.data 中存储样本来源信息的列名
sample_col = "Sample"

## samples：待分析的样本 ID 列表，用逗号分隔，用于过滤特定样本的细胞
samples = "25030508_pbmc_1_arc,25030508_pbmc_1_arc"

## use_partition：是否使用 Monocle3 的分区（partition）信息进行轨迹学习，TRUE 或 FALSE
use_partition = "TRUE"

## downsample：是否对细胞进行下采样以减少计算量，TRUE 或 FALSE
downsample = "FALSE"

## downsample_num：当下采样启用时，每个细胞类型保留的最大细胞数
downsample_num = "1000"

💡 说明：首次分析时，不建议修改 use_partition 和 downsample 参数。use_partition = TRUE 会在细胞数较多时自动按连通域分区推断轨迹，有助于提高大规模数据的学习效率；downsample = FALSE 保留全部细胞以最大化轨迹分辨率。

⚠️ 注意：
root 参数指定的细胞类型必须在 celltypes 列表中存在，否则轨迹排序将失败。
Reduction 参数的名称（如 "WNNUMAP"）必须与 Seurat 对象中 obj@reductions 的键名完全匹配（本工具支持大小写不敏感匹配）。

2.2 辅助参数

配置绘图与分析中使用的辅助参数。

核心参数说明：

n_bins：将拟时序划分为多少区间（默认 10），用于 Peak 差异分析的细胞分箱。
top_n_heatmap：热图中展示的差异特征（RNA/Motif/Peak）的总数上限。

#不建议更换下面参数的值

n_bins <- 10
top_n_heatmap <- 100
my36colors <- c(
  "#E5D2DD", "#53A85F", "#F1BB72", "#F3B1A0", "#D6E7A3", "#57C3F3", "#476D87",
  "#E95C59", "#E59CC4", "#AB3282", "#23452F", "#BD956A", "#8C549C", "#585658",
  "#9FA3A8", "#E0D4CA", "#5F3D69", "#C5DEBA", "#58A4C3", "#E4C755", "#F7F398",
  "#AA9A59", "#E63863", "#E39A35", "#C1E6F3", "#6778AE", "#91D0BE", "#B53E2B",
  "#712820", "#DCC1DD", "#CCE0F5", "#CCC9E6", "#625D9E", "#68A180", "#3A6963",
  "#968175", "#6495ED", "#FFC1C1", "#f1ac9d", "#f06966", "#dee2d1", "#6abe83",
  "#39BAE8", "#B9EDF8", "#221a12", "#b8d00a", "#74828F", "#96C0CE", "#E95D22",
  "#017890"
)

2.3 输出目录初始化与参数解析

在进入正式分析前，创建输出目录并将字符串参数转换为 R 原生类型。

核心步骤解析：

目录创建：在主输出目录下创建 output/ 子目录，所有分析产物保存于此。
字符串拆分：将逗号分隔的 celltypes 和 samples 字符串拆分为字符向量。
类型转换：将 "TRUE"/"FALSE" 字符串转换为逻辑值，将 downsample_num 转换为数值。

⚠️ 注意：参数验证会在转换后立即执行，如果 downsample_num、use_partition 或 downsample 的类型转换失败（返回 NA），将抛出明确的错误信息。

# ---- Cell 4 ----
outdir <- paste0(outdir, "/output")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

celltypes <- strsplit(celltypes, ",")[[1]]
celltypes <- trimws(celltypes)
samples <- strsplit(samples, ",")[[1]]

use_partition <- as.logical(use_partition)
downsample <- as.logical(downsample)
downsample_num <- as.numeric(downsample_num)
if (is.na(downsample_num)) stop("downsample_num must be numeric.")
if (is.na(use_partition)) stop("use_partition must be TRUE/FALSE.")
if (is.na(downsample)) stop("downsample must be TRUE/FALSE.")

2.4 辅助工具函数

在进入正式分析之前，先定义一组将在后续流程中反复使用的辅助函数。这些函数负责名称规范化、降维映射解析、特征选择、矩阵平滑以及多种可视化类型的绘制。

💡 说明：以下函数为本教程内部实现，不涉及核心统计方法的修改。每个函数均有明确的输入/输出契约，可独立复用。

（1）名称规范化函数 用于降维嵌入名称的模糊匹配（大小写不敏感）。

normalize_name <- function(x) {
  gsub("[^a-z0-9]", "", tolower(x))
}

（2）降维嵌入解析函数

从 Seurat 对象中提取所有可用的降维结果名称。
将用户输入的 reduction_name 进行规范化（去除特殊字符、统一小写）。
在可用列表中进行匹配，找到对应的原始键名。

⚠️ 注意：如果指定的降维方法不存在，函数会列出所有可用的 reduction 名称以便调试。常见的可用选项包括 umap、ATACUMAP、WNNUMAP 等。

resolve_reduction <- function(obj, reduction_name) {
  available <- names(obj@reductions)
  if (length(available) == 0) {
    stop("No reductions found in Seurat object.")
  }
  normalized_available <- normalize_name(available)
  target <- normalize_name(reduction_name)
  idx <- match(target, normalized_available)
  if (is.na(idx)) {
    stop(
      paste0(
        "Reduction '", reduction_name, "' not found. Available reductions: ",
        paste(available, collapse = ", ")
      )
    )
  }
  available[idx]
}

（3）特征列自动选择函数 在不同分析步骤中，差异结果表格的列名可能为 site_name、id 或 feature_id。本函数自动识别并返回首个匹配的列名。

pick_feature_col <- function(df, candidates) {
  existing <- intersect(candidates, colnames(df))
  if (length(existing) == 0) {
    stop(paste0("Feature column not found. Candidates: ", paste(candidates, collapse = ", ")))
  }
  existing[1]
}

（4）多模态 CDS 切换函数

Monocle3 的 CellDataSet（CDS）对象默认绑定一种模态的计数矩阵（本教程中为 ATAC 可及性）。然而，在差异基因分析阶段，我们需要将相同的轨迹拓扑结构应用到 RNA 表达矩阵上，以检测基因表达沿拟时序的动态变化。

本函数在保持主图结构（Principal Graph）和降维嵌入不变的前提下，将 CDS 的计数矩阵从 ATAC 替换为 RNA，从而实现一套轨迹框架、多种模态分析。

核心步骤解析：

提取 Seurat 对象中指定 Assay（如 "RNA"）的归一化表达矩阵。
创建新的 CDS 对象，仅替换 expression_data 字段。
完整继承原 CDS 的 reducedDims、principal_graph 和 principal_graph_aux。

💡 说明：这种"模态切换"是 Monocle3 多组学分析的标准范式：轨迹用 ATAC 学习，基因表达沿同一轨迹进行 graph_test 检测。

swap_cds_assay <- function(main_cds, seurat_obj, assay) {
  common_cells <- intersect(colnames(main_cds), colnames(seurat_obj))
  main_cds <- main_cds[, common_cells]
  mat <- GetAssayData(seurat_obj, assay = assay, layer = "data")
  mat <- mat[, common_cells]
  
  new_cds <- new_cell_data_set(
    expression_data = mat,
    cell_metadata = colData(main_cds),
    gene_metadata = data.frame(gene_short_name = rownames(mat), row.names = rownames(mat))
  )
  reducedDims(new_cds) <- reducedDims(main_cds)
  new_cds@principal_graph_aux <- main_cds@principal_graph_aux
  new_cds@principal_graph <- main_cds@principal_graph
  return(new_cds)
}

（5）拟时序平滑矩阵构建函数

核心步骤解析：

拟时序排序：将所有有效细胞按拟时序值从小到大排序。
等距采样：在拟时序范围 [min(pt), max(pt)] 上均匀生成 n_points（默认 100）个采样点。
邻域均值平滑：对每个采样点，选取距离最近的 K 个细胞（至少 10 个，或总细胞数的 5%），计算其特征表达均值。这一步有效消除了单细胞数据的稀疏噪声（Drop-out）。
Z-score 标准化：对平滑后的矩阵按行（特征）进行 Z-score 标准化，使不同特征的相对变化趋势在同一尺度上可比较。

💡 说明：平滑矩阵是热图可视化的关键中间产物。它将离散的细胞表达数据转换为沿拟时序的连续信号曲线，使得 RNA、Motif、Peak 三个模态的热图能够精确对齐比较。

build_smoothed_matrix <- function(mat, pt, features, n_points = 100) {
  features <- intersect(features, rownames(mat))
  if (length(features) == 0) return(NULL)
  
  common_cells <- intersect(colnames(mat), names(pt))
  pt <- pt[common_cells]
  valid <- is.finite(pt)
  common_cells <- common_cells[valid]
  pt <- pt[valid]
  
  ord <- order(pt)
  common_cells <- common_cells[ord]
  pt <- pt[ord]
  m <- mat[features, common_cells, drop = FALSE]
  pt_seq <- seq(min(pt), max(pt), length.out = n_points)
  
  smooth_m <- matrix(0, nrow = length(features), ncol = n_points)
  rownames(smooth_m) <- features
  colnames(smooth_m) <- paste0("P", 1:n_points)
  
  for (i in 1:n_points) {
    dist <- abs(pt - pt_seq[i])
    idx <- order(dist)[1:max(10, floor(length(pt)*0.05))]
    smooth_m[, i] <- Matrix::rowMeans(m[, idx, drop=FALSE])
  }
  
  z <- t(scale(t(smooth_m)))
  z[is.na(z)] <- 0
  return(list(z = z, pseudotime = pt_seq))
}

（6）复杂热图绘制函数

本函数基于 ComplexHeatmap 绘制沿拟时序的特征动态变化热图。

核心步骤解析：

Motif 名称映射：将内部的 Motif ID（如 MA0001.1）映射为人类可读的 TF 名称（如 SP1）。
拟时序颜色条：在热图顶部添加拟时序梯度色条（深紫 → 青绿 → 亮黄），直观标注细胞在轨迹上的位置。
关键特征标注：使用 anno_mark 自动标注效应值最大的 Top 特征名称，用连线关联到对应热图行。

💡 说明：热图行按 K-means 分为 row_km = 2 组，通常分别代表"沿拟时序上调"和"沿拟时序下调"的特征群。

plot_complex_heatmap <- function(z_obj, title_text, top_labels, motif_names = NULL) {
  if (is.null(z_obj)) return(NULL)
  z_mat <- z_obj$z
  pt_vec <- z_obj$pseudotime
  
  if (!is.null(motif_names)) {
    rn <- rownames(z_mat)
    mapped <- motif_names[rn]
    rn <- ifelse(is.na(mapped), rn, mapped)
    rownames(z_mat) <- rn
    
    mapped_top <- motif_names[top_labels]
    top_labels <- ifelse(is.na(mapped_top), top_labels, mapped_top)
  }
  
  col_fun <- colorRamp2(c(-2, 0, 2), c("#2166AC", "#F7F7F7", "#B2182B"))
  pt_range <- range(pt_vec, na.rm = TRUE)
  if (!all(is.finite(pt_range)) || diff(pt_range) == 0) {
    pt_range <- c(0, 1)
  }
  pt_col_fun <- colorRamp2(
    c(pt_range[1], mean(pt_range), pt_range[2]),
    c("#440154", "#21908C", "#FDE725")
  )
  top_anno <- HeatmapAnnotation(
    Pseudotime = pt_vec,
    col = list(Pseudotime = pt_col_fun),
    show_annotation_name = TRUE
  )
  
  idx <- which(rownames(z_mat) %in% top_labels)
  row_anno <- NULL
  if (length(idx) > 0) {
    row_anno <- rowAnnotation(
      link = anno_mark(at = idx, labels = rownames(z_mat)[idx], labels_gp = gpar(fontsize = 8))
    )
  }
  
  ht <- Heatmap(
    z_mat,
    name = "Z-score",
    column_title = title_text,
    cluster_columns = FALSE,
    cluster_rows = TRUE,
    row_km = 2,
    col = col_fun,
    show_row_names = FALSE,
    show_column_names = FALSE,
    top_annotation = top_anno,
    right_annotation = row_anno,
    use_raster = TRUE,
    raster_quality = 3
  )
  return(ht)
}

（7）Top 特征散点图绘制函数

绘制 Top 特征（RNA 基因、Motif 或 Peak）沿拟时序的单细胞表达散点图，并叠加 LOESS 平滑曲线。

核心步骤解析：

数据对齐：将特征矩阵、拟时序向量和细胞类型注释对齐到共同的细胞集合。
分面可视化：使用 facet_wrap 为每个特征生成独立的子图，便于单独评估其动态模式。
LOESS 趋势线：在每个子图上叠加 LOESS 平滑曲线（黑色），直观展示特征随拟时序的整体变化方向。

💡 说明：
散点颜色：映射自细胞类型注释，帮助识别哪些细胞亚群贡献了该特征的主要信号。
LOESS 曲线：反映了该特征的总体趋势方向。曲线向上 = 沿拟时序上调；曲线向下 = 沿拟时序下调。

plot_top_feature_scatter <- function(
  mat,
  features,
  pseudotime_vec,
  cell_types,
  ct_colors,
  title_prefix,
  y_label,
  out_prefix,
  outdir
) {
  features <- unique(features)
  features <- features[features %in% rownames(mat)]
  if (length(features) == 0) return(invisible(NULL))

  common_cells <- intersect(colnames(mat), names(pseudotime_vec))
  common_cells <- intersect(common_cells, names(cell_types))
  pt <- pseudotime_vec[common_cells]
  valid <- is.finite(pt)
  common_cells <- common_cells[valid]
  pt <- as.numeric(pt[valid])
  ct <- as.character(cell_types[common_cells])

  plot_df_list <- lapply(features, function(feat) {
    y <- as.numeric(mat[feat, common_cells])
    data.frame(
      feature = feat,
      pseudotime = pt,
      value = y,
      celltype = ct,
      stringsAsFactors = FALSE
    )
  })
  plot_df <- do.call(rbind, plot_df_list)
  if (nrow(plot_df) == 0) return(invisible(NULL))

  plot_df$feature <- factor(plot_df$feature, levels = features)
  plot_df$celltype <- factor(plot_df$celltype, levels = names(ct_colors))

  p <- ggplot(plot_df, aes(x = pseudotime, y = value, color = celltype)) +
    geom_point(size = 0.6, alpha = 0.75) +
    geom_smooth(aes(group = feature), method = "loess", se = FALSE, color = "black", linewidth = 0.6) +
    facet_wrap(~feature, scales = "free_y", ncol = 2) +
    scale_color_manual(values = ct_colors, drop = FALSE) +
    labs(x = "pseudotime", y = y_label, title = title_prefix, color = "CellType") +
    theme_bw(base_size = 11) +
    theme(
      panel.grid = element_blank(),
      strip.background = element_rect(fill = "white"),
      plot.title = element_text(hjust = 0.5)
    )
  print(p) 
  ggsave(file.path(outdir, paste0(out_prefix, ".png")), plot = p, width = 10, height = 12, dpi = 300)
  ggsave(file.path(outdir, paste0(out_prefix, ".pdf")), plot = p, width = 10, height = 12)
}

3. 数据加载与预处理

3.1 读取 Seurat 对象

从 RDS 文件加载预处理好的单细胞 ATAC-seq 数据。如果提供了额外的元数据文件（meta.tsv），将其按 barcode 合并到 Seurat 对象中。

obj <- readRDS(rds)
if (meta != "") {
  meta_df <- read.table(meta, header = TRUE, sep = "\t", check.names = FALSE)
  rownames(meta_df) <- meta_df$barcodes
  obj <- AddMetaData(obj, meta_df)
}

3.2 细胞类型与样本筛选

仅保留 clusters_col 列中与 celltypes 列表匹配的细胞类型。
仅保留 sample_col 列中与 samples 列表匹配的样本。

obj@meta.data[[clusters_col]] <- as.character(obj@meta.data[[clusters_col]])
obj <- subset(obj, subset = !!sym(clusters_col) %in% celltypes)
obj <- subset(obj, subset = !!sym(sample_col) %in% samples)

3.3 下采样平衡

如果 downsample = TRUE，按 downsample_num 对每个聚类进行随机下采样，以减少计算负担。

DefaultAssay(obj) <- "ATAC"
Idents(obj) <- obj@meta.data[[clusters_col]]
if (downsample) {
  obj <- subset(obj, cells = WhichCells(obj, downsample = downsample_num))
}

⚠️ 注意：
过滤后如果细胞数过少（如 < 50），Monocle3 的图学习算法可能无法收敛。建议在运行前检查过滤后的细胞总数。
meta.tsv 文件必须包含 barcodes 列，且其值必须与 Seurat 对象中的细胞 barcode 一致。

4. Monocle3 轨迹推断

4.1 降维嵌入映射

Monocle3 内部使用 UMAP 作为降维坐标。由于 Seurat 对象中的降维结果可能命名为 WNNUMAP 或 ATACUMAP，我们需要将其提取并重新注册为 Monocle3 可识别的 "UMAP" reduction。

核心步骤解析：

降维名称解析：调用 resolve_reduction() 函数，将用户指定的 Reduction 名称匹配到 Seurat 对象中的实际键名（大小写不敏感）。
提取嵌入坐标：获取所有细胞在该降维下的坐标矩阵。
创建 Monocle3 兼容的 reduction：将提取的坐标注册为 "umap" reduction。

# ---- Cell 6 ----
selected_reduction <- resolve_reduction(obj, Reduction)
umap_embed <- Embeddings(obj, reduction = selected_reduction)
obj[["umap"]] <- CreateDimReducObject(
  embeddings = umap_embed,
  key = "UMAP_",
  assay = DefaultAssay(obj)
)

💡 说明： Monocle3 对 reduction 的键名有严格要求（必须包含 "UMAP"）。直接使用非标准名称会导致 learn_graph() 报错。

4.2 CDS 转换与主图学习

核心步骤解析：

Seurat → CDS 转换：利用 as.cell_data_set() 将 Seurat 对象转换为 Monocle3 的 CellDataSet 对象。这是 Monocle3 所有分析的数据结构基础。
UMAP 坐标注入：将之前提取的 UMAP 嵌入坐标直接写入 CDS 的 reducedDims，确保 Monocle3 使用与我们一致的降维结果。
基因检测：调用 detect_genes() 标记在每个细胞中表达的基因/特征，为后续差异分析提供背景集。
细胞聚类：基于 UMAP 坐标进行 cluster_cells()，识别细胞状态群。
主图学习（learn_graph）：Monocle3 的核心算法——在降维空间中学习连接各细胞状态群的图结构（Principal Graph）。use_partition 参数控制是否在大型数据集中使用分区模式。
拟时序排序（order_cells）：以指定的 root 细胞群为起点，沿主图计算每个细胞的拟时序值。
结果保存：将完整的 CDS 对象保存为 RDS 文件，并将拟时序值注入回 Seurat 对象的元数据中。

💡 说明：
根细胞选择：root 参数指定的细胞类型必须是生物学上合理的发育起点（如干细胞、祖细胞）。选择错误的根细胞会导致整个拟时序方向反转。
分区模式：当细胞存在多个不连通的发育分支时，use_partition = TRUE 会自动将数据分为独立分区分别学习。如果强制 use_partition = FALSE，Monocle3 会尝试将所有细胞连入同一张图，可能导致不合理的跨谱系连接。

# ---- 手动构建 CDS（兼容 Seurat v5）----

# 1. 获取 ATAC counts 数据（Seurat V5 用 layer 替代 slot）
counts_mat <- GetAssayData(obj, assay = "ATAC", layer = "counts")

# 2. 确保是标准的 dgCMatrix 格式
counts_mat <- as(counts_mat, "dgCMatrix")

# 3. 构建 CDS
cds <- new_cell_data_set(
  expression_data = counts_mat,
  cell_metadata = obj@meta.data,
  gene_metadata = data.frame(
    gene_short_name = rownames(counts_mat),
    row.names = rownames(counts_mat)
  )
)

# 4. 注入 UMAP
reducedDims(cds)$UMAP <- umap_embed[colnames(cds), , drop = FALSE]

cds <- detect_genes(cds)
cds <- cluster_cells(cds = cds, reduction_method = "UMAP")
cds <- learn_graph(cds, use_partition = use_partition)

root_cells <- rownames(obj@meta.data[obj@meta.data[[clusters_col]] == root, ])
if (length(root_cells) == 0) {
  stop("No root cells found. Check root and clusters_col.")
}
cds <- order_cells(cds, reduction_method = "UMAP", root_cells = root_cells)
saveRDS(cds, file.path(outdir, "cds.rds"))

pseudotime_vec <- pseudotime(cds)
obj <- AddMetaData(obj, metadata = pseudotime_vec, col.name = "Monocle3_pseudotime")

output

|======================================================================| 100%

4.3 轨迹可视化

将拟时序和细胞类型分别映射到 UMAP 上，生成四格图（带/不带轨迹图 × 两种着色方案）。

输出文件：

pseudotime.png/pdf：左图显示拟时序渐变 + 轨迹骨架，右图仅显示拟时序渐变。
celltype.png/pdf：左图显示细胞类型分布 + 轨迹骨架，右图仅显示细胞类型分布。

💡 说明：对比拟时序图和细胞类型图可以帮助您验证：轨迹方向是否与各细胞类型的预期发育顺序一致（如祖细胞 → 中间态 → 终末分化细胞）。

# ---- Cell 7 ----
p1 <- plot_cells(
  cds = cds,
  color_cells_by = "pseudotime",
  label_cell_groups = FALSE,
  show_trajectory_graph = TRUE
)

p2 <- plot_cells(
  cds = cds,
  color_cells_by = "pseudotime",
  label_cell_groups = FALSE,
  show_trajectory_graph = FALSE
)
p3 <- plot_cells(
  cds = cds,
  color_cells_by = clusters_col,
  show_trajectory_graph = TRUE,
  label_cell_groups = FALSE
)
p4 <- plot_cells(
  cds = cds,
  color_cells_by = clusters_col,
  show_trajectory_graph = FALSE,
  label_cell_groups = FALSE
)
ggsave(p1 + p2, file = file.path(outdir, "pseudotime.png"), width = 12, height = 5, dpi = 300)
ggsave(p3 + p4, file = file.path(outdir, "celltype.png"), width = 12, height = 5, dpi = 300)
ggsave(p1 + p2, file = file.path(outdir, "pseudotime.pdf"), width = 12, height = 5)
ggsave(p3 + p4, file = file.path(outdir, "celltype.pdf"), width = 12, height = 5)

options(repr.plot.width = 14, repr.plot.height = 7)
p1 + p2

options(repr.plot.width = 14, repr.plot.height = 7)
p3 + p4

5. 多组学动态特征检测（沿拟时序）

利用 Monocle3 的 fit_models() 函数，以广义线性模型（GLM）评估每个 Peak 的可及性是否随拟时序显著变化。

5.1 沿拟时序的染色质可及性差异分析

核心步骤解析：

过滤无效细胞：仅保留具有有限拟时序值的细胞（排除无法排序的细胞）。
拟时序分箱：将拟时序范围等分为 n_bins 个区间，每个区间视为一个"伪样本"。
细胞聚合：调用 aggregate_by_cell_bin() 将同一分箱内的细胞聚合为一个伪 bulk 样本，降低稀疏性。
GLM 拟合：使用公式 ~Pseudotime + num_genes_expressed 拟合每个 Peak 的广义线性模型。其中：
- Pseudotime：主效应项，评估 Peak 可及性随拟时序的变化。
- num_genes_expressed：协变量，校正测序深度差异。
显著性过滤：提取 term == "Pseudotime" 且 q_value < 0.05 的 Peak，保存为 CSV 文件。

💡 说明：
q_value：经过 Benjamini-Hochberg (BH) 多重检验校正后的 P-value。q_value < 0.05 意味着在 5% 的假发现率（FDR）水平下显著。
estimate 符号：正值表示 Peak 可及性随拟时序增加（开放增强）；负值表示随拟时序降低（逐渐关闭）。

⚠️ 注意： fit_models() 的计算量与 Peak 数量成正比。对于 > 50K Peak 的数据集，建议使用 cores 参数指定多线程加速。本教程自动使用系统核心数的 80%（上限 25 核）。

# ---- Cell 8: differential peaks along pseudotime ----
cores <- ifelse(detectCores() > 32, 25, trunc(0.8 * detectCores()))

input_cds_lin <- cds[, is.finite(pseudotime_vec[colnames(cds)])]
pData(input_cds_lin)$Pseudotime <- pseudotime(input_cds_lin)
pData(input_cds_lin)$cell_subtype <- cut(pseudotime(input_cds_lin), n_bins)
binned_input_lin <- aggregate_by_cell_bin(input_cds_lin, "cell_subtype")
acc_fits <- fit_models(
  binned_input_lin,
  cores = cores,
  model_formula_str = "~Pseudotime + num_genes_expressed"
)
fit_coefs <- coefficient_table(acc_fits)
significant_sites <- subset(fit_coefs, term == "Pseudotime" & q_value < 0.05)

important_cols <- intersect(
  c("site_name", "estimate", "normalized_effect", "p_value", "q_value"),
  colnames(significant_sites)
)
important_data <- significant_sites[, important_cols, drop = FALSE]
write.csv(
  important_data,
  file.path(outdir, "pseudotime_differentially_accessible_sites.csv"),
  row.names = FALSE
)
peak <- unique(important_data$site_name)

writeLines(
  as.character(peak),
  con = file.path(outdir, "peak.txt")
)

head(important_data)

A tibble: 6 × 5
site_name	estimate	normalized_effect	p_value	q_value
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
chr14-69042674-69044104	0.1474021	0.21239396	1.579794e-07	0.021606051
chr16-30496998-30497771	1.3722741	0.00000000	2.556358e-07	0.034961780
chr16-71972739-71973671	1.6343574	0.00000000	5.002789e-08	0.006842215
chr17-76604881-76606067	1.3722741	0.00000000	2.556358e-07	0.034961780
chr2-85963580-85964391	3.7951820	0.03617204	2.921967e-08	0.003996345
chr4-106266040-106266712	1.3586657	0.00000000	9.671732e-08	0.013227641

5.2 Top 变化 Peak 可视化

将显著变化的 Peak 按效应值（estimate）排序，分别提取最显著上调和下调的 Top 5 Peak，并绘制其在拟时序中的可及性变化曲线。

输出文件：

top_increasing.png/pdf：沿拟时序逐渐开放的 Top Peak。
top_decreasing.png/pdf：沿拟时序逐渐关闭的 Top Peak。

💡 说明：这些 Peak 对应的基因组区域可能包含关键的顺式调控元件（如发育特异的增强子）。结合下游的 Motif 分析，可以进一步揭示驱动这些变化的转录因子。


# ---- Cell 9: top changing peaks ----
# Plotting top increasing/decreasing point plots
peak_feature_col <- pick_feature_col(significant_sites, c("site_name", "id", "feature_id"))
positive_sites <- significant_sites[significant_sites$estimate > 0, , drop = FALSE]
negative_sites <- significant_sites[significant_sites$estimate < 0, , drop = FALSE]
top_increasing <- head(positive_sites[order(-positive_sites$estimate), , drop = FALSE], 5)
top_decreasing <- head(negative_sites[order(-abs(negative_sites$estimate)), , drop = FALSE], 5)
top_increasing <- top_increasing[!is.na(top_increasing[[peak_feature_col]]), ]
top_decreasing <- top_decreasing[!is.na(top_decreasing[[peak_feature_col]]), ]

if (nrow(top_increasing) > 0) {
  p5 <- plot_accessibility_in_pseudotime(input_cds_lin[top_increasing[[peak_feature_col]]])
  ggsave(p5, file = file.path(outdir, "top_increasing.png"), width = 5, height = 8, dpi = 300)
  ggsave(p5, file = file.path(outdir, "top_increasing.pdf"), width = 5, height = 8)
} else {
  message("No positive pseudotime peaks found after filtering; skip top_increasing plot.")
}

if (nrow(top_decreasing) > 0) {
  p6 <- plot_accessibility_in_pseudotime(input_cds_lin[top_decreasing[[peak_feature_col]]])
  ggsave(p6, file = file.path(outdir, "top_decreasing.png"), width = 5, height = 8, dpi = 300)
  ggsave(p6, file = file.path(outdir, "top_decreasing.pdf"), width = 5, height = 8)
} else {
  message("No negative pseudotime peaks found after filtering; skip top_decreasing plot.")
}

output

No negative pseudotime peaks found after filtering; skip top_decreasing plot.

p5

💡 解读指南
全局说明： top_increasing.png 展示了 5 个随拟时序最显著上升的可及性位点，用于观察"染色质逐步开放"的动态模式。
横轴（pseudotime）： 表示细胞沿轨迹从起始到终末的相对进程，越靠右通常越接近终末状态。
纵轴（accessibility）： 表示对应 peak 的可及性强度（原始值或标准化值）；数值越高代表该位点越开放。
趋势判断： 若曲线整体随 pseudotime 上升，说明该位点在发育推进过程中逐步激活，可能参与后期状态建立。
位点间比较： 上升更陡、末端水平更高的位点，往往提示更强的阶段特异开放信号，可优先进入下游验证。

#p6

💡 解读指南
全局说明： top_decreasing.png 展示了 5 个随拟时序最显著下降的可及性位点，反映"染色质逐步关闭"的动态过程。
横轴（pseudotime）： 表示细胞从起始到终末的轨迹进程，越靠右通常越接近终末状态。
纵轴（accessibility）： 表示对应 peak 的可及性水平；数值下降说明该区域逐步关闭。
趋势判断： 若曲线整体向下，提示该位点在分化推进过程中持续失活，可能对应早期程序的退出。
位点间比较： 下降更陡、末端更低的位点通常具有更强的阶段特异关闭特征，可作为重点候选。

5.3 沿拟时序的基因表达差异分析（RNA）

基因表达水平的变化是细胞状态转变的直接体现。为了识别随拟时序显著波动的基因，程序切换到 RNA 表达矩阵进行分析：

CDS 切换：通过 swap_cds_assay 函数将 CDS 对象的表达数据替换为 RNA Assay 的数据，同时保留原有的轨迹结构和降维信息，确保分析在同一拟时序框架下进行。
低表达基因过滤：剔除在少于10个细胞中表达的基因，减少稀疏噪声对分析的干扰。
空间自相关检验：利用 graph_test 函数，以主图（principal_graph）为邻居图结构，计算每个基因的 Moran's I 统计量。Moran's I 衡量基因表达在轨迹空间中的自相关性——高 Moran's I 值意味着基因表达沿轨迹呈现连续、平滑的变化模式，而非随机分布。
显著基因筛选：以 q_value < 0.05 且 morans_I > 0.05 为双重阈值，筛选既具有统计显著性又具有生物学意义的轨迹相关基因。

分析结果保存为 pseudotime_differentially_expressed_genes.csv，基因列表同时输出为 genes.txt，便于后续功能富集分析或与其他组学数据进行交叉验证。

# ---- Cell 11: differential genes along pseudotime ----
cds_rna <- swap_cds_assay(cds, obj, "RNA")
# Filter genes expressed in at least 10 cells to speed up graph_test
cds_rna <- cds_rna[Matrix::rowSums(SingleCellExperiment::counts(cds_rna) > 0) > 10, ]
rna_graph_res <- graph_test(cds_rna, neighbor_graph="principal_graph", cores=cores)
significant_genes <- subset(rna_graph_res, q_value < 0.05 & morans_I > 0.05)

write.csv(
  as.data.frame(significant_genes),
  file.path(outdir, "pseudotime_differentially_expressed_genes.csv"),
  row.names = FALSE
)

gene_list <- unique(as.character(significant_genes$gene_short_name))
gene_list <- gene_list[!is.na(gene_list) & nzchar(gene_list)]

writeLines(gene_list, con = file.path(outdir, "genes.txt"))

5.4 沿拟时序的 Motif 活性差异分析

转录因子通过结合特异的 DNA 基序（Motif）来调控基因表达，其结合活性的动态变化是理解上游调控逻辑的关键。本教程利用 chromVAR 推断的 Motif 活性矩阵进行差异分析：

活性矩阵提取：从 Seurat 对象的 chromvar Assay 中提取 Motif 活性得分矩阵，每个值代表该 Motif 在单个细胞中的相对结合活性。
线性回归建模：由于 chromVAR 活性数据不符合计数分布假设，程序采用经典线性回归模型，对每个 Motif 独立拟合 活性 ~ 拟时间 的关系，提取回归斜率（estimate）和显著性（p_value）。
多重假设检验校正：对所有 Motif 的 P 值进行 Benjamini-Hochberg 校正，以 q_value < 0.05 为阈值筛选活性随拟时序显著变化的 Motif。
TF 名称映射：将 Motif ID（如 MA0137.3）映射为易读的转录因子名称（如 STAT1），提升结果的可解释性。映射失败的 Motif 保留原始 ID。

分析结果保存为 pseudotime_differentially_active_motifs.csv，包含 Motif ID、TF 名称、效应估计值和显著性指标。转录因子列表同时输出为 TFs.txt，可直接用于后续的调控网络构建。

# ---- Cell 10: motif activity analysis ----
genome_seq <- readDNAStringSet(ref_genome)
names(genome_seq) <- sub(" .*", "", names(genome_seq))
keep_chromosomes <- unique(as.character(seqnames(granges(obj[["ATAC"]]))))
genome_seq <- genome_seq[names(genome_seq) %in% keep_chromosomes]
obj@assays$ATAC@ranges <- keepSeqlevels(
  obj@assays$ATAC@ranges,
  keep_chromosomes,
  pruning.mode = "coarse"
)

if (species == "human") {
  pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = 9606, all_versions = TRUE))
} else if (species == "mouse") {
  pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = 10090, all_versions = TRUE))
} else {
  stop("The value of 'species' is not supported.")
}

obj <- AddMotifs(object = obj, genome = genome_seq, pfm = pfm)
obj <- RunChromVAR(object = obj, genome = genome_seq)
saveRDS(obj, file.path(outdir, "output.rds"))

# ---- Cell 12: differential motif activity along pseudotime ----
DefaultAssay(obj) <- "chromvar"
motif_mat_tmp <- GetAssayData(obj, assay = "chromvar", layer = "data")
motif_cells <- intersect(colnames(motif_mat_tmp), names(pseudotime_vec))
motif_pt <- pseudotime_vec[motif_cells]
valid_pt <- is.finite(motif_pt)
motif_cells <- motif_cells[valid_pt]
motif_pt <- as.numeric(motif_pt[valid_pt])
motif_mat_tmp <- motif_mat_tmp[, motif_cells, drop = FALSE]

fit_single_motif <- function(y, x) {
  ok <- is.finite(y) & is.finite(x)
  if (sum(ok) < 10 || sd(y[ok]) == 0) {
    return(c(estimate = NA_real_, p_value = NA_real_, normalized_effect = NA_real_))
  }
  model <- lm(y[ok] ~ x[ok])
  est <- unname(coef(model)[2])
  pval <- summary(model)$coefficients[2, 4]
  norm_eff <- est / sd(y[ok])
  c(estimate = est, p_value = pval, normalized_effect = norm_eff)
}

motif_stats <- t(vapply(
  seq_len(nrow(motif_mat_tmp)),
  function(i) fit_single_motif(as.numeric(motif_mat_tmp[i, ]), motif_pt),
  numeric(3)
))

motif_coef <- data.frame(
  id = rownames(motif_mat_tmp),
  term = "Pseudotime",
  estimate = motif_stats[, "estimate"],
  p_value = motif_stats[, "p_value"],
  normalized_effect = motif_stats[, "normalized_effect"],
  stringsAsFactors = FALSE
)
motif_coef$q_value <- p.adjust(motif_coef$p_value, method = "BH")
motif_coef$status <- ifelse(is.na(motif_coef$p_value), "FAIL", "OK")




significant_motifs <- subset(motif_coef, status == "OK" & q_value < 0.05)

if (!exists("motif_name_map") || is.null(motif_name_map)) {
  motif_obj <- Motifs(obj[["ATAC"]])
  motif_name_map <- if (!is.null(motif_obj)) unlist(motif_obj@motif.names) else NULL
}

if (!is.null(motif_name_map)) {
  mapped_tf <- unname(motif_name_map[significant_motifs$id])
  significant_motifs$TF_name <- ifelse(is.na(mapped_tf), significant_motifs$id, mapped_tf)
} else {
  significant_motifs$TF_name <- significant_motifs$id
}

write.csv(
  significant_motifs,
  file.path(outdir, "pseudotime_differentially_active_motifs.csv"),
  row.names = FALSE
)

tf_list <- unique(as.character(significant_motifs$TF_name))
tf_list <- tf_list[!is.na(tf_list) & nzchar(tf_list)]
writeLines(tf_list, con = file.path(outdir, "TFs.txt"))

6. 多组学动态特征可视化

在完成三类分子特征的差异分析之后，核心任务是将这些动态变化以直观的形式呈现出来。本教程从两个维度进行可视化：平滑热图展示全局动态模式，散点图展示单个特征的连续变化趋势。

6.1 平滑热图数据构建

为了在拟时序轴上连续展示大量特征的动态变化，程序首先构建平滑表达矩阵：

Motif ID 映射：从 Seurat 对象的 ATAC Assay 中提取 Motif 名称映射表（motif_name_map），将 MA ID（如 MA0137.3）转换为易读的转录因子名称（如 STAT1），提升热图的可解释性。
Top 标签选取：分别从三类显著特征中，按效应量绝对值排序选取 Top 50（top_n_heatmap / 2）个最具代表性的特征作为热图的行标签：
- RNA：按 Moran's I 降序排列，选取表达沿轨迹空间自相关性最强的基因。
- Motif：按线性回归估计值绝对值降序排列，选取活性变化幅度最大的 Motif。
- Peak：按 fit_models 回归估计值绝对值降序排列，选取开放性变化幅度最大的 Peak。
全量特征矩阵构建：使用所有显著特征（而非仅 Top 标签）作为热图的行数据，确保热图展示完整的动态全貌。
平滑矩阵计算：利用 build_smoothed_matrix 函数，沿拟时序方向对每个特征的表达/活性值进行局部加权平均平滑处理（100 个等距点），生成 Z-score 标准化矩阵，有效降低单细胞噪声。
细胞类型颜色映射：为所有涉及的细胞类型分配统一颜色方案，确保跨图表的一致性。

# Map motif IDs to names
motif_obj <- Motifs(obj[["ATAC"]])
if (!is.null(motif_obj)) {
  motif_name_map <- unlist(motif_obj@motif.names)
} else {
  motif_name_map <- NULL
}

# Select top labels based on absolute effect size or morans_I
top_rna_labels <- rownames(head(significant_genes[order(-significant_genes$morans_I), ], top_n_heatmap / 2))
# Motif uses LM estimate, order by absolute estimate
top_motif_labels <- significant_motifs$id[order(-abs(significant_motifs$estimate))][1:(top_n_heatmap / 2)]
top_motif_labels <- top_motif_labels[!is.na(top_motif_labels)]
# Peak uses fit_models estimate
peak_feature_col <- pick_feature_col(significant_sites, c("site_name", "id", "feature_id"))
top_atac_labels <- significant_sites[[peak_feature_col]][order(-abs(significant_sites$estimate))][1:(top_n_heatmap / 2)]
top_atac_labels <- top_atac_labels[!is.na(top_atac_labels)]

# Use ALL significant features for heatmap rows
all_sig_genes <- rownames(significant_genes)
all_sig_motifs <- significant_motifs$id
all_sig_peaks <- significant_sites[[peak_feature_col]]

rna_mat <- GetAssayData(obj, assay = "RNA", layer = "data")
motif_mat <- GetAssayData(obj, assay = "chromvar", layer = "data")
atac_mat <- GetAssayData(obj, assay = "ATAC", layer = "data")
cell_types <- obj@meta.data[[clusters_col]]
names(cell_types) <- colnames(obj)

z_rna_obj <- build_smoothed_matrix(rna_mat, pseudotime_vec, all_sig_genes, n_points = 100)
z_motif_obj <- build_smoothed_matrix(motif_mat, pseudotime_vec, all_sig_motifs, n_points = 100)
z_peak_obj <- build_smoothed_matrix(atac_mat, pseudotime_vec, all_sig_peaks, n_points = 100)

6.2 热图绘制与输出

程序利用 plot_complex_heatmap 函数生成三类热图对象，再通过 safe_draw_heatmap 函数安全地输出为文件：

plot_complex_heatmap 函数：构建带有拟时序颜色条注释和特征标签的 ComplexHeatmap 对象。热图行自动聚为 2 个簇，分别对应随拟时序上升和下降的两类动态模式。右侧用 anno_mark 标注 Top 50 特征名称，便于识别关键分子。
safe_draw_heatmap 函数：封装了 tryCatch 异常处理机制，确保热图以 PNG（300 DPI 高分辨率位图）和 PDF（矢量图）两种格式稳定输出。即使某一种格式生成失败，也不会中断整个分析流程。

最终生成三个热图文件：

trajectory_rna_heatmap.png/pdf：基因表达沿拟时序的动态变化。
trajectory_motif_heatmap.png/pdf：Motif 染色质可及性沿拟时序的动态变化。
trajectory_peak_heatmap.png/pdf：Peak 开放性沿拟时序的动态变化。

# Generate consistent colors for cell types
unique_cell_types <- sort(unique(as.character(cell_types)))
unique_cell_types <- unique_cell_types[!is.na(unique_cell_types) & nzchar(unique_cell_types)]
n_colors <- length(unique_cell_types)
if (n_colors == 0) {
  stop("No valid cell types found for heatmap annotation.")
}
ct_colors <- rep(my36colors, length.out = n_colors)
names(ct_colors) <- unique_cell_types

ht_rna <- plot_complex_heatmap(z_rna_obj, "RNA expression", top_rna_labels)
ht_motif <- plot_complex_heatmap(z_motif_obj, "Motif chromatin accessibility", top_motif_labels, motif_name_map)
ht_peak <- plot_complex_heatmap(z_peak_obj, "Peak accessibility", top_atac_labels)

# 封装一个通用的安全绘图函数，避免重复写 tryCatch
safe_draw_heatmap <- function(ht_obj, file_prefix, width = 12, height = 8, outdir = outdir) {
  # 尝试画 PDF
  pdf_file <- file.path(outdir, paste0(file_prefix, ".pdf"))
  tryCatch({
    pdf(pdf_file, width = width, height = height)
    draw(ht_obj, ht_gap = unit(1, "cm"))
    dev.off()
    print(draw(ht_obj, ht_gap = unit(1, "cm")))  
    #message("Successfully generated: ", pdf_file)
  }, error = function(e) {
    if (names(dev.cur()) != "null device") dev.off() # 确保关闭设备，防止后续绘图错乱
    warning("Failed to generate PDF for ", file_prefix, ". Reason: ", e$message)
  })
}
# RNA Panel
if (exists("ht_rna")) {
  safe_draw_heatmap(ht_rna, "trajectory_rna_heatmap", outdir = outdir)
} else {
  warning("ht_rna object does not exist. Skipping RNA heatmap.")
}

# Motif Panel
if (exists("ht_motif")) {
  safe_draw_heatmap(ht_motif, "trajectory_motif_heatmap", outdir = outdir)
} else {
  warning("ht_motif object does not exist. Skipping Motif heatmap.")
}

# Peak Panel
if (exists("ht_peak")) {
  safe_draw_heatmap(ht_peak, "trajectory_peak_heatmap", outdir = outdir)
} else {
  warning("ht_peak object does not exist. Skipping Peak heatmap.")
}

💡 解读指南
全局说明： 热图展示了所有显著特征沿拟时序的 Z-score 标准化表达/活性变化，行代表特征，列代表拟时序进程（从左到右）。
颜色映射： 红色表示高表达/高活性，蓝色表示低表达/低活性，白色表示中等水平。
顶部注释条： 展示拟时序从起始（紫色）到终末（黄色）的渐变过程。
行聚类： 热图自动将特征分为两个簇——上方簇通常随拟时序上升（晚期激活），下方簇随拟时序下降（早期激活后退出的程序）。
右侧标签： 标注了变化幅度最大的 Top 50 特征名称，便于快速识别关键调控因子。

6.3 Top 10 特征散点图

为了更细致地观察单个特征的动态变化模式，程序进一步提取每类特征中变化最显著的 Top 10，绘制沿拟时序的散点图：

特征选取：
- RNA：按 Moran's I 降序取 Top 10 基因。
- Motif：按回归估计值绝对值降序取 Top 10 Motif。
- Peak：按回归估计值绝对值降序取 Top 10 Peak。
散点图绘制：利用 plot_top_feature_scatter 函数，以拟时间为横轴、特征值为纵轴，每个点代表一个细胞，按细胞类型着色。叠加黑色 LOESS 平滑曲线展示整体趋势。
分面展示：每个特征单独一个小面板，2 列排布，便于逐一查看。
结果输出：三类特征的散点图分别保存为：
- top10_rna_scatter_pseudotime.png/pdf
- top10_motif_scatter_pseudotime.png/pdf
- top10_peak_scatter_pseudotime.png/pdf


# ---- Cell 13.5: top10 feature scatter plots along pseudotime ----
top10_rna <- rownames(head(significant_genes[order(-significant_genes$morans_I), ], 10))
top10_motif <- significant_motifs$id[order(-abs(significant_motifs$estimate))][1:min(10, nrow(significant_motifs))]
top10_peak <- significant_sites[[peak_feature_col]][order(-abs(significant_sites$estimate))][1:min(10, nrow(significant_sites))]

plot_top_feature_scatter(
  mat = rna_mat,
  features = top10_rna,
  pseudotime_vec = pseudotime_vec,
  cell_types = cell_types,
  ct_colors = ct_colors,
  title_prefix = "Top10 RNA Features Along Pseudotime",
  y_label = "Gene expression",
  out_prefix = "top10_rna_scatter_pseudotime",
  outdir = outdir
)

output

\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'

plot_top_feature_scatter(
  mat = motif_mat,
  features = top10_motif,
  pseudotime_vec = pseudotime_vec,
  cell_types = cell_types,
  ct_colors = ct_colors,
  title_prefix = "Top10 Motif Features Along Pseudotime",
  y_label = "Motif activity",
  out_prefix = "top10_motif_scatter_pseudotime",
  outdir = outdir
)

output

\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'

plot_top_feature_scatter(
  mat = atac_mat,
  features = top10_peak,
  pseudotime_vec = pseudotime_vec,
  cell_types = cell_types,
  ct_colors = ct_colors,
  title_prefix = "Top10 Peak Features Along Pseudotime",
  y_label = "Peak accessibility",
  out_prefix = "top10_peak_scatter_pseudotime",
  outdir = outdir
)

output

\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'
\`geom_smooth()\` using formula = 'y ~ x'

💡 解读指南
全局说明： 散点图展示了 Top 10 特征沿拟时序的连续变化，每个面板对应一个特征。
横轴（pseudotime）： 表示细胞从起始到终末的轨迹进程，越靠右越接近终末状态。
纵轴： 对应特征的表达/活性/开放性值。
散点颜色： 按细胞类型着色，便于观察不同细胞类型在轨迹上的分布位置。
黑色平滑曲线： LOESS 拟合的趋势线，直观展示特征的整体变化方向——上升、下降或先升后降。
应用场景： 通过观察曲线的形状和陡峭程度，可以判断该特征是早期激活、晚期激活还是瞬时表达，为后续功能验证提供线索。

6.4 TF 特异性分析：聚焦转录因子动态

在完成基因水平的差异表达分析之后，为进一步聚焦上游调控因子，程序专门从显著差异表达的基因中筛选转录因子（TF），并进行独立的可视化分析。

TF 名称列表构建：从 Seurat 对象的 Motif 注释中提取所有已知的 TF 名称。对于异源二聚体（如 FOS::JUN），按 :: 拆分为独立的 TF 名称（FOS 和 JUN），并去除版本标签（如 (var.2)）和首尾空格，确保名称干净、匹配准确。
TF 筛选：将显著差异表达基因的基因名与 TF 名称列表进行大小写不敏感匹配，筛选出其中属于转录因子的子集（significant_tfs）。这些 TF 代表了在细胞分化过程中表达水平发生显著变化的潜在上游调控因子。
结果导出：将筛选出的 TF 列表保存为 pseudotime_differentially_expressed_TFs.csv，包含每个 TF 的 Moran's I、P 值、q 值等统计信息，便于后续功能验证或文献检索。
TF 表达热图：利用与 RNA 热图相同的平滑矩阵构建流程，专门为所有显著 TF 绘制沿拟时序的表达热图。热图以 trajectory_tf_expression_heatmap.png/pdf 命名输出，行标签标注变化最显著的 Top 50 TF 名称，直观展示转录因子在分化进程中的表达动态。

if (!is.null(motif_obj)) {
  motif_names_raw <- unlist(motif_obj@motif.names)
  # Split heterodimers like FOS::JUN and remove version tags like (var.2)
  tf_names_clean <- unique(unlist(strsplit(motif_names_raw, "::")))
  tf_names_clean <- gsub("\\(var\\.[0-9]+\\)", "", tf_names_clean)
  tf_names_clean <- trimws(tf_names_clean)
  
  # Filter significant genes to only include TFs
  rna_genes <- rownames(significant_genes)
  is_tf <- toupper(rna_genes) %in% toupper(tf_names_clean)
  significant_tfs <- significant_genes[is_tf, , drop = FALSE]
  
  write.csv(
    as.data.frame(significant_tfs),
    file.path(outdir, "pseudotime_differentially_expressed_TFs.csv"),
    row.names = FALSE
  )
      
  # Plot TF heatmap
  if (nrow(significant_tfs) > 0) {
    all_sig_tfs <- rownames(significant_tfs)
    # Select top TFs for labeling
    top_tf_labels <- rownames(head(significant_tfs[order(-significant_tfs$morans_I), ], top_n_heatmap / 2))
    
    z_tf_obj <- build_smoothed_matrix(rna_mat, pseudotime_vec, all_sig_tfs, n_points = 100)
    
    if (!is.null(z_tf_obj)) {
      ht_tf <- plot_complex_heatmap(z_tf_obj, "TF RNA expression", top_tf_labels)
      pdf(file.path(outdir, "trajectory_tf_expression_heatmap.pdf"), width = 8, height = 8)
      draw(ht_tf)
      dev.off()
      draw(ht_tf)        
    }
  }
}

💡 解读指南
全局说明： TF 特异性热图聚焦于转录因子这一特殊类别，帮助研究者快速定位可能驱动细胞命运决定的上游调控因子。
筛选逻辑： 只有同时满足两个条件的基因才会出现在此热图中：（1）在 RNA 水平上沿拟时序显著差异表达；（2）基因名存在于已知的转录因子数据库中。
生物学意义： 如果一个 TF 的表达沿拟时序显著变化，提示它可能在分化过程中发挥阶段性调控作用——早期高表达的 TF 可能启动分化程序，晚期高表达的 TF 可能维持终末状态。
与 Motif 热图的区别： Motif 热图反映的是染色质层面的 TF 结合活性（由 chromVAR 推断），而 TF 热图反映的是 RNA 层面的 TF 自身表达水平。两者结合可以判断：TF 的表达变化是否与其结合活性的变化一致，从而推断调控是发生在转录水平还是翻译后水平。

ATAC + RNA 多组学：基于 ATAC 的 Monocle3 分析 ​

1. 教程简介 ​

2. 输入文件准备 ​

2.1 输入文件准备与参数配置 ​

2.2 辅助参数 ​

2.3 输出目录初始化与参数解析 ​

2.4 辅助工具函数 ​

3. 数据加载与预处理 ​

3.1 读取 Seurat 对象 ​

3.2 细胞类型与样本筛选 ​

3.3 下采样平衡 ​

4. Monocle3 轨迹推断 ​

4.1 降维嵌入映射 ​

4.2 CDS 转换与主图学习 ​

4.3 轨迹可视化 ​

5. 多组学动态特征检测（沿拟时序） ​

5.1 沿拟时序的染色质可及性差异分析 ​

5.2 Top 变化 Peak 可视化 ​

5.3 沿拟时序的基因表达差异分析（RNA） ​

5.4 沿拟时序的 Motif 活性差异分析 ​

6. 多组学动态特征可视化 ​

6.1 平滑热图数据构建 ​

6.2 热图绘制与输出 ​

6.3 Top 10 特征散点图 ​

6.4 TF 特异性分析：聚焦转录因子动态 ​

ATAC + RNA 多组学：基于 ATAC 的 Monocle3 分析

1. 教程简介

2. 输入文件准备

2.1 输入文件准备与参数配置

2.2 辅助参数

2.3 输出目录初始化与参数解析

2.4 辅助工具函数

3. 数据加载与预处理

3.1 读取 Seurat 对象

3.2 细胞类型与样本筛选

3.3 下采样平衡

4. Monocle3 轨迹推断

4.1 降维嵌入映射

4.2 CDS 转换与主图学习

4.3 轨迹可视化

5. 多组学动态特征检测（沿拟时序）

5.1 沿拟时序的染色质可及性差异分析

5.2 Top 变化 Peak 可视化

5.3 沿拟时序的基因表达差异分析（RNA）

5.4 沿拟时序的 Motif 活性差异分析

6. 多组学动态特征可视化

6.1 平滑热图数据构建

6.2 热图绘制与输出

6.3 Top 10 特征散点图

6.4 TF 特异性分析：聚焦转录因子动态