ATAC + RNA 多组学：DORC 分析

作者: SeekGene

时长: 29 分钟

字数: 7.0k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于 Seurat 和 Signac 等单细胞多组学分析工具，专门用于对 SeekArc 单细胞多组学数据（scRNA-seq + scATAC-seq） 进行深度的 DORC (Domains of Regulatory Chromatin, 调控染色质结构域) 分析。该教程面向已完成基础质控、细胞群注释及初步 Peak Calling 的 SeekArc 多组学数据，通过整合染色质可及性与基因表达水平，实现以下三大核心分析目标：

• Peak-Gene 关联网络构建 (Peak-to-Gene Linkage)：基于共变性假设，计算每个基因与其附近（通常在数百 kb 范围内）各个 Peak 的相关性，识别具有显著顺式调控关系的 Peak-Gene 调控对。
• DORC 基因鉴定 (DORC Identification)：突破单个增强子调控的局限，寻找受大量（如 ≥ 10 个）顺式调控元件协同控制的超级靶基因。这些 DORC 基因通常是决定细胞命运、维持细胞状态或驱动疾病进展的核心调控枢纽。
• DORC 活性评估与细胞特异性分析 (Activity Scoring & Differential Analysis)：将每个 DORC 基因关联的多个 Peak 的开放度信号进行整合，量化其在单细胞水平上的综合调控活性，并通过差异分析和聚类热图揭示决定特定细胞状态的深层生物学机制。

suppressPackageStartupMessages(suppressWarnings({
    library(future)
    library(Seurat)
    library(BiocGenerics)
    library(S4Vectors)
    library(IRanges)
    library(Signac)
    library(ComplexHeatmap)
    library(Biobase)
    library(AnnotationDbi)
    library(stringr)
    library(dplyr)
    library(foreach)
    library(doParallel)
}))

2. 输入文件与参数配置

2.1 输入文件要求

本分析教程基于已经初步处理好的 Seurat 多组学对象进行 DORC (Domains of Regulatory Chromatin) 分析。开始分析前，需要准备以下几类核心输入文件：

• input.rds：预处理完成的 Seurat 对象文件。该对象必须同时包含 RNA 和 ATAC 两个 Assay，并且 ATAC 模态中需要已经完成基本的 Peak Calling，以及降维聚类分析。
• meta.tsv：细胞元数据文件（制表符分隔，可选）。用于提供额外的细胞注释信息。必须包含用于定义细胞类型和样本分组的对应列，以便后续按需进行细胞亚群的精准提取和降维聚类分析。
• genome.fa：对应物种的参考基因组 FASTA 文件。主要作用是根据 ATAC peaks 的基因组坐标信息，提取出对应的真实 DNA 序列，为后续计算 Peak 的 GC 含量、长度等序列特征（RegionStats），从而构建 Peak-Gene 相关性的背景零模型提供底层序列支持。

2.2 核心参数配置

以下参数用于控制文件路径、过滤阈值及分析细节，请根据您的实际数据进行修改：

# rds: 输入的 Seurat 对象 RDS 文件路径，必须包含预处理好的 RNA 和 ATAC 数据。
rds = "/path/to/input.rds"

# meta: 细胞元数据文件路径（TSV格式），用于提供额外的分组或样本信息。若无需添加，可留空。
meta = "/path/to/meta.tsv"

# genome_file: 对应物种的参考基因组 FASTA 文件路径，用于提取序列做 Motif 分析。
genome_file = "/path/to/genome.fa"

# dorc_peak_thresh: 界定 DORC 基因的阈值（通常设定为关联 ≥ 10 个 Peak）。
dorc_peak_thresh = "10"

# clusters_col: 指定 metadata 中哪一列作为细胞亚群的分组依据（如细胞类型）。
clusters_col = "Celltype"

# celltypes: 逗号分隔的字符串，指定要纳入分析的细胞亚群名称。
celltypes = "B cells,CMP,Dividing B cells,Erythroblast,NK cells,pDC,Plasma Cells,Pro B cells,T cells"

# sample_col: 指定 metadata 中哪一列作为样本批次的区分依据。
sample_col = "Sample"

# samples: 逗号分隔的字符串，指定要纳入分析的样本名称。
samples = "25030508_pbmc_1_arc,XYRD_pbmc_2_arc"

# filter: 逻辑值字符串（"TRUE"/"FALSE"），是否在分析前对 Peak-Gene 关联对进行显著性和方差过滤。
filter = "FALSE"

# corCutOff: 相关性得分过滤阈值（绝对值）。仅在 filter 为 TRUE 时生效。
corCutOff = ""

# pCutOff: P-value 显著性过滤阈值。仅在 filter 为 TRUE 时生效。
pCutOff = ""

# downsample: 逻辑值字符串（"True"/"False"），是否在分析前对每个细胞亚群进行随机降采样以节省内存。
downsample = "True"

# downsample_num: 降采样数量，每个细胞亚群最多保留的细胞数。仅当 downsample 为 "True" 时生效。
downsample_num = "500"

