ATAC + RNA 多组学：Peak2Gene 分析
时长: 48 分钟
字数: 11.4k 字
更新: 2026-07-07
阅读: 0 次
1. 教程简介
本教程基于 Seurat 和 Signac 框架，用于对 SeekArc 单细胞多组学数据（RNA + ATAC）进行深度的 Peak-to-Gene 关联分析。通过打破线性距离的限制，基于 共变性假设，重构高分辨率的顺式调控网络。 该教程面向已完成基础质控、降维聚类及 Peak Calling 的单细胞多组学数据。核心分析目标包括：
- 计算 Peak-Gene 关联：计算染色质开放区域（Peak）与基因表达（RNA）之间的共变性，识别潜在的顺式调控关系。
- 假混合（Pseudobulk）构建：为提高单细胞数据计算相关性时的信噪比和鲁棒性，按照细胞亚群构建假混合矩阵。
- 共调控模块 K-means 聚类：将具有相似动态变化模式的 Peak-Gene 关联对划分为独立的共调控簇（Cluster）。
- 联合热图可视化：并排展示 ATAC 可及性 Z-score 与 RNA 表达 Z-score，直观呈现不同细胞群中的共调控模式。
- 转录因子 Motif 富集分析：对每个共调控模块的 Peak 序列进行 TF Motif 富集，追溯上游调控因子。
- 靶基因 GO/KEGG 功能富集：对每个模块的靶基因进行通路富集，揭示调控网络参与的核心生物学功能。
注意
💡 文档中展示的图表仅为部分代表性的可视化结果。如需查看所有完整分析结果（如GO/KEGG/Motif富集等结果），请前往
./result/目录进行查阅。
R
suppressPackageStartupMessages(suppressWarnings({
library(future)
library(SeuratObject)
library(Seurat)
library(BiocGenerics)
library(S4Vectors)
library(IRanges)
library(Signac)
library(ComplexHeatmap)
library(Biobase)
library(AnnotationDbi)
library(org.Hs.eg.db)
#library(org.Mm.eg.db)
library(JASPAR2020)
library(TFBSTools)
library(stringr)
library(dplyr)
library(DOSE)
library(clusterProfiler)
library(enrichplot)
library(foreach)
library(doParallel)
library(base64enc)
library(KEGG.db)
}))R
# 初始化并行计算框架
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%)")2. 输入文件准备
2.1 输入文件要求
本分析教程基于已经初步处理好的 Seurat 多组学对象进行 Peak2Gene 关联分析。开始分析前，需要准备以下几类核心输入文件：
input.rds：预处理完成的 Seurat 对象文件。该对象必须包含RNA和ATAC两个 Assay，并且 ATAC 模态中需要已经完成基本的 Peak 调用和 Motif 注释。meta.tsv：细胞元数据文件（制表符分隔）。必须包含用于定义细胞类型和样本分组的对应列，以便后续按需进行细胞亚群的提取和聚类分析。genome.fa：对应物种的参考基因组 FASTA 文件。主要作用是根据 ATAC peaks 的基因组坐标信息，提取出对应的真实 DNA 序列。因为后续需要扫描这些序列中是否存在特定转录因子（TF）的结合位点（Motif），所以提供与单细胞比对时一致的精确参考基因组序列是进行 Motif 分析的必要前提。
2.2 核心参数配置
以下参数用于控制文件路径、过滤阈值及分析细节，请根据您的实际数据进行修改：
R
# 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"
# species: 物种名称（"human" 或 "mouse"），用于匹配正确的基因注释库和 TF 数据库。
species = "human"
# 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 = "XYRD_pbmc_2_arc,25030508_pbmc_1_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 = "1000"
# varCutOffATAC: ATAC 数据的方差分位数阈值（0~1）。过滤掉变化极小的 Peak。仅在 filter 为 TRUE 时生效。
varCutOffATAC = ""
# varCutOffRNA: RNA 数据的方差分位数阈值（0~1）。过滤掉表达量极稳定的 Gene。仅在 filter 为 TRUE 时生效。
varCutOffRNA = ""
# k: K-means 聚类的模块数量，即预期将共变特征划分为几个簇。
k = "10"
# nPlot: 热图绘制时随机抽样的连线数量上限，避免绘制过多连线导致出图崩溃。
nPlot = "1000"R
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)
varCutOffATAC <- as.numeric(varCutOffATAC)
varCutOffRNA <- as.numeric(varCutOffRNA)
k <- as.numeric(k)
nPlot <- as.numeric(nPlot)3. 数据加载与预处理
在环境、数据和参数等都设置好后，开始加载和预处理待分析的数据。
3.1 加载 Seurat 对象与元数据
读取预处理好的多组学对象，并按需整合额外的细胞分类信息。
R
obj <- readRDS(rds)
if (meta != ""){
meta <- read.table(meta,header=T,sep='\t',check.names=F)
rownames(meta) <- meta$barcode
obj <- AddMetaData(obj, meta)
}
Idents(obj) <- obj@meta.data[[clusters_col]]3.2 细胞类型筛选与下采样平衡
为了降低分析的内存消耗并平衡各细胞亚群的数据量，如果开启了 downsample 机制，会首先对对象进行降采样（限制每个亚群的最大细胞数）。随后，根据 celltypes 和 samples 参数提取指定的细胞群和样本，并将细胞亚群的因子水平按照自然顺序进行重排。
R
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 <- stringr::str_sort(unique(as.character(v)), numeric = TRUE, na_last = TRUE)
obj@meta.data[[clusters_col]] <- factor(v, levels = sorted_levels, ordered = TRUE)3.3 参考基因组处理
为后续的 TF（转录因子）分析做准备，需要先处理基因组序列文件，读取 genome.fa 并将其处理为 DNAStringSet 对象，同时对齐对象中 ATAC peaks 的染色体名称（保留共有的染色体信息）。
R
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")3.4 Motif 数据获取
在执行转录因子 Motif 富集分析之前，我们需要先获取参考的 Motif 位置频率矩阵（PFM，Position Frequency Matrix）。然后使用 AddMotifs() 将 Motif 信息添加至 ATAC Assay
为什么要提前做 AddMotifs 分析？
- 构建 Motif x Peak 匹配矩阵：
AddMotifs会扫描我们刚刚获取的 ATAC Peak DNA 序列，寻找是否存在转录因子（TF）的结合基序（由pfm提供）。这一步会在 Seurat 对象中全局存储每个 Peak 与 Motif 的匹配关系。 - 下游模块富集的基础：后续我们将 Peak 和 Gene 聚类为多个共调控模块（Cluster）。当需要对某个模块的 Peak 集合进行 Motif 富集时（推测是哪些 TF 调控了这些靶基因），
run_Motif_enrichment会直接利用这里构建好的全局匹配信息，而无需对每个聚类簇重新扫描序列，极大地提高了计算效率。
R
species_tmp <- species
if (species_tmp == "human"){
pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = 9606, all_versions = TRUE))
species = "hsa"
db = 'org.Hs.eg.db'
spe = "human"
} else if (species_tmp == "mouse"){
pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = 10090, all_versions = TRUE))
species = "mmu"
db = 'org.Mm.eg.db'
spe = "mouse"
} else {
stop(paste0("目前选项只有 'human' 或 'mouse',","不支持该物种: '", species,
"'\n请自行安装对应物种的基因注释数据库。"))
}
obj <- AddMotifs(object = obj, genome = genome, pfm = pfm)3.5 核心辅助函数定义
在开始主体分析流程之前，我们预先定义一系列核心的辅助函数。这些函数是构建高级分析流水线的基础，它们包括：
- 日志记录与进度追踪 (
.ts_log) - 元数据的统计与聚合 (
.aggregate_meta) - 单细胞数据向 Pseudobulk（假混合）对象的转化封装 (
build_pseudobulk_seurat) - 针对共变特征的 K-means 聚类算法 (
clusterRowsKmeans) - 下游模块基因的功能通路富集分析（包括
go_enrich、kegg_enrich、run_Motif_enrichment） - 综合绘图与分析引擎 (
peak2GeneHeatmap) - 基础网络评估与可视化统计图函数 (
plotPeaksPerGene,plotPeak2GeneVolcano,plotRankCorrelation)
1. .ts_log: 带时间戳的日志输出函数。目的是在长时间运行的单细胞流程中，为控制台输出添加标准的时间前缀，方便学者监控分析进度和定位耗时步骤。
R
# 带时间戳日志
.ts_log <- function(...) {
ts <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
cat(sprintf("[%s] %s\n", ts, paste0(..., collapse = "")))
}2. .majority_value: 统计众数函数。主要服务于后续的细胞聚合（Pseudobulk）。当合并多个细胞为一个伪细胞时，数值型数据可以求平均，而对于字符/因子型（如细胞类型名称），我们需要求出占比最多的类别作为该伪细胞的属性。
R
# 统计众数（用于字符/因子/逻辑型）
.majority_value <- function(x) {
x <- x[!is.na(x)]
if (length(x) == 0) return(NA)
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}3. .aggregate_meta: 汇总元数据函数。将单个细胞级别的 metadata 按照指定分组（如 Pseudobulk 块）聚合，数值型取平均，类别型取众数，确保降维后的伪样本对象仍保留准确的注释信息。
R
.aggregate_meta <- function(meta, groups) {
stopifnot(length(groups) == nrow(meta))
groups <- as.factor(groups)
lv <- levels(groups)
res_list <- lapply(lv, function(l) {
idx <- which(groups == l)
sub <- meta[idx, , drop = FALSE]
out <- lapply(sub, function(col) {
if (is.numeric(col) || is.integer(col)) {
mean(col, na.rm = TRUE)
} else if (is.logical(col)) {
as.logical(.majority_value(col))
} else if (is.factor(col)) {
as.character(.majority_value(as.character(col)))
} else if (inherits(col, "POSIXt") || inherits(col, "Date")) {
suppressWarnings(mean(as.numeric(col), na.rm = TRUE))
} else {
.majority_value(as.character(col))
}
})
out
})
df <- as.data.frame(do.call(rbind, res_list), stringsAsFactors = FALSE)
rownames(df) <- lv
df
}4. build_pseudobulk_seurat: 构建 Pseudobulk（假混合）对象的主函数。
单细胞数据本身具有极高的稀疏性（大量 0 值），如果直接计算 Peak 与 Gene 之间的相关性，噪音极大。该函数的核心目的是将属于同一细胞亚群（由 groupBy 参数指定）的若干个细胞 Count 值进行合并累加，生成一个个“超细胞（Pseudobulk）”。这样不仅能有效克服 Drop-out 效应，提高表达和可及性矩阵的信噪比，还能成百倍地减少矩阵维度，极大加快后续分析的计算速度。
R
# 主函数：构建 pseudobulk Seurat 对象
build_pseudobulk_seurat <- function(
obj,
groupBy,
pseudobulkMem = 10,
min_pseudobulk = 100,
respect_group_boundaries = TRUE,
keep_remainders = TRUE,
seed = 1,
verbose = TRUE
) {
set.seed(seed)
if (verbose) .ts_log("开始构建 pseudobulk（order_bin）")
# 基本检查
stopifnot("RNA" %in% names(obj@assays))
stopifnot("ATAC" %in% names(obj@assays))
meta <- obj@meta.data
all_cells <- colnames(obj)
# 取 counts（稀疏）
rna_counts <- GetAssayData(obj, assay = "RNA", slot = "counts")
atac_counts <- GetAssayData(obj, assay = "ATAC", slot = "counts")
if (!inherits(rna_counts, "dgCMatrix")) rna_counts <- as(rna_counts, "dgCMatrix")
if (!inherits(atac_counts, "dgCMatrix")) atac_counts <- as(atac_counts, "dgCMatrix")
# 1) 确定分组（pseudobulk ID）
if (!(groupBy %in% colnames(meta))) {
stop("meta.data 中未找到用于排序的列：", groupBy)
}
# 仅使用在 RNA 与 ATAC counts 中同时存在的细胞，构建用于排序的元数据子集
cells_in_counts <- intersect(colnames(rna_counts), colnames(atac_counts))
meta_sub <- meta[cells_in_counts, , drop = FALSE]
if (verbose) .ts_log("使用 ", groupBy, " 进行排序")
ordered_cellID <- rownames(meta_sub[order(meta_sub[[groupBy]], na.last = NA), , drop = FALSE])
# 判断 groupBy 是否离散
x <- meta_sub[[groupBy]]
uniq_n <- length(unique(x[!is.na(x)]))
discrete_group <- is.factor(x) || is.character(x) || (is.numeric(x) && uniq_n <= max(50, ceiling(0.01 * length(x))))
make_bins <- function(mem) {
bins <- list()
if (respect_group_boundaries && discrete_group) {
if (verbose) .ts_log("启用分组内分箱，确保不跨 ", groupBy, " 边界")
# 保持组的出现顺序（按整体排序后的顺序）
grp_order <- unique(as.character(meta_sub[ordered_cellID, groupBy]))
for (g in grp_order) {
if (is.na(g)) next
gi <- ordered_cellID[as.character(meta_sub[ordered_cellID, groupBy]) == g]
if (length(gi) < mem && !keep_remainders) next
n_bin_g <- if (keep_remainders) ceiling(length(gi) / mem) else floor(length(gi) / mem)
if (n_bin_g > 0) {
use_n <- (n_bin_g - 1) * mem
if (n_bin_g > 1) {
m <- matrix(gi[seq_len(use_n)], nrow = n_bin_g - 1, ncol = mem, byrow = TRUE)
bins <- c(bins, split(m, row(m)))
}
# 剩余部分（可能 < mem）
if (keep_remainders) {
rem <- gi[(use_n + 1):length(gi)]
if (length(rem) > 0) bins <- c(bins, list(rem))
} else if (length(gi) >= mem) {
last_full <- gi[(use_n + 1):(use_n + mem)]
bins <- c(bins, list(last_full))
}
}
}
} else {
# 原顺序固定窗口分箱
total_cells <- length(ordered_cellID)
if (keep_remainders) {
n_full <- floor(total_cells / mem)
if (n_full > 0) {
use_n <- n_full * mem
m <- matrix(ordered_cellID[seq_len(use_n)], nrow = n_full, ncol = mem, byrow = TRUE)
bins <- c(bins, split(m, row(m)))
}
if (total_cells > n_full * mem) {
bins <- c(bins, list(ordered_cellID[(n_full * mem + 1):total_cells]))
}
} else {
n_bin <- floor(total_cells / mem)
if (n_bin > 0) {
use_n <- n_bin * mem
m <- matrix(ordered_cellID[seq_len(use_n)], nrow = n_bin, ncol = mem, byrow = TRUE)
bins <- split(m, row(m))
}
}
}
bins
}
bins <- make_bins(pseudobulkMem)
total_bins <- length(bins)
if (total_bins < min_pseudobulk) {
new_mem <- max(1, floor(length(ordered_cellID) / max(1, min_pseudobulk)))
if (new_mem != pseudobulkMem) {
if (verbose) .ts_log("为满足最小伪样本数，自动调整每个伪样本细胞数为 ", new_mem)
pseudobulkMem <- new_mem
bins <- make_bins(pseudobulkMem)
total_bins <- length(bins)
}
}
if (total_bins < 1) stop("可用细胞数不足以构建 pseudobulk（分组内分箱可能导致不足）")
selected <- as.character(unlist(bins, use.names = FALSE))
pb_names <- paste0("pseudobulk", seq_len(total_bins))
counts_per_bin <- vapply(bins, length, integer(1))
pb_id <- factor(rep(pb_names, times = counts_per_bin), levels = pb_names)
names(pb_id) <- selected
# 记录部分伪样本是否小于 mem
if (verbose) {
n_partial <- sum(counts_per_bin < pseudobulkMem)
if (n_partial > 0) .ts_log("存在 ", n_partial, " 个伪样本细胞数小于 ", pseudobulkMem)
}
# 与矩阵列对齐
selected <- intersect(selected, colnames(rna_counts))
selected <- intersect(selected, colnames(atac_counts))
if (length(selected) == 0) stop("没有匹配的细胞用于聚合")
pb_id <- pb_id[selected]
# 稀疏分配矩阵（n_selected x n_pb）
assign_mat <- Matrix::sparse.model.matrix(~ pb_id - 1)
colnames(assign_mat) <- gsub("^pb_id", "", colnames(assign_mat))
# 2) 聚合计数：genes x cells %*% cells x pseudobulk
if (verbose) .ts_log("聚合 RNA 计数")
rna_pb <- rna_counts[, selected, drop = FALSE] %*% assign_mat
if (verbose) .ts_log("聚合 ATAC 计数")
atac_pb <- atac_counts[, selected, drop = FALSE] %*% assign_mat
colnames(rna_pb) <- colnames(atac_pb) <- colnames(assign_mat)
# 3) 汇总元数据
if (verbose) .ts_log("汇总 pseudobulk 元数据")
meta_sel <- meta[selected, , drop = FALSE]
pb_meta <- .aggregate_meta(meta_sel, groups = pb_id)
pb_meta$pb_id <- rownames(pb_meta)
pb_meta$pb_index <- seq_len(nrow(pb_meta))
# 记录每个伪样本的细胞数量（可用于权重/过滤）
pb_counts <- table(pb_id)
pb_meta$pb_ncells <- as.integer(pb_counts[rownames(pb_meta)])
if (groupBy %in% colnames(meta_sel)) {
suppressWarnings(pb_meta[[paste0("mean_", groupBy)]] <- as.numeric(pb_meta[[groupBy]]))
}
# 对齐列名
pb_meta <- pb_meta[colnames(rna_pb), , drop = FALSE]
# 将非数值型元数据统一转换为字符，避免 list/因子在下游分组时出错
pb_meta <- {
.old <- pb_meta
.new <- lapply(.old, function(col) {
if (is.list(col)) {
vapply(col, function(v) {
if (length(v) == 0) return(NA_character_)
paste(as.character(v), collapse = ";")
}, character(1))
} else if (is.numeric(col) || is.integer(col)) {
col
} else if (is.logical(col)) {
as.character(col)
} else {
as.character(col)
}
})
.df <- as.data.frame(.new, stringsAsFactors = FALSE)
rownames(.df) <- rownames(.old)
.df
}
# 强制确保 groupBy 列为字符，便于 heatmap 分组
if (groupBy %in% colnames(pb_meta)) {
pb_meta[[groupBy]] <- as.character(pb_meta[[groupBy]])
}
# 4) ATAC 注释与基因组信息
if (verbose) .ts_log("准备 ATAC 注释与基因组信息")
annotations <- tryCatch(Annotation(obj[["ATAC"]]), error = function(e) NULL)
genome_info <- NULL
if (!is.null(annotations) && inherits(obj[["ATAC"]], "ChromatinAssay")) {
genome_info <- tryCatch({
gi <- seqinfo(annotations)
if (!is.null(gi)) gi@genome[1] else NA
}, error = function(e) NA)
}
if (is.null(genome_info) || length(genome_info) == 0 || is.na(genome_info)) {
genome_info <- tryCatch({ metadata(annotations)$genome }, error = function(e) NA)
}
if (is.null(genome_info) || length(genome_info) == 0 || is.na(genome_info)) {
genome_info <- NULL
if (verbose) .ts_log("未检测到 genome 信息，将留空")
}
# 5) 构建 Seurat 对象（RNA）并标准化/缩放
if (verbose) .ts_log("构建 Seurat pseudobulk 对象（RNA）并 Normalize/Scale")
pb_seurat <- CreateSeuratObject(counts = rna_pb, assay = "RNA", meta.data = pb_meta)
# 强制确保聚合后的所有元数据列写回（有些版本/路径可能只保留了默认列）
md_existing <- pb_seurat@meta.data
extra_cols <- setdiff(colnames(pb_meta), colnames(md_existing))
if (length(extra_cols) > 0) {
pb_seurat@meta.data <- cbind(
md_existing,
pb_meta[rownames(md_existing), extra_cols, drop = FALSE]
)
}
# 将细胞身份设置为伪样本ID，避免因其它元数据列名造成混淆
try({ Idents(pb_seurat) <- pb_seurat$pb_id }, silent = TRUE)
pb_seurat <- NormalizeData(pb_seurat, assay = "RNA")
pb_seurat <- ScaleData(pb_seurat, assay = "RNA", features = rownames(pb_seurat[["RNA"]]))
# 6) 添加 ATAC assay 并 TF-IDF/Scale
if (verbose) .ts_log("添加 ATAC assay 并运行 TF-IDF/Scale")
pb_seurat[["ATAC"]] <- CreateChromatinAssay(
counts = atac_pb,
annotation = annotations,
genome = genome_info
)
DefaultAssay(pb_seurat) <- "ATAC"
pb_seurat <- RunTFIDF(pb_seurat)
pb_seurat <- ScaleData(pb_seurat, assay = "ATAC", features = rownames(pb_seurat[["ATAC"]]))
if (verbose) .ts_log("pseudobulk 构建完成：列数=", ncol(pb_seurat))
mapping <- split(names(pb_id), f = as.character(pb_id))
return(list(
seurat = pb_seurat,
mapping = mapping,
meta = pb_meta
))
}5. clusterRowsKmeans: K-means 聚类算法封装。 在 Peak-to-Gene 分析中，通常会鉴定出成千上万对显著共变的 Peak-Gene 对。该函数的目的是基于 ATAC 和 RNA 的 Z-score 矩阵，将具有相似跨样本变化模式（如都在某特定亚群中高表达/高开放）的特征聚为一类（共调控模块）。它不仅完成聚类，还通过对各组均值的层次聚类，优化返回的 Cluster 顺序，使得最终展示在热图上的模块视觉效果更为平滑。
R
# 1. 提取的聚类函数
# ==========================================
clusterRowsKmeans <- function(atac_zscores, rna_zscores, cell_groups, k = 10, seed = 1) {
set.seed(seed)
m <- nrow(atac_zscores)
if(k > m) {
message(sprintf("要求的聚类数(k=%d)大于数据点数量(m=%d)", k, m))
k <- 1
message(sprintf("自动调整聚类数为: k=%d", k))
}
k_clusters <- kmeans(atac_zscores, centers = k)
cluster_ids <- k_clusters$cluster
while(1 %in% table(cluster_ids)){
k <- ceiling(k/2)
message(sprintf("聚类数过多，自动调整聚类数为: k=%d", k))
k_clusters <- kmeans(atac_zscores, centers = k)
cluster_ids <- k_clusters$cluster
}
group_means_atac <- matrix(0, nrow = k, ncol = length(levels(as.factor(cell_groups))))
colnames(group_means_atac) <- levels(as.factor(cell_groups))
for (i in 1:k) {
cluster_rows <- which(cluster_ids == i)
if (length(cluster_rows) > 0) {
cluster_atac <- atac_zscores[cluster_rows, , drop = FALSE]
for (g in levels(as.factor(cell_groups))) {
group_cells <- names(cell_groups)[cell_groups == g]
if (length(group_cells) > 0) {
group_means_atac[i, g] <- mean(rowMeans(cluster_atac[, group_cells, drop = FALSE]), na.rm = TRUE)
}
}
}
}
hc <- hclust(dist(group_means_atac))
cluster_order <- hc$order
return(list(
cluster_ids = cluster_ids,
cluster_order = cluster_order,
k = k
))
}6. 下游模块基因的功能通路富集分析 (go_enrich, kegg_enrich)： 当我们把特征聚类成多个共调控模块后，了解每个模块（一堆相关的靶基因）主要参与了哪些生物学功能至关重要。
go_enrich: 调用clusterProfiler包对模块基因进行 GO（Gene Ontology）富集分析，自动输出条形图、气泡图以及详细的表格数据，用于热图右侧的注释。kegg_enrich: 对模块基因进行 KEGG 通路富集分析，并输出相应的图表，揭示代谢或信号通路的富集情况。
R
go_enrich <- function(eg,db,outdir,prefix) {
if (is.null(eg) || nrow(eg) == 0) return(NULL)
key_type <- if ("SYMBOL" %in% colnames(eg)) {
"SYMBOL"
} else if ("ENSEMBL" %in% colnames(eg)) {
"ENSEMBL"
} else if ("ENTREZID" %in% colnames(eg)) {
"ENTREZID"
} else {
return(NULL)
}
genelist <- unique(as.character(eg[[key_type]]))
genelist <- genelist[!is.na(genelist) & genelist != ""]
if (length(genelist) == 0) return(NULL)
go <- enrichGO(genelist, OrgDb=db, ont='ALL',pAdjustMethod = 'BH',qvalueCutoff = 1,pvalueCutoff = 1,keyType = key_type)
go_df <- as.data.frame(go)
if (is.null(go_df) || nrow(go_df) == 0) return(NULL)
go1 <- data.frame(cluster=prefix, go_df)
write.table(go1,file=paste(outdir,'/',prefix,'_GOenrich.xls',sep=''),sep='\t',quote=F,row.names=F)
pdf(paste0(outdir,'/',prefix,'_GO_bar.pdf',sep=''),width = 10, height = 9)
print(barplot(go,showCategory=20,drop=T)+ ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
pdf(paste0(outdir,'/',prefix,'_GO_dot.pdf',sep=''),width = 10, height = 8)
print(dotplot(go,showCategory=20)+ ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
png(paste0(outdir,'/',prefix,'_GO_bar.png',sep=''),width = 780, height = 580)
print(barplot(go,showCategory=20,drop=T)+ ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
png(paste0(outdir,'/',prefix,'_GO_dot.png',sep=''),width = 780, height = 580)
print(dotplot(go,showCategory=20)+ ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
return(go1)
}R
kegg_enrich <- function(eg, kegg_species, outdir, prefix) {
if (is.null(eg) || nrow(eg) == 0) return(NULL)
if (!("ENTREZID" %in% colnames(eg))) return(NULL)
genenames <- if ("SYMBOL" %in% colnames(eg)) as.character(eg$SYMBOL) else as.character(eg$ENTREZID)
names(genenames) <- as.character(eg$ENTREZID)
genelist <- as.character(unique(eg$ENTREZID))
# 定义使用内部KEGG数据的物种（KEGG物种缩写）
internal_species_abbr <- c("hsa", "mmu", "rno", "gga", "ssc")
# 判断是否使用内部数据
if (species %in% internal_species_abbr) {
use_internal <- TRUE
kegg_species <- species
} else {
use_internal <- FALSE
kegg_species <- species
}
kegg <- enrichKEGG(
gene = genelist,
organism = kegg_species,
keyType = "kegg",
pAdjustMethod = "BH",
pvalueCutoff = 1,
qvalueCutoff = 1,
use_internal_data = use_internal
)
kegg_df <- as.data.frame(kegg)
if (is.null(kegg_df) || nrow(kegg_df) == 0) {
message(sprintf("[跳过] %s: KEGG无显著条目", prefix))
return(NULL)
}
gene_list <- strsplit(as.character(kegg_df$geneID), split = "/")
geneName <- vapply(
gene_list,
function(x) paste(genenames[x], collapse = "/"),
character(1)
)
kegg_df$Description <- sub("\\s+-\\s+[^\\(]+\\s*\\([^)]+\\)\\s*$", "",as.character(kegg_df$Description))
KEGGenrich <- data.frame(
cluster = prefix,
kegg_df,
geneName = geneName,
stringsAsFactors = FALSE
)
write.table(
KEGGenrich,
file = paste0(outdir, "/", prefix, "_KEGGenrich.xls"),
sep = "\t",
quote = FALSE,
row.names = FALSE
)
pdf(paste0(outdir, "/", prefix, "_KEGG_dot.pdf"), width = 8, height = 8)
print(dotplot(kegg, showCategory = 20) +
ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
pdf(paste0(outdir, "/", prefix, "_KEGG_bar.pdf"), width = 8, height = 8)
print(barplot(kegg, showCategory = 20) +
ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
png(paste0(outdir, "/", prefix, "_KEGG_dot.png"), width = 780, height = 580)
print(dotplot(kegg, showCategory = 20) +
ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
png(paste0(outdir, "/", prefix, "_KEGG_bar.png"), width = 780, height = 580)
print(barplot(kegg, showCategory = 20) +
ggplot2::scale_y_discrete(labels = function(x) stringr::str_wrap(x, width = 50)))
dev.off()
return(KEGGenrich)
}7. run_Motif_enrichment: 模块 Peak 序列的 Motif 富集函数。 对于每一个聚类簇（Cluster），我们提取其中包含的 Peak 集合，调用 Signac::FindMotifs 进行超几何检验。该函数的作用是回答“在这些具有相似动态变化的染色质开放区域中，哪些转录因子的结合位点出现了显著富集？”。结果将作为热图左侧的核心注释。
R
run_Motif_enrichment <- function(obj, peaks) {
if (is.null(Signac::Motifs(obj[["ATAC"]])) || length(peaks) == 0) return(NULL)
tryCatch({
motif_res <- Signac::FindMotifs(
object = obj,
assay = "ATAC",
features = peaks
)
if (!is.null(motif_res) && nrow(motif_res) > 0) {
return(as.data.frame(motif_res))
}
}, error = function(e) {
message("Motif Error: ", e$message)
NULL
})
return(NULL)
}8. peak2GeneHeatmap: Peak2Gene 分析流程的综合绘图与分析引擎。 这个庞大的函数是本脚本的核心执行者，其目的是：
- 接收
Signac::Links的结果和 Pseudobulk 对象。 - 根据用户参数过滤低置信度的连接。
- 对筛选后的 ATAC 和 RNA 数据进行 Z-score 标准化。
- 调用
clusterRowsKmeans进行共调控模块聚类。 - 对每一个 Cluster，自动化调用富集函数（Motif, GO, KEGG）并保存结果表格。
- 使用
ComplexHeatmap包构建双组学数据、聚类结果以及富集文本的注释对象，并返回热图对象供外部渲染出出版级的高清全景热图。
R
peak2GeneHeatmap <- function(
seurat_obj,
links,
filter = FALSE,
corCutOff = 0.45,
pCutOff = 0.0001,
varCutOffATAC = 0.25,
varCutOffRNA = 0.25,
k = 10,
nPlot = 5000,
groupBy = "seurat_clusters",
palGroup = NULL,
palATAC = ArchR::paletteContinuous("solarExtra"),
palRNA = ArchR::paletteContinuous("blueYellow"),
seed = 1,
run_enrichment = TRUE
) {
set.seed(seed)
peak_matrix <- GetAssayData(seurat_obj, assay = "ATAC", slot = "data")
peak_vars <- MatrixGenerics::rowVars(peak_matrix)
peak_var_quantiles <- rank(peak_vars) / length(peak_vars)
names(peak_var_quantiles) <- rownames(peak_matrix)
rna_matrix <- GetAssayData(seurat_obj, assay = "RNA", slot = "data")
rna_vars <- MatrixGenerics::rowVars(rna_matrix)
rna_var_quantiles <- rank(rna_vars) / length(rna_vars)
names(rna_var_quantiles) <- rownames(rna_matrix)
links_data <- data.frame(
peak = links$peak, gene = links$gene, score = links$score, pvalue = links$pvalue
)
links_data <- subset(links_data, !is.na(peak) & !is.na(gene) & !is.na(score) & !is.na(pvalue))
if (filter){
links_filtered <- links_data[abs(links_data$score) >= corCutOff & links_data$pvalue <= pCutOff, ]
links_filtered$peak_id <- match(links_filtered$peak, rownames(peak_matrix))
links_filtered$gene_id <- match(links_filtered$gene, rownames(rna_matrix))
links_filtered$VarQATAC <- peak_var_quantiles[rownames(peak_matrix)[links_filtered$peak_id]]
links_filtered$VarQRNA <- rna_var_quantiles[rownames(rna_matrix)[links_filtered$gene_id]]
links_filtered <- links_filtered[links_filtered$VarQATAC > varCutOffATAC, ]
links_filtered <- links_filtered[links_filtered$VarQRNA > varCutOffRNA, ]
} else {
links_filtered <- links_data
# 即使不过滤，也必须强制剔除方差为 0 的行，否则后续 scale() 会产生 NA 导致聚类崩溃
links_filtered$peak_id <- match(links_filtered$peak, rownames(peak_matrix))
links_filtered$gene_id <- match(links_filtered$gene, rownames(rna_matrix))
# 提取对应的方差真实值
links_filtered$VarATAC <- peak_vars[links_filtered$peak_id]
links_filtered$VarRNA <- rna_vars[links_filtered$gene_id]
# 剔除方差为 0 的特征
links_filtered <- links_filtered[links_filtered$VarATAC > 0 & links_filtered$VarRNA > 0, ]
# 清理临时列，保持和原逻辑一致
links_filtered$VarATAC <- NULL
links_filtered$VarRNA <- NULL
}
if (nrow(links_filtered) == 0) stop("过滤后links为空，无法进行聚类与富集")
peak_ids <- match(links_filtered$peak, rownames(peak_matrix))
gene_ids <- match(links_filtered$gene, rownames(rna_matrix))
atac_data <- peak_matrix[peak_ids, ]
rna_data <- rna_matrix[gene_ids, ]
rowZscores <- function(matrix1, matrix2, min = NULL, max = NULL, limit = TRUE) {
z1 <- t(scale(t(matrix1)))
z2 <- t(scale(t(matrix2)))
if(limit) {
if(is.null(min) || is.null(max)) {
middle_range1 <- quantile(abs(z1), probs = 0.875, na.rm = TRUE)
if (middle_range1 < 1) { min1 <- -1; max1 <- 1
} else if (middle_range1 < 1.5) { min1 <- -1.5; max1 <- 1.5
} else { min1 <- -2; max1 <- 2 }
middle_range2 <- quantile(abs(z2), probs = 0.875, na.rm = TRUE)
if (middle_range2 < 1) { min2 <- -1; max2 <- 1
} else if (middle_range2 < 1.5) { min2 <- -1.5; max2 <- 1.5
} else { min2 <- -2; max2 <- 2 }
min <- max(min1, min2)
max <- min(max1, max2)
}
z1[z1 > max] <- max; z1[z1 < min] <- min
z2[z2 > max] <- max; z2[z2 < min] <- min
}
return(list(matrix1 = z1, matrix2 = z2))
}
zs <- rowZscores(as.matrix(atac_data), as.matrix(rna_data))
atac_zscores <- zs$matrix1
rna_zscores <- zs$matrix2
cell_groups <- seurat_obj@meta.data[, groupBy]
names(cell_groups) <- colnames(seurat_obj)
v <- seurat_obj@meta.data[[groupBy]]
sorted_levels <- stringr::str_sort(unique(as.character(v)), numeric = TRUE, na_last = TRUE)
cell_groups <- factor(v, levels = sorted_levels, ordered = TRUE)
# 1. 独立调用聚类函数
clust_res <- clusterRowsKmeans(atac_zscores, rna_zscores, cell_groups, k = k, seed = seed)
cluster_ids <- clust_res$cluster_ids
cluster_order <- clust_res$cluster_order
k <- clust_res$k
ordered_links <- list()
ordered_atac <- list()
ordered_rna <- list()
cluster_info_list <- list()
# 行分组标签（用于热图分割和注释）
row_split_factor <- c()
for (i in seq_along(cluster_order)) {
orig_cluster_id <- cluster_order[i]
cluster_rows <- which(cluster_ids == orig_cluster_id)
if (length(cluster_rows) > 0) {
ordered_links <- c(ordered_links, list(links_filtered[cluster_rows, ]))
ordered_atac <- c(ordered_atac, list(atac_zscores[cluster_rows, ]))
ordered_rna <- c(ordered_rna, list(rna_zscores[cluster_rows, ]))
cluster_links <- links_filtered[cluster_rows, ]
cluster_links$kmeans_cluster <- i
cluster_info_list[[i]] <- cluster_links
row_split_factor <- c(row_split_factor, rep(i, length(cluster_rows)))
}
}
ordered_links_df <- do.call(rbind, ordered_links)
ordered_atac_mat <- do.call(rbind, ordered_atac)
ordered_rna_mat <- do.call(rbind, ordered_rna)
cluster_info_df <- do.call(rbind, cluster_info_list)
# 保存聚类结果
write.table(cluster_info_df, file = file.path(dirname(outdir), "result", "Peak2Gene_Kmeans_Clusters.tsv"),
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
# ----------------------------------------------------
# 2. 执行富集并构建行注释 (Motif & GO)
# ----------------------------------------------------
atac_row_anno_text <- rep("", k)
rna_row_anno_text <- rep("", k)
if (run_enrichment) {
enrich_outdir <- file.path(dirname(outdir), "result", "Cluster_Enrichment")
dir.create(enrich_outdir, showWarnings = FALSE, recursive = TRUE)
has_motif <- !is.null(Signac::Motifs(obj[["ATAC"]]))
if (!has_motif) {
.ts_log("ATAC assay 未检测到 Motifs；若需进行 Motif 富集，请先添加 Motif 注释")
}
all_motif_res <- list()
all_go_res <- list()
all_kegg_res <- list()
for (cid in seq_len(k)) {
sub_df <- cluster_info_df[cluster_info_df$kmeans_cluster == cid, ]
if(nrow(sub_df) == 0) next
.ts_log("正在对 Cluster ", cid, " 进行富集分析计算用于注释...")
motif_res <- run_Motif_enrichment(obj, unique(sub_df$peak))
if (!is.null(motif_res)) {
motif_res$cluster <- as.character(cid)
all_motif_res[[as.character(cid)]] <- motif_res
write.table(motif_res, file.path(enrich_outdir, paste0("Cluster_", cid, "_Motif.tsv")), sep="\t", row.names=F, quote=F)
top6_motif <- head(motif_res[order(motif_res$pvalue), "motif.name"], 6)
top6_motif <- gsub("\\s*\\(.*\\)", "", top6_motif)
if(length(top6_motif) > 0) {
line1 <- paste(top6_motif[seq_len(min(3, length(top6_motif)))], collapse = " ")
line2 <- if (length(top6_motif) > 3) paste(top6_motif[4:min(6, length(top6_motif))], collapse = " ") else NULL
atac_row_anno_text[cid] <- paste(c(line1, line2), collapse = "\n")
}
} else {
.ts_log("Cluster ", cid, " 未富集到显著 Motif")
}
genes_raw <- unique(as.character(sub_df$gene))
genes_raw <- genes_raw[!is.na(genes_raw) & genes_raw != ""]
genes_clean <- genes_raw
.ts_log("Cluster ", cid, " 输入基因数: ", length(genes_clean))
eg <- NULL
used_from <- "SYMBOL"
eg_try <- tryCatch({
clusterProfiler::bitr(
genes_clean,
fromType = used_from,
toType = unique(c("SYMBOL", "ENSEMBL", "ENTREZID")),
OrgDb = get(db)
)
}, error = function(e) NULL)
if (!is.null(eg_try) && nrow(eg_try) > 0) {
eg <- unique(eg_try)
}
if (!is.null(eg) && nrow(eg) > 0) {
keep_cols <- intersect(c("SYMBOL", "ENSEMBL", "ENTREZID"), colnames(eg))
eg <- unique(eg[, keep_cols, drop = FALSE])
.ts_log("Cluster ", cid, " ID映射成功: ", nrow(eg), " 条，fromType=", used_from)
} else {
.ts_log("Cluster ", cid, " 基因ID映射失败，输入的 SYMBOL 未能匹配到数据库")
}
go_res <- go_enrich(
eg = eg,
db = get(db),
outdir = enrich_outdir,
prefix = paste0("Cluster_", cid)
)
if (!is.null(go_res)) {
go_res_df <- as.data.frame(go_res)
if (nrow(go_res_df) > 0) {
go_res_df$cluster <- as.character(cid)
all_go_res[[as.character(cid)]] <- go_res_df
}
top5_go <- head(go_res_df[order(go_res_df$p.adjust), "Description"], 5)
if (length(top5_go) > 0) {
go_text <- stringr::str_wrap(paste(top5_go, collapse = "; "), width = 100)
go_lines <- strsplit(go_text, "\n", fixed = TRUE)[[1]]
rna_row_anno_text[cid] <- paste(head(go_lines, 2), collapse = "\n")
}
} else {
.ts_log("Cluster ", cid, " 未富集到显著 GO Term")
}
# 使用 tryCatch 捕获网络报错，防止程序崩溃
kegg_res <- tryCatch({
kegg_enrich(
eg = eg,
kegg_species = kegg_species,
outdir = enrich_outdir,
prefix = paste0("Cluster_", cid)
)
}, error = function(e) {
.ts_log("Cluster ", cid, " KEGG富集请求失败 (网络原因): ", e$message)
return(NULL)
})
if (!is.null(kegg_res)) {
# 如果有需要处理 kegg_res 的逻辑可以写在这里
kegg_res_df <- as.data.frame(kegg_res)
if (nrow(kegg_res_df) > 0) {
kegg_res_df$cluster <- as.character(cid)
all_kegg_res[[as.character(cid)]] <- kegg_res_df
} else {
.ts_log("Cluster ", cid, " 未富集到显著 KEGG Pathway 或因网络跳过")
}
}
# --- 汇总所有 Cluster 的 Motif 结果并保存 ---
if (length(all_motif_res) > 0) {
combined_motif <- do.call(rbind, all_motif_res)
write.table(combined_motif, file.path(enrich_outdir, "All_Cluster_Motif.tsv"), sep="\t", row.names=F, quote=F)
}
# --- 汇总所有 Cluster 的 GO 结果并保存 ---
if (length(all_go_res) > 0) {
combined_go <- do.call(rbind, all_go_res)
write.table(combined_go, file.path(enrich_outdir, "All_Cluster_GOenrich.xls"), sep="\t", row.names=F, quote=F)
.ts_log("所有 Cluster 的 GO 结果已合并保存为 All_Cluster_GOenrich.xls。")
}
# --- 汇总所有 Cluster 的 KEGG 结果并保存 ---
if (length(all_kegg_res) > 0) {
combined_kegg <- do.call(rbind, all_kegg_res)
write.table(combined_kegg, file.path(enrich_outdir, "All_Cluster_KEGGenrich.xls"), sep="\t", row.names=F, quote=F)
.ts_log("所有 Cluster 的 KEGG 结果已合并保存为 All_Cluster_KEGGenrich.xls。")
}
}
if(nrow(cluster_info_df) < 20) stop("link过少，建议放宽过滤阈值")
plot_idx <- seq_len(nrow(cluster_info_df))
if (nrow(cluster_info_df) > nPlot) {
# 为了保证聚类的相对顺序不被打乱（这会导致 row_split 和注释错位）
# 必须对抽样后的 idx 进行排序
plot_idx <- sort(sample(plot_idx, nPlot))
}
links_filtered <- ordered_links_df[plot_idx, , drop = FALSE]
ordered_atac_mat <- ordered_atac_mat[plot_idx, , drop = FALSE]
ordered_rna_mat <- ordered_rna_mat[plot_idx, , drop = FALSE]
# 注意：在对矩阵进行 downsample 时，row_split_factor 也必须对应地被截取，
# 但之前计算出的 mark_at 依赖于被截取后的 row_split_factor
row_split_factor <- row_split_factor[plot_idx]
# 准备分组颜色
if (is.null(palGroup)) {
group_levels <- levels(as.factor(cell_groups))
palGroup <- setNames(scales::hue_pal()(length(group_levels)), group_levels)
}
# 动态计算图例行数，防止类别过多时图例溢出或遮挡
n_groups <- length(unique(cell_groups))
legend_nrow <- ifelse(n_groups > 8, ceiling(n_groups / 8), 1)
column_atac <- ComplexHeatmap::HeatmapAnnotation(Group = cell_groups, col = list(Group = palGroup), show_legend = FALSE, show_annotation_name = FALSE)
column_rna <- ComplexHeatmap::HeatmapAnnotation(
Cell = cell_groups,
col = list(Cell = palGroup),
show_annotation_name = FALSE,
annotation_legend_param = list(
Cell = list(nrow = legend_nrow, title_position = "topcenter") # 动态控制离散型图例行数
)
)
# ----------------------------------------------------
# 3. 统一使用 anno_mark 构建左右两侧注释（完全废弃 row_split）
# ----------------------------------------------------
split_fac <- factor(row_split_factor, levels = seq_len(k))
# 预先计算每个 Cluster 在绝对行坐标系下的中间位置
cluster_mid_row <- vapply(seq_len(k), function(cid) {
idx <- which(split_fac == cid)
if (length(idx) == 0) NA_integer_ else as.integer(round(mean(range(idx))))
}, integer(1))
# --- 3.1 构建左侧行注释 (Motif) ---
valid_cid_atac <- which(!is.na(cluster_mid_row) & nzchar(atac_row_anno_text))
mark_at_atac <- cluster_mid_row[valid_cid_atac]
mark_labels_atac <- atac_row_anno_text[valid_cid_atac]
row_anno_atac <- ComplexHeatmap::rowAnnotation(
Motif = ComplexHeatmap::anno_mark(
at = mark_at_atac,
labels = mark_labels_atac,
which = "row",
side = "left", # 关键：由于在左侧，引线必须向左画
labels_gp = grid::gpar(fontsize = 10, fontface = "bold", col = "#0c0c0cff"),
lines_gp = grid::gpar(lty = 2, col = "black"),
link_width = unit(5, "mm"),
padding = unit(1.5, "mm"),
extend = unit(c(1, 10), "mm")
),
width = unit(1, "cm")
)
# --- 3.2 构建右侧行注释 (GO Enrichment) ---
valid_cid_rna <- which(!is.na(cluster_mid_row) & nzchar(rna_row_anno_text))
mark_at_rna <- cluster_mid_row[valid_cid_rna]
mark_labels_rna <- stringr::str_wrap(rna_row_anno_text[valid_cid_rna], width = 65)
row_anno_rna <- ComplexHeatmap::rowAnnotation(
GO = ComplexHeatmap::anno_mark(
at = mark_at_rna,
labels = mark_labels_rna,
which = "row",
side = "right", # 引线向右画
labels_gp = grid::gpar(fontsize = 10, fontface = "bold", col = "#0c0c0cff"),
lines_gp = grid::gpar(lty = 2, col = "black"),
link_width = unit(8, "mm"),
padding = unit(0.005, "mm"),
extend = unit(c(0.001, 1), "mm")
),
width = unit(1, "cm")
)
panel_width <- unit(max(10, ncol(ordered_rna_mat) * 0.02), "cm")
# ----------------------------------------------------
# 4. 创建并绘制热图
# ----------------------------------------------------
atac_heatmap <- ComplexHeatmap::Heatmap(
ordered_atac_mat,
name = "ATAC Z-score",
col = palATAC,
cluster_rows = FALSE, cluster_columns = FALSE,
column_order = order(colMeans(ordered_rna_mat), decreasing = TRUE),
show_row_names = FALSE, show_column_names = FALSE,
top_annotation = column_atac,
left_annotation = row_anno_atac,
column_split = cell_groups, column_gap = unit(0, "mm"), column_title = NULL,
# 完全注释掉热图物理切割功能，维持全局绝对坐标系
# row_split = factor(row_split_factor, levels=1:k),
# row_title_rot = 0,
row_order = seq_len(nrow(ordered_atac_mat)),
cluster_row_slices = FALSE,
width = panel_width,
use_raster = TRUE, raster_quality = 1,
heatmap_legend_param = list(direction = "horizontal", title_position = "topcenter")
)
rna_heatmap <- ComplexHeatmap::Heatmap(
ordered_rna_mat,
name = "RNA Z-score",
col = palRNA,
cluster_rows = FALSE, cluster_columns = FALSE,
column_order = order(colMeans(ordered_rna_mat), decreasing = TRUE),
show_row_names = FALSE, show_column_names = FALSE,
top_annotation = column_rna,
right_annotation = row_anno_rna,
column_split = cell_groups, column_gap = unit(0, "mm"), column_title = NULL,
# 同样注释掉物理切割功能
# row_split = factor(row_split_factor, levels=1:k), row_title = NULL,
width = panel_width,
use_raster = TRUE, raster_quality = 1,
heatmap_legend_param = list(direction = "horizontal", title_position = "topcenter")
)
return(list(
links = cluster_info_df, #返回的 cluster_info_df 是经过了前面的 NA 去除、阈值过滤（或零方差过滤）、以及 K-means 聚类排序后的 全量数据
atac_heatmap = atac_heatmap,
rna_heatmap = rna_heatmap
))
}
}9. plotPeaksPerGene: 基础统计图函数（基因受控度分布）。 一个基因可能被周围的多个 Enhancer（Peak）共同调控。该函数绘制直方图，直观展示全局网络中“每个基因平均被多少个 Peak 关联”，帮助评估调控网络的复杂度和连通性。
R
plotPeaksPerGene <- function(links) {
# 统计每个基因连接的peak数量
gene_counts <- table(links$gene)
gene_counts_df <- data.frame(
gene = names(gene_counts),
count = as.numeric(gene_counts)
)
# 计算中位数
median_linked_Peaks <- median(gene_counts_df$count)
# 绘制直方图
p <- ggplot2::ggplot(gene_counts_df, ggplot2::aes(x = count)) +
ggplot2::geom_histogram(binwidth = 1, fill = "black", color = "white") +
ggplot2::geom_vline(xintercept = median_linked_Peaks, linetype = "dashed", color = "gray") +
ggplot2::annotate("text", x = 10, y = 1000,
label = paste0("Median Linked Peaks = ", median_linked_Peaks)) +
ggplot2::labs(x = "Number of linked peaks", y = "Number of genes") +
ggplot2::theme_classic() +
ggplot2::xlim(0, 25) #+
# theme(text = element_text(family = "sans"))
return(p)
}10. plotPeak2GeneVolcano: 基础统计图函数（关联显著性火山图）。 类似于差异表达分析的火山图，这里展示的是所有 Peak-Gene 对的 Correlation（X轴）与 P-value（Y轴）。帮助学者一眼识别出那些既具有极强相关性、又具有极高统计学显著性的“Top 核心驱动基因”。
R
# 添加火山图：展示峰区域与基因连接的相关性和显著性
plotPeak2GeneVolcano <- function(links, topN = 20, p = 0.05, score = 0.05) {
# 准备数据
volcano_data <- data.frame(
peak = links$peak,
gene = links$gene,
correlation = links$score,
pvalue = links$pvalue,
logPvalue = -log10(links$pvalue + 1e-300) # 避免log(0)
)
# 直接找出最显著的topN个基因
top_genes_indices <- order(volcano_data$logPvalue, decreasing = TRUE)[1:min(topN, nrow(volcano_data))]
top_genes <- volcano_data[top_genes_indices, ]
# 标记哪些点需要高亮显示（被标注的点）
volcano_data$highlighted <- "No"
volcano_data$highlighted[top_genes_indices] <- "Yes"
top_genes$highlighted <- "Yes"
# 绘制火山图
p <- ggplot2::ggplot(volcano_data, ggplot2::aes(x = correlation, y = logPvalue, color = highlighted)) +
ggplot2::geom_point(alpha = 0.6, size = 1) +
ggplot2::scale_color_manual(values = c("No" = "gray", "Yes" = "red")) +
ggplot2::labs(
# title = "Peak-to-Gene Linkage Volcano Plot",
x = "Correlation Score",
y = "-log10(p-value)",
color = "Top Genes"
) +
ggplot2::theme_classic() +
ggplot2::theme(
legend.position = "none",
# text = element_text(family = "sans")
) +
# 添加虚线
ggplot2::geom_vline(xintercept = -score, linetype = "dashed", color = "black", alpha = 0.7) +
ggplot2::geom_vline(xintercept = score, linetype = "dashed", color = "black", alpha = 0.7) +
ggplot2::geom_hline(yintercept = -log10(p), linetype = "dashed", color = "black", alpha = 0.7)
# 标记顶部基因
if(nrow(top_genes) > 0) {
p <- p + ggrepel::geom_text_repel(
data = top_genes,
# family = "sans",
ggplot2::aes(label = gene),
size = 3,
box.padding = 0.5,
point.padding = 0.2,
force = 2,
max.overlaps = 20
)
}
return(p)
}11. plotRankCorrelation: 基础统计图函数（相关性排序图）。 将所有关联对的相关系数从大到小排序并绘制曲线。可以用于评估整体网络中强相关与弱相关对的比例分布趋势。
R
# 使用links中的所有score值绘制相关性排序图
plotRankCorrelation <- function(links) {
# 提取所有相关性值
all_correlations <- links$score
# 按相关性从高到低排序
sorted_correlations <- sort(all_correlations, decreasing = TRUE)
correlation_df <- data.frame(
rank = 1:length(sorted_correlations),
correlation = sorted_correlations
)
# 绘制排序图
p <- ggplot2::ggplot(correlation_df, ggplot2::aes(x = rank, y = correlation)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::geom_hline(yintercept = 0, linetype = "solid", color = "black") +
ggplot2::labs(
# title = "Peak-Gene correlation",
x = "Peak-Gene Linkage",
y = "Correlation Score"
) +
ggplot2::theme_classic() +
ggplot2::theme(
# 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轴刻度线
)
return(p)
}4. Peak-to-Gene 关联
Peak-to-Gene 分析是本教程的核心步骤。由于增强子（Enhancer）等调控元件可能通过染色质空间折叠，跨越很长的线性距离来调控靶基因的表达，单纯依赖线性距离（如 promoter 区域）往往无法准确捕获这种调控关系。在这里，我们首先检查对象中是否已存在关联信息；如果不存在，则先计算 GC 含量等区域统计信息 (RegionStats)，再使用 Signac 提供的 LinkPeaks 函数，计算每一个染色质可及性 Peak 与周围基因表达量之间的统计学共变性（Correlation）。如果同一细胞中某个 Peak 开放度的变化始终伴随着特定基因表达量的同步变化，那么我们可以假设这两者之间存在顺式调控关系。
R
# 检查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
)
}))
}))R
saveRDS(obj, file.path(outdir, "output.rds"))5. 构建 Pseudobulk (假混合) 矩阵
在获得了基础的 Peak-to-Gene 链接后，直接利用稀疏的单细胞矩阵去计算组间（如各 Cluster）的精细动态差异往往会受到严重的 Drop-out（测序信号丢失）干扰。 为了提高信噪比、提升相关性计算的稳健性并大幅降低后续聚类热图的计算开销，我们通过 build_pseudobulk_seurat 函数，按照细胞群标识（clusters_col）将同属于一类的数百个细胞的 Count 矩阵进行累加（Aggregation），生成低维度、高信噪比的假混合（Pseudobulk）对象。
R
invisible(capture.output({
suppressMessages(suppressWarnings({
res <- build_pseudobulk_seurat(obj = obj, groupBy = clusters_col)
pseudo_obj <- res$seurat
}))
}))6. 共调控模块 K-means 聚类与富集分析
这是流程中最为复杂和关键的数据整合与可视化步骤，由我们封装的 peak2GeneHeatmap 函数自动完成。它实现了以下完整的流水线：
- 数据过滤与缩放：根据设定的阈值（
corCutOff,pCutOff）筛除低置信度的关联对，并将 Pseudobulk 矩阵中的表达和可及性数据转换为标准化的 Z-score 矩阵。 - K-means 模块聚类：将具有类似动态变化趋势（即在相同的细胞群中同时上调或下调）的 Peak-Gene 组合，划分为
k个独立的共调控簇（Cluster）。 - 功能富集推断：
- 上游调控追溯 (Motif)：对每个簇的 Peak 序列执行转录因子 Motif 富集分析，推测是哪些 TF 驱动了该模块的开放。
- 下游功能揭示 (GO/KEGG)：对每个簇关联的靶基因执行通路富集，解释这些基因主要参与了哪些生物学过程。
- 联合热图对象构建：使用
ComplexHeatmap构建 ATAC 信号、RNA 信号、Motif 注释（左侧）、GO 注释（右侧）的热图对象并返回，供外部渲染为一张直观的多组学调控网络全景图。
R
result <- peak2GeneHeatmap(
pseudo_obj,
links = links_all,
filter = filter,
corCutOff = corCutOff,
pCutOff = pCutOff,
varCutOffATAC = varCutOffATAC,
varCutOffRNA = varCutOffRNA,
k = k,
nPlot = nPlot,
groupBy = clusters_col
)R
# ComplexHeatmap 需要 draw 时计算实际尺寸，为了确保注释显示，可以尝试在 draw 里显式指定
options(repr.plot.width = 20, repr.plot.height = 10)
draw((result$atac_heatmap + result$rna_heatmap),
heatmap_legend_side = "bottom",
annotation_legend_side = "bottom",
ht_gap = unit(1, "mm"),
auto_adjust = TRUE,
padding = unit(c(4, 3, 20, 15), "mm"))
💡 Peak-to-Gene 共调控模块联合热图
该热图展示了染色质开放区域与基因表达之间的相关性模式。热图中的每一行代表一个 Peak-Gene 对，每一列代表一个细胞，且整体行数据已根据共调控模式进行了 K-means 聚类划分。
- 顶部颜色条：标识了不同的细胞类型或群体，帮助直观比对特定细胞亚群中 Peak-Gene 关联的活跃模式。
- ATAC 热图（主图左侧）：显示 ATAC-seq 信号强度（Z-score 归一化可及性）。颜色越蓝表示可及性越低，颜色越红表示可及性越高。
- RNA 热图（主图右侧）：显示相应靶基因的表达水平（Z-score 归一化表达量）。颜色越蓝表示表达量越低，颜色越黄表示表达量越高。
- Motif 注释（最左侧文本）：标识了每个 K-means 聚类簇中显著富集的转录因子结合基序（Top Motifs）。这些 TF 往往是驱动该簇基因表达的上游核心调控因子。
- GO 功能注释（最右侧文本）：标识了每个 K-means 聚类簇中靶基因显著富集的基因本体论条目（Top GO Terms），直观揭示了该共调控模块所参与的核心生物学功能或细胞通路。
7. 结果展示与统计
除了最核心的聚类热图外，您还可以在指定的 outdir 目录中找到所有的富集结果表格和连线详细信息表。 此外，为了从不同维度评估本次分析的全局质量，您还可以调用脚本中预定义的以下几个作图函数：
plotPeaksPerGene(links_all)：绘制直方图，展示每个基因平均被多少个 Peak 所调控，评估调控网络的连通度。plotPeak2GeneVolcano(links_all)：绘制火山图，直观展现显著性（P-value）与相关性（Correlation）的分布，突出核心驱动基因。plotRankCorrelation(links_all)：绘制所有关联对的相关性打分排序图，评估整体关联分布趋势。
R
options(repr.plot.width = 10, repr.plot.height = 5)
plotPeaksPerGene(result$links)
💡 基因-Peak 关联数量分布直方图
该直方图主要用于展示靶基因与其显著关联的调控元件（Peak）数量的群体分布情况（默认截取展示 1 ~ 25 个关联 Peak 的区间）。
- 横轴（X轴）：代表单个基因所关联的 Peak 数量。
- 纵轴（Y轴）：代表拥有对应数量关联 Peak 的基因总数。柱子越高，说明拥有该数量调控元件的基因群体越庞大。
- 虚线与标注：图中的垂直虚线指示了所有靶基因关联 Peak 数量的中位数所在位置，旁边的文本直接标注了中位数的具体数值。这有助于快速评估全局顺式调控网络的复杂程度底线。
R
if (filter == TRUE) {
print(plotPeak2GeneVolcano(result$links, p = pCutOff, score = corCutOff))
} else {
print(plotPeak2GeneVolcano(result$links))
}
💡 Peak-to-Gene 关联火山图
该火山图用于直观展示所有 Peak-Gene 关联对的相关性强度与统计学显著性，帮助快速定位最关键的调控关系。
- 横轴（X轴）：代表相关性得分。它表示 Peak 的可及性与靶基因表达量之间的相关程度，数值绝对值越大说明共变性越强。
- 纵轴（Y轴）：代表
-log10(p-value)。数值越高表示该 Peak-Gene 关联的统计学显著性越强、假阳性率越低。- 文本标注：图中标注了 Top 10 最显著的 Peak-Gene 关联对应的基因名称，这些基因通常受极其强烈的顺式元件调控。
R
print(plotRankCorrelation(result$links))
💡 Peak-to-Gene 相关性排序分布图 (Rank Correlation Plot)
该折线图全景式地展示了所有检测到的 Peak-Gene 关联对的相关性强度分布，用于直观评估全局顺式调控网络的统计学强度和置信度阈值。
- 横轴（X轴，Rank）：代表所有 Peak-Gene 关联对的排序（Rank）。系统将所有的调控对按照其“相关性得分（Correlation Score）”从高到低进行降序排列，排名越靠前（X 轴数值越小）代表相关性越强。
- 纵轴（Y轴，Correlation）：代表每个 Peak-Gene 关联对实际的 Pearson/Spearman 相关系数得分。正值表示该 Peak 的开放伴随着靶基因表达的上调（经典的增强子激活作用），负值则可能暗示抑制作用。
- 虚线与标注（阈值截断）：图中的水平虚线通常表示用于后续下游分析的“强关联置信度截断阈值（Cutoff）”。通过观察曲线的拐点（Knee Point）与虚线的交汇处，可以快速评估在当前过滤标准下，保留了多少具有强劲生物学意义的高质量调控对，以及整体背景噪声的水平。