# 这里是对上面的参数做一些处理，不需要做任何改动
dorc_peak_thresh <- as.numeric(dorc_peak_thresh)

outdir <- "./result/"
dir.create(outdir, recursive = TRUE)

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

filter <- as.logical(filter)
corCutOff <- as.numeric(corCutOff)
pCutOff <- as.numeric(pCutOff)

3. 数据加载与预处理

3.1 初始化环境与并行计算

在正式加载数据前，设置并行的核心数和内存限制，以加速后续的分析步骤。

# 初始化并行计算框架
if (requireNamespace("future", quietly = TRUE)) {
    cores <- parallelly::availableCores()
    workers <- as.integer(floor(cores * 0.6))
    
    if (.Platform$OS.type == "unix") {
        future::plan(future::multicore, workers = workers)
    } else {
        future::plan(future::multisession, workers = workers)
    }
    
    message("Footprint 并行计算已启动，总核心数=", cores, "，使用 ", workers, " 个工作线程")
}

# 动态设置 future 全局变量大小限制为系统内存的 60%
if (requireNamespace("memuse", quietly = TRUE)) {
    total_mem <- memuse::Sys.meminfo()$totalram
    total_mem_bytes <- as.numeric(total_mem)  # 已经是字节了！
    future_mem <- total_mem_bytes * 0.6
    future_mem_gb <- future_mem / 1024^3
    total_mem_gb <- total_mem_bytes / 1024^3
} else {
    # 如果没有 memuse 包，使用默认 50GB
    future_mem <- 50 * 1024^3
    future_mem_gb <- 50
    total_mem_gb <- 50 / 0.6
}
options(future.globals.maxSize = future_mem)
message("系统总内存: ", round(total_mem_gb, 1), " GB，future 全局变量内存限制设置为: ", round(future_mem_gb, 1), " GB (60%)")

3.2 加载 Seurat 对象与元数据

读取预处理好的多组学对象，并按需整合额外的细胞分类信息。

obj <- readRDS(rds)
if (meta != ""){
  meta <- read.table(meta,header=T,sep='\t',check.names=F)
  #rownames(meta) <- meta$barcodes
  rownames(meta) <- meta$barcode
  obj <- AddMetaData(obj, meta)
}
Idents(obj) <- obj@meta.data[[clusters_col]]

3.3 细胞类型筛选与下采样平衡

为了降低由于某些细胞类型数量极其庞大而导致的假阳性相关性，以及控制整体分析的内存消耗，我们可以根据 celltypes 参数提取特定细胞群，并开启 downsample 机制对过大的细胞群进行随机抽样。

if (exists("downsample") && downsample != "" && !is.na(as.logical(downsample)) && as.logical(downsample)) {
  if (exists("downsample_num") && downsample_num != "" && !is.na(as.numeric(downsample_num)) && as.numeric(downsample_num) > 1) {
    obj <- subset(obj, downsample = as.numeric(downsample_num))
  }
}

obj <- subset(obj, subset = !!sym(clusters_col) %in% celltypes)
obj <- subset(obj, subset = !!sym(sample_col) %in% samples)

v <- obj@meta.data[[clusters_col]]
sorted_levels <- str_sort(unique(as.character(v)), numeric = TRUE, na_last = TRUE)
obj@meta.data[[clusters_col]] <- factor(v, levels = sorted_levels, ordered = TRUE)

3.4 参考基因组处理

为后续的 Peak-to-Gene 关联分析中的序列特征提取（RegionStats）做准备，需要先处理基因组序列文件，读取 genome.fa 并将其处理为 DNAStringSet 对象，同时对齐对象中 ATAC peaks 的染色体名称（保留共有的染色体信息）。

DefaultAssay(obj) <- "ATAC"
Idents(obj) <- clusters_col

library(Biostrings)
genome <- readDNAStringSet(genome_file)
names(genome) <- sub(" .*", "", names(genome)) 
keep_chromosomes <- unique(as.character(seqnames(granges(obj[["ATAC"]]))))
genome <- genome[which(names(genome) %in% keep_chromosomes)]
obj@assays$ATAC@ranges <- keepSeqlevels(granges(obj),keep_chromosomes, pruning.mode="coarse")

4. 细胞群间 Peak-to-Gene 关联分析

Peak-to-Gene (LinkPeaks) 关联分析是进行 DORC 鉴定的核心前置步骤。

• 为什么 DORC 分析必须先做 LinkPeaks？

DORC (调控染色质结构域) 的本质是“受极其密集的顺式调控元件协同控制的基因”。因此，要找到这些超级调控枢纽，我们必须先在全基因组范围内穷举并验证每一个 Peak 与每一个潜在靶基因之间是否真的存在调控关系。

传统的线性距离假设（例如简单认为距离 TSS 最近的 Peak 就是调控该基因的增强子）存在巨大局限性，因为染色质三维折叠使得远端增强子也能跨越数百 kb 调控靶基因。因此，本步骤采用 “共变性驱动” 的策略：利用 Signac 提供的 LinkPeaks 算法，在跨细胞（Pseudobulk）层面上，计算每个 Peak 的可及性与附近基因表达量之间的统计学相关性（Correlation）。

• 核心逻辑： 如果某个 Peak 开放度的变化始终伴随着特定靶基因表达量的同步波动，我们就认为该 Peak 是驱动该基因转录的真实增强子。只有当这层最基础的“1对1”调控连接网络构建完成后，我们才能在后续步骤中去统计“哪些基因关联了异常庞大（≥10）的调控元件”，从而完成最终的 DORC 基因鉴定。

# 检查peak-to-gene链接 
invisible(capture.output({
  suppressMessages(suppressWarnings({ 
    links_all <- Signac::Links(obj[["ATAC"]]) 
    
    if (length(links_all) == 0) { 
    
      obj <- RegionStats( 
        object = obj, 
        genome = genome, 
        assay = "ATAC" 
      ) 
      
      obj <- LinkPeaks( 
          object = obj, 
          peak.assay = "ATAC", 
          expression.assay = "RNA" 
      ) 
      
      links_all <- Signac::Links(obj[["ATAC"]]) 
      
      if (length(links_all) == 0) { 
        stop("未找到peak-to-gene链接") 
      } 
    } 
    
    write.table( 
      as.data.frame(links_all), 
      file = file.path(outdir, "Gene2PeakLinks.xls"), 
      sep = "\t", 
      row.names = FALSE, 
      col.names = TRUE, 
      quote = FALSE 
    )          
  }))
}))

saveRDS(obj, file.path(outdir, "output.rds"))

5. DORC 基因鉴定与分析

这是本教程中最核心的分析步骤。在完成全局 Peak-to-Gene 关联分析的基础上，我们将进一步引入 DORC (Domains of Regulatory Chromatin, 调控染色质结构域) 的概念，对基因组进行超级调控枢纽的扫描。以实现DORC的鉴定与下游分析和可视化。

## DORC分析
#.ts_log("开始DORC分析（使用 links_all）")
dorc_use_filter <- filter

links_dorc <- as.data.frame(links_all)
links_dorc <- links_dorc[, c("peak", "gene", "score", "pvalue")]

5.1 DORC 基因鉴定

对所有显著的 Peak-Gene 连接进行系统性盘点，统计每个靶基因在全基因组范围内所关联的顺式调控元件（Peak）总数。设定严格的背景阈值（由 dorc_peak_thresh 参数控制，通常设定为关联 ≥ 10 个 Peak）。当特定基因关联的调控元件数量显著超越该阈值时，该基因即被定义为 DORC 基因。

#添加DORC识别图
plotGeneLinkCountRank <- function(links, dorc_peak_thresh = 10, label_n = 10, seed = 1) {
  set.seed(seed)
  link_df <- data.frame(
    gene = links$gene,
    peak = links$peak,
    stringsAsFactors = FALSE
  )
  link_df <- link_df[!is.na(link_df$gene) & !is.na(link_df$peak), , drop = FALSE]
  
  # 统计并按 link_count 升序排列 (rank 由小到大)
  gene_link_df <- link_df %>%
    dplyr::distinct(gene, peak) %>%
    dplyr::count(gene, name = "link_count") %>%
    dplyr::arrange(link_count) %>%
    dplyr::mutate(
      rank = dplyr::row_number(),
      is_dorc = link_count > dorc_peak_thresh
    )
    
  n_dorc <- sum(gene_link_df$is_dorc)
  # 因为现在是升序排列，DORC 都在图的右侧
  # 分界线应该在总基因数减去 DORC 数量的位置
  total_genes <- nrow(gene_link_df)
  split_x <- total_genes - n_dorc + 0.5
  ymax <- max(gene_link_df$link_count, na.rm = TRUE)
  
  # 选择 rank 最大（link_count 最大）的 top 基因用于加标签
  dorc_genes_df <- gene_link_df %>% dplyr::filter(is_dorc == TRUE)
  label_data <- data.frame()
  if (nrow(dorc_genes_df) > 0) {
    n_sample <- min(label_n, nrow(dorc_genes_df))
    label_data <- dorc_genes_df %>%
      dplyr::arrange(dplyr::desc(link_count)) %>%
      dplyr::slice_head(n = n_sample)
  }
  
  # 定义中等蓝色
  dorc_color <- "#4A90E2"
  
label_x <- total_genes / 2

p <- ggplot2::ggplot(gene_link_df, ggplot2::aes(x = rank, y = link_count)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_point(ggplot2::aes(color = is_dorc), size = 1.2, alpha = 0.9) +
  ggplot2::geom_vline(xintercept = split_x, linetype = "dashed", color = dorc_color) +
  ggplot2::geom_hline(yintercept = dorc_peak_thresh, linetype = "dashed", color = "grey50") +
  ggplot2::scale_color_manual(values = c("TRUE" = dorc_color, "FALSE" = "grey70"), guide = "none") +
  ggplot2::annotate(
    "text",
    x = label_x,
    y = ymax,
    label = paste0("DORC n=", n_dorc),
    vjust = -0.6,
    hjust = 0.5,
    color = dorc_color,
    size = 4
  ) +
  ggplot2::labs(
    x = "Rank of genes by linked peak count (Ascending)",
    y = "Number of linked peaks per gene"
  ) +
  ggplot2::theme_classic() +
  ggplot2::theme(
    legend.position = "none",
  # text = element_text(family = "sans"),
    axis.title = ggplot2::element_text(size = 14),
    axis.text = ggplot2::element_text(size = 12),
    plot.title = ggplot2::element_text(hjust = 0.5, size = 14),
    axis.text.x = ggplot2::element_blank(),  # 隐藏x轴刻度标签
    axis.ticks.x = ggplot2::element_blank()  # 隐藏x轴刻度线
  )
    
  # 添加 top rank 的 DORC 基因标签
  if (nrow(label_data) > 0) {
    p <- p + ggrepel::geom_text_repel(
      data = label_data,
      ggplot2::aes(label = gene),
      color = "black",
      size = 3,
      box.padding = 0.5,
      point.padding = 0.3,
      force = 0.03,
      max.overlaps = 50,
      segment.color = "grey50",
      nudge_x = -0.1 * total_genes, # 尽量把标签往左边推一点，避免出界
      direction = "y"
    )
  }
  
  return(p)
}

options(repr.plot.width = 10, repr.plot.height = 5)
print(plotGeneLinkCountRank(links_all, dorc_peak_thresh = dorc_peak_thresh))

💡 图表解读DORC基因识别排序图 (Peak2GeneGeneLinkCountRank)：

• 横轴（X轴）：代表所有靶基因。基因按照其关联的 Peak 数量由少到多（升序）排列。
• 纵轴（Y轴）：代表每个靶基因在全基因组范围内显著关联的 Peak 绝对数量。
• 虚线含义：水平灰色虚线代表界定 DORC 基因的截断阈值（默认关联 ≥ 10 个 Peak）；垂直蓝色虚线将右侧满足阈值条件的超级调控基因（蓝色点）与左侧的普通基因（灰色点）区分开来。
• 文本标注：图表顶部标注了识别出的 DORC 基因总数，同时在右侧密集区随机抽取了部分具有代表性的 DORC 基因名称进行加粗标记，帮助研究者快速锁定潜在的关键核心驱动基因。

5.2 DORC 活性评估

计算每个细胞中靶向同一 DORC 基因的所有 Peak 的累积可及性，生成单细胞水平的 DORC 活性评分（Accessibility Score）矩阵。 这不仅仅是一个简单的计数过程。我们会将属于同一个 DORC 基因的所有相关 Peak 的 ATAC 可及性信号（Pseudobulk Counts）进行加权或合并计算，从而为每个细胞群生成一个综合的、代表该 DORC 基因整体上游调控强度的“DORC Score”。

invisible(capture.output({
  suppressMessages(suppressWarnings({
    
    links_dorc <- links_dorc[!is.na(links_dorc$peak) & !is.na(links_dorc$gene), , drop = FALSE] 
    
    if (dorc_use_filter) { 
      links_dorc <- links_dorc[abs(links_dorc$score) >= corCutOff & links_dorc$pvalue <= pCutOff, , drop = FALSE] 
    } 
    
    if (nrow(links_dorc) == 0) { 
      stop("DORC分析失败：links_all 在当前阈值下为空") 
    } 
    
    atac_counts <- GetAssayData(obj, assay = "ATAC", slot = "counts") 
    if (!inherits(atac_counts, "dgCMatrix")) atac_counts <- as(atac_counts, "dgCMatrix") 
    
    peak_idx <- match(links_dorc$peak, rownames(atac_counts)) 
    if (any(is.na(peak_idx))) { 
      stop("DORC分析失败：links_all 中存在不在 ATAC counts 中的 peak") 
    } 
    
    peak_idx_by_gene <- split(peak_idx, links_dorc$gene) 
    peak_idx_by_gene <- lapply(peak_idx_by_gene, unique) 
    peak_count <- vapply(peak_idx_by_gene, length, integer(1)) 
    peak_idx_by_gene <- peak_idx_by_gene[peak_count > dorc_peak_thresh] 
    
    if (length(peak_idx_by_gene) == 0) { 
      stop(paste0("DORC分析失败：没有基因满足链接peak数 >", dorc_peak_thresh)) 
    } 
    
    peak_count <- vapply(peak_idx_by_gene, length, integer(1)) 
    gene_names <- names(peak_idx_by_gene) 
    
    cell_sums <- Matrix::colSums(atac_counts) 
    cell_sums[cell_sums == 0] <- 1 
    scale_diag <- Matrix::Diagonal(x = 1e6 / cell_sums) 
    atac_cpm <- atac_counts %*% scale_diag 
    
    dorc_mat <- matrix(0, nrow = ncol(atac_cpm), ncol = length(gene_names)) 
    rownames(dorc_mat) <- colnames(atac_counts) 
    colnames(dorc_mat) <- gene_names 
    
    for (i in seq_along(gene_names)) { 
      idx <- peak_idx_by_gene[[i]] 
      if (length(idx) == 1) { 
        dorc_mat[, i] <- as.numeric(atac_cpm[idx, ]) 
      } else { 
        dorc_mat[, i] <- as.numeric(Matrix::colSums(atac_cpm[idx, , drop = FALSE])) 
      } 
    } 
    
    peak_sum_df <- data.frame( 
      peak_count = peak_count, 
      row.names = gene_names 
    ) 
    peak_sum_sorted <- peak_sum_df[order(peak_sum_df$peak_count, decreasing = TRUE), , drop = FALSE] 
    
    output_dorc_peaks <- file.path(outdir, paste0("DORC_peaks_thresh", dorc_peak_thresh, "_sum_sorted.tsv")) 
    write.table(peak_sum_sorted, output_dorc_peaks, sep = "\t", row.names = TRUE, col.names = FALSE, quote = FALSE) 
    
    output_dorc <- file.path(outdir, paste0("DORC_scores_thresh", dorc_peak_thresh, ".tsv")) 
    write.table(dorc_mat, output_dorc, sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE) 
    
    saveRDS(dorc_mat, file.path(outdir, paste0("DORC_scores_thresh", dorc_peak_thresh, ".rds"))) 
    
  }))
}))

5.3 DORC 活性差异分析

基于前一步构建的单细胞 DORC 活性评分矩阵，利用 Seurat 的 FindAllMarkers 函数执行差异活性分析。该步骤旨在寻找在各个细胞亚群中特异性表现出高调控活性的 DORC 基因，从而鉴定出维持特定细胞群稳态或驱动其分化的标志性调控枢纽。

invisible(capture.output({
  suppressMessages(suppressWarnings({
    
    dorc_df <- as.data.frame(dorc_mat, check.names = FALSE) 
    if (!all(colnames(obj) %in% rownames(dorc_df))) { 
      stop("DORC下游分析失败：dorc_mat 与 Seurat 对象细胞不一致") 
    } 
    dorc_df <- dorc_df[colnames(obj), , drop = FALSE] 
    
    dorc_meta <- dorc_df 
    colnames(dorc_meta) <- paste0(colnames(dorc_meta), "_dorc") 
    existing_dorc_meta <- intersect(colnames(obj@meta.data), colnames(dorc_meta)) 
    if (length(existing_dorc_meta) > 0) { 
      obj@meta.data <- obj@meta.data[, setdiff(colnames(obj@meta.data), existing_dorc_meta), drop = FALSE] 
    } 
    obj@meta.data <- cbind(obj@meta.data, dorc_meta[rownames(obj@meta.data), , drop = FALSE]) 
    
    dorc_assay_mat <- t(as.matrix(dorc_df)) 
    if (!inherits(dorc_assay_mat, "dgCMatrix")) { 
      dorc_assay_mat <- as(dorc_assay_mat, "dgCMatrix") 
    } 
    obj[["DORC"]] <- CreateAssayObject(counts = dorc_assay_mat) 
    
    if (!(clusters_col %in% colnames(obj@meta.data))) { 
      stop("DORC下游分析失败：meta.data中不存在细胞类型列 ", clusters_col) 
    } 
    Idents(obj) <- as.factor(as.character(obj@meta.data[[clusters_col]])) 
    DefaultAssay(obj) <- "DORC" 
    
    # 这里是产生红色刷屏的地方，现在被外层的 suppressMessages 吃掉了
    markers_dorc <- FindAllMarkers( 
      object = obj, 
      assay = "DORC", 
      only.pos = TRUE, 
      logfc.threshold = 0, 
      min.pct = 0.05 
    ) 
    
    if (!"gene" %in% colnames(markers_dorc)) { 
      markers_dorc$gene <- rownames(markers_dorc) 
    } 
    output_dorc_markers <- file.path(outdir, "DORC_FindAllMarkers.tsv") 
    write.table(markers_dorc, output_dorc_markers, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)     
    markers_sig <- markers_dorc[!is.na(markers_dorc$p_val_adj) & markers_dorc$p_val_adj < 0.05 & markers_dorc$avg_log2FC > 0, , drop = FALSE] 
    if (nrow(markers_sig) == 0) { 
      stop("DORC下游分析失败：没有显著差异DORC（p_val_adj < 0.05）") 
    } 
    
    # 提取每个 cluster 中 avg_log2FC 最高的 top 10 基因 
    markers_top10 <- markers_sig %>% 
      group_by(cluster) %>% 
      slice_max(n = 10, order_by = avg_log2FC, with_ties = FALSE) %>% 
      ungroup() 
    
    output_top10 <- file.path(outdir, "DORC_top10_per_cluster_markers.tsv") 
    write.table(markers_top10, output_top10, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE) 

  }))
}))

head(markers_sig)

A data.frame: 6 × 7

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<fct>	<chr>
DTX1	2.397218e-189	2.463573	0.924	0.462	4.569098e-186	B cells	DTX1
RASAL1	1.771071e-171	2.109955	0.942	0.588	3.375661e-168	B cells	RASAL1
CFAP73	2.060654e-165	2.016315	0.944	0.598	3.927607e-162	B cells	CFAP73
IGLC1	2.272933e-163	1.590690	0.994	0.796	4.332209e-160	B cells	IGLC1
TCL6	8.499601e-160	1.902479	0.944	0.530	1.620024e-156	B cells	TCL6
SLC8B1	5.623266e-158	1.932736	0.942	0.637	1.071795e-154	B cells	SLC8B1

💡 解读指南：DORC (Domains of Regulatory Chromatin) 差异分析结果表

该表展示了通过 Seurat FindAllMarkers 鉴定出的在各个细胞亚群中具有特异性高调控活性的 DORC 基因。DORC 得分反映了基因附近多个开放染色质区域信号的累加综合效应，是强增强子调控基因的重要标志。建议联合 p_val_adj 和 avg_log2FC 进行优先级排序。

• gene：DORC 基因的名称，即那些受多个顺式调控元件强力驱动的靶基因。
• cluster：该 DORC 基因表现出显著高调控活性的目标细胞亚群名称或编号。
• p_val：差异分析的原始显著性 P 值。
• avg_log2FC：平均 Log2 差异倍数。正值表示该 DORC 基因在当前 cluster 中的综合调控活性高于其他细胞群。数值越大，代表该组的特异性调控强度越高。
• pct.1：该基因作为 DORC 在当前 cluster 中的检出（激活）比例，范围在 0~1 之间。
• pct.2：该基因作为 DORC 在其他所有细胞亚群背景中的检出（激活）比例，范围在 0~1 之间。
• p_val_adj：经多重检验校正后的 P 值，是判断该基因是否为该亚群显著 DORC Marker 的金标准。

5.4 DORC 活性 UMAP 空间可视化

基于前面提取的各个细胞亚群的特异性 Top 10 DORC 基因，通过 FeaturePlot 在 UMAP 降维空间中直观呈现其单细胞水平的活性分布。由于 DORC 评分差异可能极大，这里预先进行了 log10 转换和分位数截断以平滑极值。

DefaultAssay(obj) <- "DORC"

# 你要展示的DORC列表（例如 top10）
features_to_show <- markers_top10$gene
features_to_show <- unique(features_to_show)
features_to_show <- features_to_show[features_to_show %in% rownames(obj[["DORC"]])]

plot_output_dir <- file.path(outdir, "DORC_FeaturePlot_5in1")
dir.create(plot_output_dir, recursive = TRUE, showWarnings = FALSE)

plot_dorc_feature <- function(gene, obj) {
  # 提取原始数值
  values <- as.numeric(GetAssayData(obj, assay = "DORC", slot = "data")[gene, ])
  values <- values[is.finite(values)]
  if (length(values) == 0) return(NULL)

  # 关键修复：对极大的 DORC 得分进行 log10(x+1) 缩放平滑
  scaled_values <- log10(values + 1)
  
  # 将平滑后的数值临时放回 Seurat 对象的 meta.data 中用于画图
  plot_col_name <- paste0("Log10_", gsub("-", "_", gene))
  obj@meta.data[[plot_col_name]] <- scaled_values

  p1 <- as.numeric(quantile(scaled_values, 0.01, na.rm = TRUE))
  p99 <- as.numeric(quantile(scaled_values, 0.99, na.rm = TRUE))

  # 定义经典的蓝-黄-红渐变色带
  custom_colors <- c("#313695", "#4575b4", "#74add1", "#abd9e9", "#e0f3f8", 
                     "#ffffbf", "#fee090", "#fdae61", "#f46d43", "#d73027", "#a50026")

  FeaturePlot(
    object = obj,
    features = plot_col_name,  # 使用新生成的 Log10 列画图
    reduction = "umap",
    min.cutoff = p1,
    max.cutoff = p99,
    pt.size = 0.3
  ) +
    ggplot2::scale_color_gradientn(
      colors = custom_colors,
      name = "log10(DORC+1)" # 修改图例名称让其更专业
    ) +
    ggplot2::ggtitle(gene) +
    ggplot2::theme(plot.title = ggplot2::element_text(size = 10, hjust = 0.5))
}

save_and_encode_plot <- function(plot, filename_base, width = 26, height = 3.8) {
  full_path_pdf <- paste0(filename_base, ".pdf")
  full_path_png <- paste0(filename_base, ".png") # 直接生成同名PNG文件路径

  ggplot2::ggsave(full_path_pdf, plot, width = width, height = height)
  ggplot2::ggsave(full_path_png, plot, width = width, height = height, dpi = 100)

  # 检查PNG是否生成成功，直接返回文件路径而不是转码Base64
  if (file.exists(full_path_png)) {
    return(list(png_path = full_path_png, pdf_path = full_path_pdf))
  }
  NULL
}

suppressMessages(suppressWarnings({
# 每5个基因分组
feature_groups <- split(
  features_to_show,
  ceiling(seq_along(features_to_show) / 5)
)

plot_list <- list()

for (i in seq_along(feature_groups)) {
  genes <- feature_groups[[i]]

  p_list <- lapply(genes, function(g) plot_dorc_feature(g, obj))
  p_list <- p_list[!vapply(p_list, is.null, logical(1))]
  if (length(p_list) == 0) next

  # 补齐到5格，保持固定1行5列
  if (length(p_list) < 5) {
    p_list <- c(p_list, replicate(5 - length(p_list), patchwork::plot_spacer(), simplify = FALSE))
  }

  combined_plot <- patchwork::wrap_plots(p_list, ncol = 5) +
    patchwork::plot_annotation(title = paste0("DORC FeaturePlot Group ", i))

  file_base <- file.path(plot_output_dir, paste0("DORC_group_", sprintf("%03d", i)))
  res <- save_and_encode_plot(combined_plot, file_base, width = 26, height = 3.8)
}

}))

options(repr.plot.width = 23, repr.plot.height = 4)
combined_plot

💡 解读指南：差异 DORC 基因 UMAP 特征映射图 (FeaturePlot)

该图示例展示各个核心驱动基因（DORC 基因）在不同单细胞群体中的空间分布与相对调控活性。其他全部活性差异的DORC放在.result/DORC_FeaturePlot_5in1目录下。

• 散点（细胞）：图中的每一个散点代表一个单细胞，其在二维空间中的相对位置反映了细胞群体之间的整体转录组或染色质状态相似性。
• 颜色梯度（活性得分）：颜色反映了特定 DORC 基因在该细胞中累积可及性评分的对数平滑值（log10(DORC+1)）。
  • 深蓝色：代表未激活或极低活性。
  • 黄绿色：代表中等过渡态。
  • 橙色至深红色：代表该基因关联的多个增强子元件在该细胞中高度开放，调控活性极强。
• 数值平滑与截断：由于原始的 DORC 峰计数往往存在极端高值，分析中预先对其进行了 Log10 对数转换以平滑极值差异；同时，颜色映射的上下限采用了平滑后数据的 1% 和 99% 分位数进行截断，以确保微小的生物学变异模式得以清晰展现。
• 组别展示：图片下方的数据表以每组 5 个基因的格式进行分页展示。点击缩略图即可在新标签页中查看高分辨率的 PDF 矢量图。

5.5 DORC 细胞群水平活性热图 (Pseudobulk Heatmap)

为了从全局视角展示各个细胞亚群的核心 DORC 基因活性模式，我们将单细胞级别的 DORC 矩阵按细胞亚群进行汇总拟合（Pseudobulk），并计算平均活性。然后提取每个亚群的 Top 10 特异性 DORC 基因，进行 Z-score 标准化后绘制层次聚类热图。

# ---- Code cell 61 ----
markers_top10_by_cluster <- markers_sig %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE) %>%
  slice_head(n = 10) %>%
  ungroup()

output_top10_by_cluster <- file.path(outdir, "DORC_top10_by_cluster.tsv")
write.table(markers_top10_by_cluster, output_top10_by_cluster, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)

selected_dorc_genes <- unique(markers_top10_by_cluster$gene)
selected_dorc_genes <- selected_dorc_genes[selected_dorc_genes %in% rownames(obj[["DORC"]])]
if (length(selected_dorc_genes) == 0) {
  stop("DORC下游分析失败：用于拟bulk绘图的DORC基因为空")
}

dorc_selected_cells_by_gene <- t(as.matrix(GetAssayData(obj, assay = "DORC", slot = "data")[selected_dorc_genes, , drop = FALSE]))
cell_groups <- as.character(obj@meta.data[[clusters_col]])
names(cell_groups) <- rownames(obj@meta.data)
cell_groups <- cell_groups[rownames(dorc_selected_cells_by_gene)]
valid_cells <- !is.na(cell_groups) & cell_groups != ""
dorc_selected_cells_by_gene <- dorc_selected_cells_by_gene[valid_cells, , drop = FALSE]
cell_groups <- cell_groups[valid_cells]

dorc_pb_sum <- rowsum(dorc_selected_cells_by_gene, group = cell_groups, reorder = FALSE)
dorc_pb_n <- as.numeric(table(cell_groups)[rownames(dorc_pb_sum)])
dorc_pb_mean <- dorc_pb_sum / dorc_pb_n
dorc_pb_gene_by_celltype <- t(dorc_pb_mean)

output_pb <- file.path(outdir, "DORC_pseudobulk_top10ByCluster.tsv")
write.table(dorc_pb_gene_by_celltype, output_pb, sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)

dorc_pb_scaled <- t(scale(t(dorc_pb_gene_by_celltype)))
dorc_pb_scaled[is.na(dorc_pb_scaled)] <- 0



# 仅用于行注释：根据细胞类型数量动态决定每个 cluster 标注的特异性基因个数
n_celltypes <- ncol(dorc_pb_scaled)
labels_per_cluster <- if (n_celltypes <= 5) {
  5L
} else if (n_celltypes < 8) {
  4L
} else if (n_celltypes <= 10) {
  3L
} else if (n_celltypes <= 15) {
  2L
} else {
  1L
}

markers_ranked_by_cluster <- markers_top10_by_cluster %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE) %>%
  ungroup() %>%
  filter(gene %in% rownames(dorc_pb_scaled))

clusters <- unique(as.character(markers_ranked_by_cluster$cluster))
used_genes <- character()
anno_labels <- character()
for (cl in clusters) {
  genes <- as.character(markers_ranked_by_cluster$gene[markers_ranked_by_cluster$cluster == cl])
  genes <- genes[!is.na(genes) & genes != ""]
  genes <- genes[genes %in% rownames(dorc_pb_scaled)]
  preferred <- genes[!genes %in% used_genes]
  picked <- head(preferred, labels_per_cluster)
  if (length(picked) < labels_per_cluster) {
    picked <- c(picked, head(setdiff(genes, picked), labels_per_cluster - length(picked)))
  }
  used_genes <- c(used_genes, picked)
  anno_labels <- c(anno_labels, picked)
}

row_hc <- stats::hclust(stats::dist(dorc_pb_scaled))
anno_df <- data.frame(
  label = anno_labels,
  row_index = match(anno_labels, rownames(dorc_pb_scaled)),
  stringsAsFactors = FALSE
)
# anno_mark() 的 at 需要原始矩阵的行索引，不能传聚类后的显示位置，
# 否则在 cluster_rows 生效后会发生二次错位。
anno_df <- anno_df[!is.na(anno_df$row_index), , drop = FALSE]
anno_df <- anno_df[!duplicated(anno_df$row_index), , drop = FALSE]
anno_at <- anno_df$row_index
anno_labels <- anno_df$label

row_anno <- NULL
anno_width_cm <- 0
if (length(anno_at) > 0) {
  anno_width_cm <- max(4.5, min(8, max(nchar(anno_labels)) * 0.25))
  row_anno <- ComplexHeatmap::rowAnnotation(
    mark = ComplexHeatmap::anno_mark(
      at = anno_at,
      labels = anno_labels,
      which = "row",
      side = "right",
      labels_gp = grid::gpar(fontsize = 12),
      lines_gp = grid::gpar(col = "grey50", lwd = 0.8)
    ),
    width = grid::unit(anno_width_cm, "cm")
  )
}

output_pb_pdf <- file.path(outdir, "DORC_pseudobulk_top10ByCluster_heatmap.pdf")
pdf(output_pb_pdf, width = 12, height = 9)
ht <- ComplexHeatmap::Heatmap(
  dorc_pb_scaled,
  name = "DORC Z-score",
  col = circlize::colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  cluster_rows = row_hc,
  cluster_columns = TRUE,
  show_row_names = FALSE,
  show_row_dend = FALSE,
  show_column_dend = FALSE,  # 不显示列聚类树
  show_column_names = TRUE,
  right_annotation = row_anno,
  column_names_gp = grid::gpar(fontsize = 12)
)
draw(ht)
invisible(dev.off())

saveRDS(obj, file.path(outdir, "output_with_dorc_downstream.rds"))

options(repr.plot.width = 15, repr.plot.height = 8)
draw(ht)

💡 解读指南：DORC 基因活性聚类热图

该热图展示了各细胞亚群特异性核心 DORC 基因的全局活性图谱。数据已在细胞群体层面进行了拟合（Pseudobulk 平均化），并对基因活性进行了 Z-score 行标准化处理。

• 纵轴（行）：代表从每个细胞亚群中提取的基于 avg_log2FC 排序的 Top 10 最显著的差异 DORC 基因。左侧的树状图反映了这些基因在全局调控模式上的相似性聚类。
• 横轴（列）：代表数据中定义的不同细胞亚群。列上的树状图反映了各细胞亚群在核心调控网络层面的亲缘关系。
• 颜色梯度（Z-score）：反映了特定基因在不同细胞群中的相对调控活性。
  • 红色：表示该 DORC 基因在该细胞群中的活性显著高于平均水平。
  • 蓝色：表示活性显著低于平均水平。
  • 白色：表示接近平均值。
• 核心发现：通过观察热图上的“红色高亮色块矩阵”，可以直观地锁定驱动特定细胞类群发育或维持其生物学特性的核心转录调控枢纽。