ATAC + RNA 多组学：Motif 分析
时长: 67 分钟
字数: 15.7k 字
更新: 2026-07-07
阅读: 0 次
1. 教程简介
本教程基于 Seurat、Signac 以及 JASPAR2020 等单细胞多组学分析工具，专门用于对 SeekArc 单细胞多组学数据（scRNA-seq + scATAC-seq） 进行深入的表观遗传调控网络探索。
该教程面向已完成基础质控与细胞群注释的 SeekArc 多组学数据，通过识别细胞群或不同状态下的特异性开放区域，结合底层基因调控机制，实现以下三大核心分析目标：
- Motif 活性分析 (Motif Activity)：利用 ChromVAR 算法，计算单细胞全基因组范围内针对各个转录因子 Motif 的全局开放程度评分，从而评估不同转录因子的潜在整体活性。
- Motif 活性差异分析 (Differential Motif Activity)：通过比较不同细胞群（Cluster vs All）或同一细胞群在不同分组下（Case vs Control）的 Motif 活性得分，精准锁定在特定细胞类型或特定状态下异常活跃的关键转录因子。
- Motif 富集分析 (Motif Enrichment)：基于前期鉴定出的差异 Peak 集合（特异性染色质开放区域），利用超几何检验等统计方法寻找显著富集的转录因子结合位点（TFBS），以此识别驱动这些特异性开放区域的核心调控因子。
注意
💡 文档中展示的图表仅为部分代表性的可视化结果。如需查看所有细胞群或对比组的完整分析结果（包含详细的数据表格与所有高清绘图），请前往
./result/目录进行查阅。
R
suppressPackageStartupMessages(suppressWarnings({
library(presto)
library(ggrepel)
library(Seurat)
library(BiocGenerics)
library(S4Vectors)
library(IRanges)
library(Signac)
library(ComplexHeatmap)
library(Biobase)
library(AnnotationDbi)
#library(org.Hs.eg.dbb) #人类基因注释数据库
library(org.Mm.eg.db) #小鼠基因注释数据库
library(GenomicRanges)
library(JASPAR2020)
library(TFBSTools)
library(memuse)
library(stringr)
library(dplyr)
library(DOSE)
library(clusterProfiler)
library(enrichplot)
library(foreach)
library(doParallel)
library(base64enc)
library(DT)
library(KEGG.db)
}))R
# 初始化并行计算框架
if (requireNamespace("future", quietly = TRUE)) {
cores <- parallelly::availableCores()
workers <- as.integer(floor(cores * 0.7))
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.7
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%)")output
Footprint 并行计算已启动，总核心数=32，使用 22 个工作线程
2. 输入文件准备
2.1 输入文件要求
本教程主要基于已经初步处理好的 Seurat 多组学对象进行 Motif 分析。开始分析前，需要准备以下几类核心输入文件：
input.rds：预处理完成的 Seurat 对象文件。该对象必须包含RNA和ATAC两个 Assay，并且 ATAC 模态中需要已经完成基本的 Peak 调用。meta.tsv：细胞元数据文件（制表符分隔）。必须包含barcode列，以及用于定义细胞类型（如Celltype）和样本分组（如Sample）的对应列。genome.fa：对应物种的参考基因组 FASTA 文件。用于后续在 ATAC peaks 中提取 DNA 序列，以进行 Motif 匹配和 ChromVAR TF 活性打分。ATAC 片段文件及索引：构建 ATAC Assay 所依赖的fragments.tsv.gz及其.tbi索引文件。请确保这些文件存放在指定的输出目录（或与脚本中配置的路径保持一致），以便 Signac 能够正确读取。
目录结构示例
建议按照以下结构组织你的输入数据和元数据：
text
PBMC_demo/
├── input.rds # 包含 scRNA 和 scATAC 数据的 Seurat 对象
├── meta.tsv # 包含细胞类型、分组信息的元数据表
├── genome.fa # 参考基因组序列文件 (如 hg38.fa)
└── fragments/ # ATAC 片段文件存放位置
├── sample1_fragments.tsv.gz
└── sample1_fragments.tsv.gz.tbiR
# --- 输入参数配置 ---此处参数按照自己需求填写参数
## fpath：fragment文件所在的目录
fpath = "/path/to/PBMC_demo"
## rds：输入的 Seurat 对象 RDS 文件路径，包含已预处理的单细胞多组学数据
rds = "/path/to/input.rds"
## meta：元数据文件路径（TSV格式），包含样本信息、细胞类型注释等元数据
meta = "/path/to/meta.tsv"
## species：物种信息，"others"表示非标准模式物种（如人/小鼠以外的物种）
species = "human"
## ref_genome：参考基因组 FASTA 文件路径，用于基因组序列比对和注释
ref_genome = "/path/to/genome.fa"
## clusters_col：元数据中用于标识细胞类型/聚类的列名
clusters_col = "Celltype"
## celltypes：需要分析的细胞类型列表，多个细胞类型用逗号分隔
celltypes = "B cells,CMP,Dividing B cells,Erythroblast,NK cells,pDC,Plasma Cells,Pro B cells,T cells"
## downsample：是否进行细胞下采样，TRUE表示启用下采样以平衡各细胞类型数量
downsample = "TRUE"
## downsample_num：下采样后的目标细胞数，每个细胞类型最多保留的细胞数量
downsample_num = "3000"
## group_comparison：元数据中用于分组比较的列名，通常为样本分组或处理条件
group_comparison = "Sample"
## case_group：实验组/病例组的分组名称，用于差异富集分析
case_group = "XYRD_pbmc_2_arc"
## control_group：对照组的分组名称，用于差异富集分析
control_group = "25030508_pbmc_1_arc"
## logFC：对数倍变化（log fold change）阈值，用于筛选显著差异的特征
logFC = "0.25"
## pval_adj：调整后 p 值阈值，用于统计显著性检验
pval_adj = "0.05"
## cor：相关性阈值，用于特征间相关性分析
cor = "0.6"
## fe_cut：富集倍数（fold enrichment）阈值，用于筛选显著富集的区域
fe_cut = "1.3"R
# 参数处理，无需改动
group_colors <- c(
"#F7E55C",
"#F08080"
)
names(group_colors) <- c(case_group, control_group)
pseudobulk_size <- 5
output <- "./result"
celltypes <- strsplit(celltypes,",")[[1]]
downsample_num <- as.numeric(downsample_num)
logFC <- as.numeric(logFC)
pval_adj <- as.numeric(pval_adj)R
#创建必要的结果文件夹无需改动
dir.create(paste0(output,'/DIFF/01_diffPeak'), recursive = TRUE)
dir.create(paste0(output,'/DIFF/02_motif'), recursive = TRUE)2.2 核心辅助函数定义
为了提高代码的模块化和复用性，我们在主流程开始前预先定义了以下几个核心函数。这些函数涵盖了数据聚合、名称标准化、统计计算以及图表排版，它们将在后续的不同分析阶段被频繁调用。
make_pseudobulk- 功能：将高度稀疏的单细胞矩阵按指定的细胞分组和合并大小（
pb_size）压缩为拟批量（pseudobulk）矩阵。 - 目的：克服单细胞 ATAC/RNA 数据的稀疏性（Dropout），从而提升后续计算（如表达量与活性的相关性分析、热图绘制）的统计稳健性。
- 调用位置：在 4.4 节（计算全局 TF 相关性）以及 5.2 节（绘制差异 Peak 热图）中处理矩阵时调用。
- 功能：将高度稀疏的单细胞矩阵按指定的细胞分组和合并大小（
normalize_tf_names&extract_motif_tfs功能：
normalize_tf_names负责清洗 TF 名称（去除括号及冗余字符），并根据物种（人/鼠）转换大小写；extract_motif_tfs进一步用于拆解由::连接的复合转录因子名称。- 目的：解决 Motif 数据库命名与单细胞 RNA 表达矩阵中标准基因命名不一致的问题，确保后续能将表观活性与基因表达准确匹配。
- 调用位置：在 4.4 节的全局相关性计算准备阶段，以及 5.2 节中执行
FindMotifs过滤表达的 TF 时调用。
compute_global_tf_correlation- 功能：整合 RNA 表达矩阵与 ChromVAR 活性矩阵，计算每个转录因子在全局所有细胞群中的表达量与染色质开放打分之间的 Pearson 相关系数及 P 值。
- 目的：提供一个全局视角的参考矩阵，用于评估一个 TF 是否具备“表达水平越高，其调控的表观区域越开放”的正向驱动特征。
- 调用位置：在 4.4 节被独立调用生成全局背景表格，其结果随后在 5.2 节中被用于注释和过滤差异 Peak 富集出的 Motif。
get_group_motif_text&format_motif_grid- 功能：
get_group_motif_text读取 Motif 富集结果表，按倍数和显著性提取 Top 列表；format_motif_grid则将这些字符串格式化为整齐的多行多列网格文本。 - 目的：专为复杂的可视化图表设计，将冗长的 Motif 列表转换为排版规整的注释文本，避免标签在图表中相互重叠。
- 调用位置：在 5.2 节绘制包含右侧 Motif 连线注释的 Heatmap 热图时调用。
- 功能：
sort_clusters- 功能：智能识别包含数字的字符串向量，并按数字的自然大小进行排序。
- 目的：确保细胞群的排序符合人类直觉（例如
Cluster 2排在Cluster 10前面，而不是字母表顺序的后面），从而保证输出表格和图表的坐标轴顺序合理。 - 调用位置：在 5.1 节提取全部分组以及贯穿全文的
for (i in cluster)循环前调用。
R
# 把单细胞矩阵按分组做“伪批量（pseudobulk）聚合”
make_pseudobulk <- function(mat, groups, pb_size = 10) {
mat <- as.matrix(mat)
groups <- as.character(groups)
# 基础检查
stopifnot(ncol(mat) == length(groups))
if (!is.numeric(pb_size) || length(pb_size) != 1 || is.na(pb_size)) {
stop("pb_size 必须是单个数值且不能为 NA")
}
pb_size <- as.integer(pb_size)
if (pb_size < 1) pb_size <- 1L
pb_cols <- list()
pb_group <- character(0)
pb_names <- character(0)
grp_levels <- unique(groups)
for (g in grp_levels) {
idx <- which(groups == g)
if (length(idx) == 0) next
# 用整除分桶，避免 cut(..., breaks<=1) 报错
# 例如 pb_size=10: 1-10 -> pb1, 11-20 -> pb2 ...
split_id <- ((seq_along(idx) - 1L) %/% pb_size) + 1L
chunk_list <- split(idx, split_id)
for (j in seq_along(chunk_list)) {
cidx <- chunk_list[[j]]
if (length(cidx) == 1) {
# 保证返回向量
pb_vec <- as.numeric(mat[, cidx, drop = TRUE])
} else {
pb_vec <- rowMeans(mat[, cidx, drop = FALSE])
}
pb_cols[[length(pb_cols) + 1L]] <- pb_vec
pb_group <- c(pb_group, g)
pb_names <- c(pb_names, paste0(g, "_pb", j))
}
}
# 全空保护（理论上少见，但更稳）
if (length(pb_cols) == 0) {
pb_mat <- matrix(numeric(0), nrow = nrow(mat), ncol = 0)
rownames(pb_mat) <- rownames(mat)
return(list(mat = pb_mat, group = character(0), pb_names = character(0)))
}
pb_mat <- do.call(cbind, pb_cols)
colnames(pb_mat) <- make.unique(pb_names)
rownames(pb_mat) <- rownames(mat)
list(mat = pb_mat, group = pb_group, pb_names = colnames(pb_mat))
}R
normalize_tf_names <- function(tf_vec, species_name) {
tf_vec <- as.character(tf_vec)
tf_vec <- trimws(gsub("\\s*\\(.*\\)", "", tf_vec))
tf_vec <- tf_vec[!is.na(tf_vec) & tf_vec != ""]
if (length(tf_vec) == 0) return(character(0))
sp <- tolower(species_name)
if (sp %in% c("mouse", "rat")) {
tf_vec <- paste0(
toupper(substr(tolower(tf_vec), 1, 1)),
substr(tolower(tf_vec), 2, nchar(tolower(tf_vec)))
)
} else if (sp == "zebrafish") {
tf_vec <- tolower(tf_vec)
} else {
tf_vec <- toupper(tf_vec)
}
unique(tf_vec)
}R
#处理TF名称的函数
extract_motif_tfs <- function(motif_name, species_name) {
motif_name <- as.character(motif_name)
parts <- unlist(strsplit(motif_name, "::", fixed = TRUE))
normalize_tf_names(parts, species_name = species_name)
}R
#TF活性和TF表达相关性（全局计算，使用拟bulk）
compute_global_tf_correlation <- function(obj, expr_assay, spe, output_dir, group_vec, pb_size = 10) {
if (is.na(expr_assay) || !("chromvar" %in% Assays(obj))) return(NULL)
# 提取全集矩阵
expr_mat <- GetAssayData(obj, assay = expr_assay, slot = "data")
cv_all <- GetAssayData(obj, assay = "chromvar", slot = "data")
# 过滤在少于 2% 的细胞中表达的基因
expr_pct <- rowMeans(expr_mat > 0)
expressed_tfs <- normalize_tf_names(rownames(expr_mat)[expr_pct >= 0.01], species_name = spe)
# 伪批量聚合
pb_expr <- make_pseudobulk(expr_mat, group_vec, pb_size = pseudobulk_size)$mat
pb_cv <- make_pseudobulk(cv_all, group_vec, pb_size = pseudobulk_size)$mat
# 获取所有的 motif 信息
motif_obj <- Motifs(obj[["ATAC"]])
all_motif_ids <- names(motif_obj@motif.names)
all_motif_names <- as.character(motif_obj@motif.names)
out_list <- vector("list", length(all_motif_ids))
for (idx in seq_along(all_motif_ids)) {
motif_id <- all_motif_ids[idx]
motif_name <- all_motif_names[idx]
if (!(motif_id %in% rownames(pb_cv))) next
# 提取并清洗该 Motif 对应的 TF 名称
tf_vec <- extract_motif_tfs(motif_name, species_name = spe)
tf_vec <- unique(tf_vec[tf_vec %in% expressed_tfs])
if (length(tf_vec) == 0) next
act_vec <- as.numeric(pb_cv[motif_id, ])
one_motif <- lapply(tf_vec, function(tf_name) {
if (!(tf_name %in% rownames(pb_expr))) return(NULL)
expr_vec <- as.numeric(pb_expr[tf_name, ])
if (sd(expr_vec) == 0 || sd(act_vec) == 0) {
r_val <- NA_real_
p_val <- NA_real_
} else {
ct <- suppressWarnings(cor.test(expr_vec, act_vec, method = "pearson"))
r_val <- unname(ct$estimate)
p_val <- ct$p.value
}
data.frame(
motif_id = motif_id,
motif_name = motif_name,
TF = tf_name,
expr_pct = unname(expr_pct[tf_name]),
expr_mean = mean(expr_vec),
activity_mean = mean(act_vec),
pearson_cor = r_val,
p_value = p_val,
stringsAsFactors = FALSE
)
})
out_list[[idx]] <- bind_rows(one_motif)
}
tf_corr <- bind_rows(out_list)
if (is.null(tf_corr) || nrow(tf_corr) == 0) return(NULL)
# 按照相关性排序并保存全局表格
tf_corr <- tf_corr[order(tf_corr$pearson_cor, decreasing = TRUE), , drop = FALSE]
dir.create(paste0(output_dir, '/DIFF/02_motif'), recursive = TRUE, showWarnings = FALSE)
write.table(tf_corr, paste0(output_dir, "/DIFF/02_motif/Global_TF_ActExpr_Correlation.xls"),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = "\t")
return(tf_corr)
}R
#热图TF标签排版
format_motif_grid <- function(motif_names, n_col = 1, n_row = 20) {
motif_names <- motif_names[!is.na(motif_names) & motif_names != ""]
if (length(motif_names) == 0) return("")
max_n <- n_col * n_row
motif_names <- motif_names[seq_len(min(length(motif_names), max_n))]
if (length(motif_names) < max_n) {
motif_names <- c(motif_names, rep("", max_n - length(motif_names)))
}
mat <- matrix(motif_names, ncol = n_col, byrow = TRUE)
line_txt <- apply(mat, 1, function(x) {
paste(x, collapse = " ")
})
line_txt <- gsub("\\s+$", "", line_txt)
paste(line_txt, collapse = "\n")
}
expr_assay <- "RNA"R
#提取motif富集文件， 按照fold.enrichement排序，去括号，去空格，去重
get_group_motif_text <- function(motif_file, fe_cut = 1.5, p_adj_cut = 0.05, top_n = 20) {
if (!file.exists(motif_file)) return("")
motifs_df <- tryCatch({
read.table(motif_file, sep = "\t", header = TRUE, stringsAsFactors = FALSE)
}, error = function(e) NULL)
if (is.null(motifs_df) || nrow(motifs_df) == 0) return("")
if (!all(c("motif.name", "fold.enrichment", "p.adjust") %in% colnames(motifs_df))) return("")
motifs_df <- motifs_df[
!is.na(motifs_df$fold.enrichment) &
!is.na(motifs_df$p.adjust) &
motifs_df$fold.enrichment > fe_cut &
motifs_df$p.adjust < p_adj_cut,
, drop = FALSE
]
if (nrow(motifs_df) == 0) return("")
motifs_df <- motifs_df[order(motifs_df$fold.enrichment, decreasing = TRUE), , drop = FALSE]
top_motifs <- head(motifs_df$motif.name, top_n)
top_motifs <- gsub("\\s*\\(.*\\)", "", top_motifs)
top_motifs <- trimws(top_motifs)
top_motifs <- top_motifs[top_motifs != ""]
top_motifs <- unique(top_motifs)
top_motifs <- head(top_motifs, top_n)
format_motif_grid(top_motifs, n_col = 1, n_row = 20)
}3. 数据加载与预处理
在进行差异 Peak 分析之前，需要先加载 Seurat 对象及元数据，并根据分析需求对目标细胞群体进行提取、过滤和下采样，以确保后续比较的准确性和统计学有效性。
3.1 加载 Seurat 对象与元数据
首先，读取输入的 .rds 对象。如果提供了独立的元数据文件（meta.tsv），程序会将其读取并合并到 Seurat 对象中，为后续的分组比较（如 case_group vs control_group）和细胞类型鉴定（clusters_col）提供依据。
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)
}
obj <- subset(obj, subset = !!sym(clusters_col) %in% celltypes)💡 Note
注意事项:meta.tsv中的barcode必须与input.rds中 Seurat 对象的细胞名（colnames）完全一致。
3.2 细胞类型筛选与下采样平衡
为避免某些数量极大的细胞类型在分析中产生统计学偏倚，程序提供了细胞类型的提取和下采样功能：
- 筛选目标细胞：根据
celltypes参数，仅保留感兴趣的细胞类型（如B cells, T cells等）。 - 执行下采样：当启用
downsample = TRUE时，程序会对每个细胞类型进行随机抽样，将细胞数量限制在downsample_num（如 3000 个）以内，从而平衡不同细胞类型的细胞数。
R
if (downsample == 'TRUE'){
print(levels(Idents(obj)))
Idents(obj) <- clusters_col
obj <-subset(obj, downsample = downsample_num)
print('downsample')
print(dim(obj@meta.data))
}3.3 修改 fragment 文件路径
在单细胞对象转移或跨服务器读取时，ATAC 模态下的片段文件路径经常会失效。程序会遍历 Seurat 对象中 ATAC Assay 的所有 Fragment 记录，提取文件名，并将其路径重定向到当前的输出目录（outdir），确保后续计算 TSS 富集、CoveragePlot 画图等需要底层序列的操作能够正常寻址。
R
n_fragments <- length(obj@assays$ATAC@fragments)
for (i in 1:n_fragments) {
original_path <- obj@assays$ATAC@fragments[[i]]@path
filename <- basename(original_path)
obj@assays$ATAC@fragments[[i]]@path <- paste0(fpath, '/', filename)
}4. Motif 活性分析
4.1 参考基因组处理
为后续的 TF（转录因子）分析做准备，需要先处理基因组序列文件，读取 genome.fa 并将其处理为 DNAStringSet 对象，同时对齐对象中 ATAC peaks 的染色体名称（保留共有的染色体信息）。
R
DefaultAssay(obj) <- "ATAC"
Idents(obj) <- clusters_col
library(Biostrings)
genome <- readDNAStringSet(ref_genome)
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.2 Motif 数据获取
在执行转录因子 Motif 富集分析之前，我们需要先获取参考的 Motif 位置频率矩阵（PFM，Position Frequency Matrix）。然后使用 AddMotifs() 将 Motif 信息添加至 ATAC Assay
R
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(paste0("目前选项只有 'human' 或 'mouse',","不支持该物种: '", species,
"'\n请自行安装对应物种的基因注释数据库。"))
}
#将 Motif 信息添加至 ATAC Assay
obj <- AddMotifs(object = obj, genome = genome, pfm = pfm)💡 Motif 数据库选择指南
在进行 Motif 分析时，选择合适的参考数据库至关重要，建议根据您的研究物种遵循以下原则：
- 人、小鼠、大鼠、鸡等常见模式动物：JASPAR 数据库中注释极其完善，直接在代码中使用
JASPAR2020提取指定物种的矩阵即可。- 其他脊椎动物（如某些家畜或野生动物）：如果在 JASPAR 中指定该物种时发现可用的 Motif 数量极少（只有几个），为了保证分析的覆盖度，建议直接获取整个“脊椎动物”大类的 Motif 集合。修改代码如下：
pfm <- getMatrixSet(x = JASPAR2020, opts = list(tax_group = "vertebrates", all_versions = TRUE))- 植物物种（如水稻、拟南芥等）：由于 JASPAR 对植物的收录有限，强烈建议前往 PlantTFDB 数据库下载该物种专属的
.meme文件，并参考下方的代码框将其转换为PFMatrixList后导入使用。R# 如果还未安装相关包，请先取消下方注释进行安装： # BiocManager::install("universalmotif") library(universalmotif) # 1. 读取从 PlantTFDB 下载的 meme 文件 meme_path <- "/您的路径/Osi_TF_binding_motifs.meme" raw_motifs <- read_meme(meme_path) # 2. 将 motifs 转换为 TFBSTools 的 PFMatrix 对象 pfm <- convert_motifs(raw_motifs, class = "TFBSTools-PFMatrix") # 3. 确保最终对象是 PFMatrixList 类型 (Signac 严格要求的格式) if (is.list(pfm)) { pfm <- do.call(PFMatrixList, pfm) } else if (is(pfm, "PFMatrix")) { pfm <- PFMatrixList(pfm) }
4.3 Motif活性分析
利用 RunChromVAR 函数，计算每个单细胞在全基因组范围内针对各个转录因子 Motif 的开放程度评分，从而评估不同转录因子的潜在全局活性变化活性。
R
obj <- RunChromVAR(object = obj, genome = genome)4.4 全局 TF 活性与表达相关性分析
一个真实发挥作用的“核心 TF ”，不仅需要其结合靶点的染色质处于异常开放状态，其本身的基因也必须被积极转录翻译成蛋白。因此，为精准锁定核心 TF ，本步骤执行了以下两项严格的筛选与评估：
- 条件一 真实转录支撑：提取所有 Motif 对应的潜在 TF，并过滤掉在单细胞 RNA 数据中不表达，即检测率低于 1%的假阳性因子。
- 条件二 强烈的协同联动：在全细胞层面上，进行拟bulk处理，计算每个 TF 的 RNA 表达量与其对应 Motif 结合位点活性的 Pearson 相关系数。
在下方的可视化结果中，我们通过 全局 TF 活性与表达相关性散点图展示了所有参评 TF 的相关系数排序全貌。图表两端高亮标注了排名前 15 的最强正相关和最强负相关的 TF 。
R
# 立即执行全局 TF 相关性计算并保存
global_group_vec <- if (group_comparison != "") {
paste0(obj@meta.data[[clusters_col]], "_", obj@meta.data[[group_comparison]])
} else {
as.character(obj@meta.data[[clusters_col]])
}
global_tf_corr <- compute_global_tf_correlation(obj, expr_assay = "RNA", spe = species, output_dir = output, group_vec = global_group_vec, pb_size = pseudobulk_size)R
motif_dir <- file.path(output, "DIFF", "02_motif")
if (is.na(expr_assay) || !(expr_assay %in% names(obj@assays))) expr_assay <- "RNA"
cv_all <- GetAssayData(obj, assay = "chromvar", slot = "data")
expr_all <- GetAssayData(obj, assay = expr_assay, slot = "data")
# 1. 提取所有Motif和其对应的TF name
mm <- tryCatch(obj@assays$ATAC@motifs@motif.names, error = function(e) NULL)
if (!is.null(mm)) {
all_motifs <- data.frame(
motif_id = names(mm),
motif_name = as.character(mm),
stringsAsFactors = FALSE
)
# 筛选出表达的TF
expr_pct <- rowMeans(expr_all > 0)
expressed_tfs <- normalize_tf_names(rownames(expr_all)[expr_pct >= 0.01], species_name = species)
# 将TF name规范为表达矩阵里面有的基因名称
motif_tf_list <- lapply(all_motifs$motif_name, extract_motif_tfs, species_name = species)
all_motifs$tf_in_rna <- vapply(motif_tf_list, function(x) any(x %in% expressed_tfs), logical(1))
all_motifs <- all_motifs[all_motifs$tf_in_rna, , drop = FALSE]
# 2. 全细胞范围内容计算TF的活性和表达的相关性
# 保留拟bulk处理的相关参数
group_vec <- as.character(obj@meta.data[[clusters_col]])
pb_cv <- make_pseudobulk(cv_all, group_vec, pb_size = pseudobulk_size)
pb_expr <- make_pseudobulk(expr_all, group_vec, pb_size = pseudobulk_size)
common_cols <- intersect(colnames(pb_cv$mat), colnames(pb_expr$mat))
pb_cv_mat <- pb_cv$mat[, common_cols, drop = FALSE]
pb_expr_mat <- pb_expr$mat[, common_cols, drop = FALSE]
out_list <- list()
for (i in seq_len(nrow(all_motifs))) {
motif_id <- all_motifs$motif_id[i]
motif_name <- all_motifs$motif_name[i]
if (!(motif_id %in% rownames(pb_cv_mat))) next
tf_vec <- extract_motif_tfs(motif_name, species_name = species)
tf_vec <- unique(tf_vec[tf_vec %in% expressed_tfs])
if (length(tf_vec) == 0) next
act_vec <- as.numeric(pb_cv_mat[motif_id, ])
for (tf_name in tf_vec) {
if (!(tf_name %in% rownames(pb_expr_mat))) next
expr_vec <- as.numeric(pb_expr_mat[tf_name, ])
if (sd(expr_vec) == 0 || sd(act_vec) == 0) {
r_val <- NA_real_
p_val <- NA_real_
} else {
ct <- suppressWarnings(cor.test(expr_vec, act_vec, method = "pearson"))
r_val <- unname(ct$estimate)
p_val <- ct$p.value
}
out_list[[length(out_list) + 1]] <- data.frame(
motif_id = motif_id,
motif_name = motif_name,
TF = tf_name,
pearson_cor = r_val,
p_value = p_val,
stringsAsFactors = FALSE
)
}
}
if (length(out_list) > 0) {
global_tf_corr <- do.call(rbind, out_list)
global_tf_corr <- global_tf_corr[!is.na(global_tf_corr$pearson_cor), , drop = FALSE]
global_tf_corr <- global_tf_corr[order(global_tf_corr$pearson_cor, decreasing = TRUE), , drop = FALSE]
# 3. 将相关性结果表格保存在02_motif目录下
global_tf_corr$rank <- seq_len(nrow(global_tf_corr))
write.table(global_tf_corr, file.path(motif_dir, "All_TF_Global_Correlation.xls"), sep = "\t", quote = FALSE, row.names = FALSE)
message("完成：全局TF活性和表达相关性计算，结果输出到 ", file.path(motif_dir, "All_TF_Global_Correlation.xls"))
# 4. 画一个整体的相关性散点图
library(ggrepel)
# 为了打标签，我们先对每个唯一的 TF 找出其相关性的绝对最大值所在的行
# （防止同一个基因因为对应多个相似的 Motif ID 而在图中被标注多次）
unique_tf_corr <- global_tf_corr
unique_tf_corr <- unique_tf_corr[order(abs(unique_tf_corr$pearson_cor), decreasing = TRUE), ]
unique_tf_corr <- unique_tf_corr[!duplicated(unique_tf_corr$TF), ]
# 恢复按相关系数排序
unique_tf_corr <- unique_tf_corr[order(unique_tf_corr$pearson_cor, decreasing = TRUE), ]
# 挑选最正相关的唯一 TF Top 15 和 最负相关的唯一 TF Top 15
top_pos <- head(unique_tf_corr, 15)
top_neg <- tail(unique_tf_corr, 15)
# 构造用于打标签的数据框
label_data <- rbind(top_pos, top_neg)
label_data$label_text <- label_data$TF
p_scatter <- ggplot(global_tf_corr, aes(x = rank, y = pearson_cor)) +
geom_point(color = "#d7191c", alpha = 1, size = 0.5) +
geom_text_repel(
data = label_data,
aes(label = label_text),
size = 2,
box.padding = 0.4,
point.padding = 0.4,
force = 0.1,
max.iter = 20000,
max.time = 2,
max.overlaps = Inf,
min.segment.length = 0,
color = "black",
# fontface = "bold",
segment.color = "grey50",
segment.size = 0.2,
seed = 1
) +
theme_classic() +
theme(
legend.position = "none",
plot.margin = margin(5.5, 40, 5.5, 5.5),
axis.text.x = element_text(size = 8, color = "black"), # 横坐标刻度标签
axis.title.x = element_text(size = 8),
axis.text.y = element_text(size = 8, color = "black"), # 纵坐标刻度标签
axis.title.y = element_text(size = 8)
) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.15))) +
coord_cartesian(clip = "off") +
labs(x = "TF Rank (by Correlation)", y = "Pearson Correlation")
ggplot2::ggsave(file.path(motif_dir, "All_TF_Global_Correlation_Scatter.pdf"), p_scatter, width = 6, height = 4)
ggplot2::ggsave(file.path(motif_dir, "All_TF_Global_Correlation_Scatter.png"), p_scatter, width = 6, height = 4, dpi = 300, bg = "white")
} else {
message("未能计算出任何TF的全局相关性。")
}
}R
data_dir <- paste0(output,"/DIFF/02_motif")
global_corr_file <- file.path(data_dir, "All_TF_Global_Correlation.xls")
if (file.exists(global_corr_file)) {
global_corr <- read.delim(global_corr_file, sep = "\t", stringsAsFactors = FALSE)
if (nrow(global_corr) > 0) {
head(global_corr)
}
}| motif_id | motif_name | TF | pearson_cor | p_value | rank | |
|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <dbl> | <dbl> | <int> | |
| 1 | MA0154.2 | EBF1 | EBF1 | 0.9153930 | 0 | 1 |
| 2 | MA0014.2 | PAX5 | PAX5 | 0.8893159 | 0 | 2 |
| 3 | MA0140.2 | GATA1::TAL1 | GATA1 | 0.8841162 | 0 | 3 |
| 4 | MA0035.4 | GATA1 | GATA1 | 0.8747053 | 0 | 4 |
| 5 | MA0140.2 | GATA1::TAL1 | TAL1 | 0.8287265 | 0 | 5 |
| 6 | MA0091.1 | TAL1::TCF3 | TCF3 | 0.8281944 | 0 | 6 |
💡 解读指南
该表详细记录了在全局所有细胞范围内，各 TF 的 RNA 表达量与其对应 Motif 结合位点活性之间的定量关联分析结果。表格已默认按相关系数降序排列，排在前面的即为最有可能驱动染色质开放的核心激活因子。
- motif_id / motif_name：来自 JASPAR 等数据库的 Motif 编号及对应的可读名称。
- TF：与该 Motif 序列相匹配，且在当前单细胞数据中检测到实际表达的具体 TF 基因名称，如存在多个同家族基因，表格中会分行独立展示。
- pearson_cor：核心指标，皮尔逊相关系数。正值表示表达量升高伴随活性增强，负值表示表达量升高伴随活性减弱。绝对值越大相关性越强。
- p_value：相关性检验的显著性 P 值。
- rank：基于
pearson_cor从大到小排列的全局名次。排名越靠前，即 Rank 越小，代表该 TF 作为正向驱动因子的潜力越大。
R
options(repr.plot.width = 10, repr.plot.height = 5)
p_scatter
💡 解读结论
散点图两端（左上角和右下角）的这些核心 TF，是结合了表观开放与基因表达双重特性的主导调控枢纽。特别是处于强正相关区域的 TF，具备极强的证据表明其基因表达的上升切实主导了下游靶点染色质图谱的开放，建议将其作为后续深入分析的首选靶标。
4.5 细胞群间/组间 TF 活性差异分析
在完成了全基因组范围的 Motif 活性计算后，我们同样利用 presto 包的 wilcoxauc 函数对 chromvar Assay 的打分矩阵进行快速的 Wilcoxon 秩和检验，以鉴定出在特定细胞群或特定状态下异常活跃的转录因子（TF）。
与基因表达差异分析的逻辑一致，这里的分析同样分为两种模式：
- 单聚类 TF Marker 鉴定模式（Cluster vs All）
- 触发条件：当
group_comparison参数为空（""）时。 - 分析逻辑：程序计算每个细胞群相较于其他所有细胞的 TF 活性差异。过滤出显著上调/下调的转录因子（
p_val_adj < 0.05且绝对logFC > logFC阈值）后，保存为总表All_Clusters_TF_Markers.xls。 - 可视化：为了直观展示细胞群的特异性调控网络，程序会提取每个细胞群排名前 10 的核心 TF，将单细胞数据进行伪批量化（Pseudobulk）平滑处理后，绘制一张包含所有细胞群的全局 TF 活性 Z-score 热图（
All_Clusters_TF_Activity_Heatmap.pdf）。
- 触发条件：当
- 细胞类型内组间 TF 差异模式（Case vs Control）
- 触发条件：当
group_comparison参数不为空时。 - 分析逻辑：程序会遍历每一种细胞类型（如 B cells），在同一种细胞内部对比实验组（
case_group）与对照组（control_group）的 TF 活性差异。过滤出的显著结果将按细胞群单独保存（如c_Bcells_TF_Markers.xls）。 - 可视化：针对这一模式，程序会自动为每一个细胞群生成对应的火山图（Volcano Plot）。火山图中会用不同颜色高亮显著变化（Up/Down）的转录因子，并智能地在图上标注出
logFC变化最剧烈的前 15 个上调与下调 TF 的真实名称（如c_Bcells_Case_vs_Control_TF_Activity_Volcano.pdf），帮助研究者快速锁定驱动疾病或状态改变的关键调控因子。
- 触发条件：当
R
if (group_comparison == "") {
tf_markers <- wilcoxauc(
obj,
group_by = clusters_col,
assay = "data",
seurat_assay = "chromvar"
)
colnames(tf_markers)[2] <- "cluster"
tf_markers$gene <- tf_markers$feature
tf_markers$p_val_adj <- tf_markers$padj
tf_markers$avg_log2FC <- tf_markers$logFC
if (nrow(tf_markers) > 0) {
# 1. 过滤显著差异数据（根据 logFC 绝对值 和 padj）
sig_tf <- tf_markers[tf_markers$p_val_adj < 0.05 & abs(tf_markers$logFC) > logFC, , drop = FALSE]
# 2. 覆盖保存过滤后的结果
write.table(sig_tf, paste0(output, "/DIFF/02_motif/All_Clusters_TF_Markers.xls"),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = "\t")
# 3. 恢复使用热图绘制 (细胞群间)
top_tfs_df <- sig_tf %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
top_tfs <- unique(top_tfs_df$gene)
if (length(top_tfs) > 0) {
cell_order <- order(factor(obj@meta.data[[clusters_col]], levels = cluster))
cells_use <- rownames(obj@meta.data)[cell_order]
mat <- GetAssayData(obj, assay = "chromvar", slot = "data")[top_tfs, cells_use, drop = FALSE]
mat <- as.matrix(mat)
cluster_vec <- as.character(obj@meta.data[cells_use, clusters_col, drop = TRUE])
pb <- make_pseudobulk(mat, cluster_vec, pb_size = pseudobulk_size)
mat_pb <- pb$mat
mat_scaled <- t(scale(t(mat_pb)))
mat_scaled[is.na(mat_scaled)] <- 0
mat_scaled[mat_scaled > 2] <- 2
mat_scaled[mat_scaled < -2] <- -2
motif_names <- tryCatch({
motif_obj <- obj@assays$ATAC@motifs
motif_info <- motif_obj@motif.names
sapply(rownames(mat_scaled), function(x) {
if (x %in% names(motif_info)) return(motif_info[[x]])
x
})
}, error = function(e) rownames(mat_scaled))
rownames(mat_scaled) <- make.unique(as.character(motif_names))
col_df <- data.frame(Cluster = pb$group, row.names = colnames(mat_scaled), stringsAsFactors = FALSE)
col_anno <- HeatmapAnnotation(df = col_df, show_annotation_name = FALSE)
ht_tf <- Heatmap(
mat_scaled,
cluster_rows = TRUE,
cluster_columns = FALSE,
show_row_names = TRUE,
row_names_gp = gpar(fontsize = 8),
show_column_names = FALSE,
show_row_dend = FALSE,
top_annotation = col_anno,
name = "TF_Activity_Z",
use_raster = TRUE
)
pdf(paste0(output, "/DIFF/02_motif/All_Clusters_TF_Activity_Heatmap.pdf"),
width = 12, height = max(6, nrow(mat_scaled) * 0.2))
draw(ht_tf)
dev.off()
}
}
} else {
obj_group <- obj
Idents(obj_group) <- group_comparison
for (i in unique(obj_group@meta.data[[clusters_col]])) {
cells_in_cluster <- rownames(obj_group@meta.data)[obj_group@meta.data[[clusters_col]] == i]
if (length(cells_in_cluster) < 10) next
sub_obj <- subset(obj_group, cells = cells_in_cluster)
valid_groups <- unique(sub_obj@meta.data[[group_comparison]])
if (!(case_group %in% valid_groups) || !(control_group %in% valid_groups)) next
tf_markers_all <- wilcoxauc(
sub_obj,
group_by = group_comparison,
assay = "data",
seurat_assay = "chromvar",
groups_use = c(case_group, control_group)
)
tf_markers_all$gene <- tf_markers_all$feature
tf_markers_all$p_val_adj <- tf_markers_all$padj
tf_markers_all$avg_log2FC <- tf_markers_all$logFC
if (nrow(tf_markers_all) > 0) {
# 1. 过滤显著差异数据 (保留两个组的数据)
sig_tf <- tf_markers_all[tf_markers_all$p_val_adj < 0.05 & abs(tf_markers_all$logFC) > logFC, , drop = FALSE]
# 2. 保存过滤后的数据
write.table(sig_tf, paste0(output, "/DIFF/02_motif/c_", gsub(" ", "", i), "_TF_Markers.xls"),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = "\t")
# 3. 画火山图 (火山图通常只看 case vs control 的 logFC，所以提取 case 组的数据来画)
volcano_data <- tf_markers_all[tf_markers_all$group == case_group, , drop = FALSE]
volcano_data$Status <- "Not Significant"
volcano_data$Status[volcano_data$padj < 0.05 & volcano_data$logFC > logFC] <- "Up"
volcano_data$Status[volcano_data$padj < 0.05 & volcano_data$logFC < -logFC] <- "Down"
library(ggrepel)
# 由于 wilcoxauc 已经限制为 case_group 的结果，这里只画一个火山图
vd_g <- volcano_data
up_genes <- vd_g[vd_g$Status == "Up", ]
down_genes <- vd_g[vd_g$Status == "Down", ]
top_up <- up_genes[order(up_genes$logFC, decreasing = TRUE), ]
if (nrow(top_up) > 15) top_up <- top_up[1:15, ]
top_down <- down_genes[order(down_genes$logFC, decreasing = FALSE), ]
if (nrow(top_down) > 15) top_down <- top_down[1:15, ]
label_data <- rbind(top_up, top_down)
# 将 motif_id 转换为实际名字
label_data$label <- tryCatch({
motif_info <- obj@assays$ATAC@motifs@motif.names
sapply(label_data$feature, function(x) if (x %in% names(motif_info)) motif_info[[x]] else x)
}, error = function(e) label_data$feature)
title_txt <- paste0(i, "(", case_group, " vs ", control_group, ")")
out_prefix <- paste0(output, "/DIFF/02_motif/c_", gsub(" ", "", i), "_", gsub(" ", "", case_group), "_vs_", gsub(" ", "", control_group), "_TF_Activity_Volcano")
p_volcano <- ggplot(vd_g, aes(x = logFC, y = -log10(padj), color = Status)) +
geom_point(alpha = 0.8, size = 1.5) +
scale_color_manual(values = c("Up" = "#d73027", "Down" = "#4575b4", "Not Significant" = "#e0e0e0")) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black", alpha = 0.5) +
geom_vline(xintercept = c(-logFC, logFC), linetype = "dashed", color = "black", alpha = 0.5) +
geom_text_repel(data = label_data, aes(label = label), size = 3, color = "black",
max.overlaps = 50, segment.color = "grey50", segment.size = 0.5) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right",
legend.title = element_blank(),
plot.title = element_text(hjust = 0.5, face = "bold", size = 14)
) +
labs(title = title_txt, x = "log2(Fold Change)", y = "-log10(padj)")
ggplot2::ggsave(paste0(out_prefix, ".pdf"), plot = p_volcano, width = 7, height = 6)
ggplot2::ggsave(paste0(out_prefix, ".png"), plot = p_volcano, width = 7, height = 6, dpi = 300)
}
}
}R
options(repr.plot.width = 10, repr.plot.height = 7)
if (group_comparison == "") {
draw(ht_tf)
}else{
p_volcano
}
💡 生物学意义
该步骤的核心目的在于从海量的“开放信号变化”中，提炼出具有生物学解释力的调控因子线索。活性显著上调的 TF 通常暗示其在目标细胞状态中正积极参与转录激活或维持细胞表型；而活性下调则提示其调控网络的沉默。最终输出的差异 TF 列表，是后续构建核心调控网络与设计机制验证实验的最佳候选集合。
R
data_dir <- paste0(output,"/DIFF/02_motif")
# 1. 收集并合并 TF 活性差异表格 (TF_Markers.xls)
if (group_comparison == "") {
marker_files <- list.files(path = data_dir, pattern = "TF_Markers\\.xls$", full.names = TRUE)
} else {
# 组间差异模式：排除掉单聚类的 All_Clusters_TF_Markers.xls 等不相关的文件，仅保留各个 cluster 分组比较生成的 marker
marker_files <- list.files(path = data_dir, pattern = "^c_.*_TF_Markers\\.xls$", full.names = TRUE)
}
merged_markers <- data.frame()
for (file in marker_files) {
df <- tryCatch(read.table(file, sep="\t", header=TRUE, stringsAsFactors = FALSE), error = function(e) NULL)
if (!is.null(df) && nrow(df) > 0) {
if ("cluster" %in% colnames(df)) df$cluster <- as.character(df$cluster)
# 组间模式下补充 cluster 名称（如果是按照 c_XXX_TF_Markers 命名的）
if (group_comparison != "" && !("cluster" %in% colnames(df))) {
# 尝试从文件名解析 cluster name: c_Bcells_TF_Markers.xls -> Bcells
base_name <- basename(file)
cl_name <- gsub("^c_", "", gsub("_TF_Markers\\.xls$", "", base_name))
df$cluster <- cl_name
}
# 确保当前已有的 merged_markers 也不是空的，或者这是第一次合并
if (nrow(merged_markers) == 0) {
merged_markers <- df
} else {
merged_markers <- bind_rows(merged_markers, df)
}
}
}
if (group_comparison != "" && nrow(merged_markers) > 0) {
# 如果是组间差异模式，将其保存到 all_TF_Markers.xls 供交互展示和下载
write.table(merged_markers[which(merged_markers$logFC > 0),],
file = paste0(output, '/DIFF/02_motif/all_TF_Markers.xls'),
quote = F, row.names = F, col.names = T, sep = '\t')
}
head(merged_markers)| feature | group | avgExpr | logFC | statistic | auc | pval | padj | pct_in | pct_out | gene | p_val_adj | avg_log2FC | cluster | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | <int> | <chr> | <dbl> | <dbl> | <chr> | |
| 1 | MA0002.1 | 25030508_pbmc_1_arc | -1.40086958 | -0.4710262 | 150256 | 0.3696662 | 2.112365e-09 | 1.645207e-08 | 100 | 100 | MA0002.1 | 1.645207e-08 | -0.4710262 | Bcells |
| 2 | MA0003.1 | 25030508_pbmc_1_arc | 0.33990583 | -0.4448016 | 155355 | 0.3822110 | 6.215919e-08 | 3.596353e-07 | 100 | 100 | MA0003.1 | 3.596353e-07 | -0.4448016 | Bcells |
| 3 | MA0048.1 | 25030508_pbmc_1_arc | 2.92851842 | -0.9484831 | 160852 | 0.3957349 | 1.659225e-06 | 6.822194e-06 | 100 | 100 | MA0048.1 | 6.822194e-06 | -0.9484831 | Bcells |
| 4 | MA0050.1 | 25030508_pbmc_1_arc | 0.01236445 | -1.1000062 | 107759 | 0.2651133 | 3.704878e-27 | 1.364069e-25 | 100 | 100 | MA0050.1 | 1.364069e-25 | -1.1000062 | Bcells |
| 5 | MA0051.1 | 25030508_pbmc_1_arc | -0.73876620 | -0.9485148 | 116735 | 0.2871964 | 1.392626e-22 | 4.338566e-21 | 100 | 100 | MA0051.1 | 4.338566e-21 | -0.9485148 | Bcells |
| 6 | MA0052.1 | 25030508_pbmc_1_arc | 0.82173660 | -0.4136321 | 155320 | 0.3821249 | 6.079983e-08 | 3.568686e-07 | 100 | 100 | MA0052.1 | 3.568686e-07 | -0.4136321 | Bcells |
💡 解读指南
该表展示了基于 chromVAR 计算的 TF 结合位点可及性的差异分析结果。它可以用来鉴定在特定细胞亚群或疾病组中活性显著改变的核心调控因子。 核心关注
padj显著性和logFC差异倍数两个指标。
- feature：Motif ID，如
MA0494.1。- cluster / group：该 TF 活性发生显著变化的细胞群或者对应的分组。
- auc：曲线下面积 (Area Under the Curve)，用于评估该 TF 活性作为区分该组细胞 marker 的分类效能。数值越接近 1 或 0 区分度越强。
- pval / padj：原始 P 值与多重检验校正后的 P 值。通常以
padj < 0.05作为显著改变的标准。- logFC：平均对数折叠变化。正值表示该 TF 在当前目标组中活性增强，开放度升高；负值表示活性减弱。
- pct_in / pct_out：该 TF 活性在目标组细胞和背景组细胞中检测到的比例。
5. Motif 富集分析
Motif 富集分析与前面的 Motif 活性分析完全不是一个概念：
- Motif 活性分析是计算每个单细胞在全基因组所有区域中针对某个 Motif 的整体开放程度打分（基于背景变异）。
- Motif 富集分析则是基于给定的特定序列集合（如鉴定出的差异 Peak 区域），利用超几何检验等统计学方法，去寻找哪些转录因子 Motif 在这些特定序列中出现的频率显著高于背景序列。这通常用于寻找调控特定生物学过程的核心转录因子。
5.1 差异 Peak 数据集准备
在进行 Motif 富集之前，首先需要利用 wilcoxauc 获取包含目标区域的差异 Peak 集合。与前面的分析逻辑一致，这里同样会根据 group_comparison 参数的值自动切换为两种模式：
- 单聚类 Marker 模式（
group_comparison == ""）：寻找每个细胞群相较于其他所有细胞的特征性差异 Peak。 - 细胞类型内组间差异模式（
group_comparison != ""）：遍历每一种细胞群，提取其在实验组（Case）与对照组（Control）之间发生显著改变的差异 Peak，并将所有细胞群的结果合并导出为all_diffPeak.xls以备后续分析使用。
R
if (group_comparison == "") {
da_peaks <- wilcoxauc(obj, clusters_col, assay='data', seurat_assay='ATAC')
colnames(da_peaks)[2] <- 'cluster'
} else {
all_da_peaks <- list()
obj <- subset(obj, subset = !!sym(group_comparison) %in% c(case_group, control_group))
for (cell_type_analyse in unique(obj@meta.data[[clusters_col]])) {
tryCatch({
sub_obj <- subset(obj, subset = !!sym(clusters_col) %in% cell_type_analyse)
da_peaks <- wilcoxauc(sub_obj,
group_by = group_comparison,
assay = 'data',
seurat_assay = 'ATAC',
groups_use = c(case_group, control_group))
da_peaks$cluster <- as.character(cell_type_analyse)
all_da_peaks[[as.character(cell_type_analyse)]] <- da_peaks
}, error = function(e) {
message(sprintf("Error in cluster %s: %s", cell_type_analyse, e$message))
})
}
da_peaks <- do.call(rbind, all_da_peaks)
write.table(da_peaks,
file = paste0(output, '/DIFF/all_diffPeak.xls'),
quote = F, row.names = F, col.names = T, sep = '\t')
}
saveRDS(obj,paste0(output, '/output.rds'))
write.table(da_peaks[da_peaks$padj < pval_adj & da_peaks$logFC > logFC, ],paste0(output,'/DIFF/01_diffPeak/all_diffPeak_significant.xls'),quote=F,row.names=F,col.names=T,sep='\t')5.2 差异 Peaks 的 Motif 富集分析与可视化
在成功鉴定出各细胞群或实验组间的显著差异 Peaks 后，我们进而对这些特异性开放的染色质区域执行 Motif 富集分析。这一步的核心目的在于：从表观遗传现象（开放）溯源到上游的调控机制，找出真正驱动这些区域开放的潜在主导转录因子（TF）。
Motif 富集与多维证据整合： 我们利用
FindMotifs函数，检测特定转录因子的结合基序在差异 Peaks 中是否被显著富集。为了提高结果的置信度，程序会自动将该富集结果与前面计算的全局 TF 表达-活性相关性矩阵进行关联整合。通过这种多维度的交叉验证，我们能够精准筛选出那些不仅在局部差异 Peak 中高度富集，且在全局层面上其 RNA 表达量与染色质开放活性呈高度正相关的“真实”驱动因子（分析结果将统一输出至02_motif目录）。富集特征的直观可视化展示： 为了帮助研究者更直观地解析调控网络，程序会自动提取统计显著性最高的关键特征进行多维度的可视化：
- Motif 热图联构图：以热图形式展示差异 Peak 的开放信号分布，并在右侧通过连线精准标记这些区域内显著富集的核心 Motif。
- 序列徽标图 (Sequence Logo)：直观展示 Top 富集 Motif 的碱基偏好性。
- 单细胞活性投射图 (UMAP)：将这些 Top Motif 的 ChromVAR 活性打分映射到 UMAP 降维空间，展示其在不同单细胞群体中的异质性分布。
R
cluster <- unique(da_peaks$cluster)
# 通用排序函数
sort_clusters <- function(x) {
if(all(!is.na(suppressWarnings(as.numeric(x))))) {
return(x[order(as.numeric(x))])
} else {
nums <- suppressWarnings(as.numeric(gsub("[^0-9]", "", x)))
if(all(!is.na(nums))) {
return(x[order(nums)])
} else {
return(sort(x))
}
}
}
cluster <- sort_clusters(cluster)
for (i in cluster) {
tryCatch({
DefaultAssay(obj) <- "ATAC"
data <- da_peaks[da_peaks$cluster == i, ]
# 使用正则表达式检查 feature 列的格式 染色体-起始-终止
valid_rows <- grepl("^([^-]+)-(\\d+)-(\\d+)$", data$feature)
if (sum(!valid_rows) > 0) {
message(sprintf("Cluster %s: Removed %d malformed feature rows.", i, sum(!valid_rows)))
data <- data[valid_rows, ]
}
# Determine the groups to iterate over based on group_comparison
if (group_comparison == "") {
# Standard single cluster comparison
analysis_groups <- c("cluster_marker")
} else {
# Two group comparison: case and control
analysis_groups <- c(case_group, control_group)
}
for (grp in analysis_groups) {
print(paste0("正在分析：Cluster ", i, " - Group ", grp))
tryCatch({
if (grp == "cluster_marker") {
# For single cluster, we use padj < pval_adj and logFC > logFC
data.sig <- data[data$padj < pval_adj & data$logFC > logFC, ]
prefix <- paste0('c_', gsub(" ", "", i))
} else {
# For group comparison, filter by group and thresholds
data.sig <- data[data$group == grp & data$padj < pval_adj & data$logFC > logFC, ]
prefix <- paste0('c_', gsub(" ", "", i), '_', gsub(" ", "", grp))
}
if (nrow(data.sig) == 0) {
message(sprintf("No significant peaks found for %s", prefix))
next
}
peaks <- data.sig$feature
# --- FindMotifs 分析核心开始 ---
enriched.motifs <- FindMotifs(object = obj, features = peaks)
enriched.motifs$cluster <- as.character(i)
if (grp != "cluster_marker") enriched.motifs$group <- as.character(grp)
if (!is.na(expr_assay) && nrow(enriched.motifs) > 0 && "motif.name" %in% colnames(enriched.motifs)) {
if (grp == "cluster_marker") {
cells_expr <- rownames(obj@meta.data)[obj@meta.data[[clusters_col]] == i]
} else {
cells_expr <- rownames(obj@meta.data)[obj@meta.data[[clusters_col]] == i & obj@meta.data[[group_comparison]] == grp]
}
if (length(cells_expr) > 0) {
expr_mat <- GetAssayData(obj, assay = expr_assay, slot = "data")[, cells_expr, drop = FALSE]
expr_pct <- rowMeans(expr_mat > 0)
# 修正：将原代码的 spe 改为环境中已有的 species 变量
expressed_tfs <- normalize_tf_names(rownames(expr_mat)[expr_pct >= 0.01], species_name = species)
motif_tf_list <- lapply(enriched.motifs$motif.name, extract_motif_tfs, species_name = species)
enriched.motifs$motif_tf <- vapply(motif_tf_list, function(x) paste(x, collapse = "::"), character(1))
enriched.motifs$tf_in_rna <- vapply(motif_tf_list, function(x) any(x %in% expressed_tfs), logical(1))
enriched.motifs <- enriched.motifs[enriched.motifs$tf_in_rna, , drop = FALSE]
} else {
enriched.motifs <- enriched.motifs[0, , drop = FALSE]
}
}
enriched.motifs <- enriched.motifs[which(enriched.motifs$p.adjust < 0.05 & enriched.motifs$fold.enrichment > 1), , drop = FALSE]
# 修正：将 03_motif 改为 02_motif，解决与后续合并代码的目录冲突问题
write.table(enriched.motifs, paste0(output,'/DIFF/02_motif/', prefix, '_diffPeak_motifs.xls'), quote=F, row.names=F, col.names=T, sep='\t')
if (!is.na(expr_assay) && nrow(enriched.motifs) > 0) {
if (!is.null(global_tf_corr)) {
tf_corr_subset <- global_tf_corr[global_tf_corr$motif_id %in% rownames(enriched.motifs), , drop = FALSE]
if (nrow(tf_corr_subset) > 0) {
motif_meta <- enriched.motifs
motif_meta$motif_id <- rownames(motif_meta)
keep_cols <- c("motif_id", "fold.enrichment", "cluster", "group", "p.adjust")
keep_cols <- keep_cols[keep_cols %in% colnames(motif_meta)]
motif_meta <- motif_meta[, keep_cols, drop = FALSE]
tf_corr_subset <- dplyr::left_join(tf_corr_subset, motif_meta, by = "motif_id")
if (!("fold.enrichment" %in% colnames(tf_corr_subset))) tf_corr_subset$fold.enrichment <- NA_real_
if (!("cluster" %in% colnames(tf_corr_subset))) tf_corr_subset$cluster <- as.character(i)
if (!("group" %in% colnames(tf_corr_subset))) tf_corr_subset$group <- if (grp == "cluster_marker") "" else as.character(grp)
if (!("p.adjust" %in% colnames(tf_corr_subset))) tf_corr_subset$p.adjust <- NA_real_
tf_corr_subset$high_cor <- tf_corr_subset$pearson_cor > cor
tf_corr_subset <- tf_corr_subset[, c(
"motif_id", "motif_name", "TF", "fold.enrichment", "p.adjust", "cluster", "group",
"expr_pct", "expr_mean", "activity_mean", "pearson_cor", "p_value", "high_cor"
), drop = FALSE]
# 修正：将 03_motif 改为 02_motif
write.table(tf_corr_subset, paste0(output, "/DIFF/02_motif/", prefix, "_TF_ActExpr_Correlation.xls"), quote = FALSE, row.names = FALSE, col.names = TRUE, sep = "\t")
}
}
}
}, error = function(e) {
message(sprintf("Error in cluster %s, group %s: %s", i, grp, e$message))
})
}
print(paste0("Cluster ", i, " 分析完成"))
}, error = function(e) {
message(sprintf("Error in cluster %s: %s", i, e$message))
})
}R
tryCatch({
ensure_color_map <- function(values, color_map) {
values <- unique(as.character(values))
values <- values[!is.na(values)]
if (!all(values %in% names(color_map))) {
miss <- setdiff(values, names(color_map))
if (length(miss) > 0) {
extra_cols <- setNames(scales::hue_pal()(length(miss)), miss)
color_map <- c(color_map, extra_cols)
}
}
color_map
}
build_mark <- function(group_vec, label_map) {
lv <- unique(as.character(group_vec))
lv <- lv[!is.na(lv)]
at <- unname(vapply(lv, function(g) {
idx <- which(as.character(group_vec) == g)
idx[(length(idx) + 1) %/% 2]
}, integer(1)))
labels <- unname(label_map[lv])
# 去除尾部可能的多余空格
labels <- gsub("\\s+$", "", labels)
# 新增：找出不是空字符串的索引，只保留有文本的标记
valid_idx <- which(labels != "")
# 只返回有效的 at, labels 和 levels
list(at = at[valid_idx], labels = labels[valid_idx], levels = lv[valid_idx])
}
if (group_comparison == "") {
all_top_peaks <- c()
peak_cluster_anno <- c()
motif_text_anno <- c()
for (i in cluster) {
data.sig <- da_peaks[da_peaks$cluster == i & da_peaks$padj < pval_adj & da_peaks$logFC > logFC, ]
if (nrow(data.sig) > 0) {
data.sig <- data.sig[order(data.sig$logFC, decreasing = TRUE), ]
top_n <- min(nrow(data.sig), 500)
sel_peaks <- data.sig$feature[1:top_n]
all_top_peaks <- c(all_top_peaks, sel_peaks)
peak_cluster_anno <- c(peak_cluster_anno, rep(i, top_n))
prefix <- paste0("c_", gsub(" ", "", i))
motif_file <- paste0(output, "/DIFF/02_motif/", prefix, "_diffPeak_motifs.xls")
motif_text <- get_group_motif_text(motif_file, fe_cut = fe_cut, top_n = 6, n_col = 3, n_row = 2)
motif_text_anno <- c(motif_text_anno, motif_text)
}
}
if (length(all_top_peaks) > 0) {
active_clusters <- unique(peak_cluster_anno)
names(motif_text_anno) <- active_clusters
cell_order <- order(factor(obj@meta.data[[clusters_col]], levels = cluster))
cells_use <- rownames(obj@meta.data)[cell_order]
DefaultAssay(obj) <- "ATAC"
mat <- GetAssayData(obj, assay = "ATAC", slot = "data")[all_top_peaks, cells_use, drop = FALSE]
mat <- as.matrix(mat)
cluster_vec <- as.character(obj@meta.data[cells_use, clusters_col, drop = TRUE])
pb <- make_pseudobulk(mat = mat, groups = cluster_vec, pb_size = pseudobulk_size)
mat <- pb$mat
col_df <- data.frame(Cluster = pb$group, row.names = colnames(mat), stringsAsFactors = FALSE)
mat_scaled <- t(scale(t(mat)))
mat_scaled[is.na(mat_scaled)] <- 0
mat_scaled[mat_scaled > 2] <- 2
mat_scaled[mat_scaled < -2] <- -2
cols_use <- ensure_color_map(col_df$Cluster, group_colors)
col_anno <- HeatmapAnnotation(
df = col_df,
col = list(Cluster = cols_use),
show_annotation_name = FALSE
)
mark <- build_mark(peak_cluster_anno, motif_text_anno)
right_anno <- rowAnnotation(
Motifs = anno_mark(
at = mark$at, # 仍然是你算好的目标行中心
labels = mark$labels, # 文本可长，可多行
which = "row",
side = "right",
labels_gp = gpar(fontsize = 8, lineheight = 1.2,fontface = "bold",col = "black"),
lines_gp = gpar(col = "black", lty = 1), # 线样式
link_width = unit(4, "mm"),
# just = c("right","top"),# 横向“引出段”，越大越明显折线感
padding = unit(0.1, "mm"), # 标签间距
extend = unit(c(1, 1), "mm") # 把标签整体往右推，制造更多折线空
)#,
#width = unit(5, "cm")
)
col_fun <- circlize::colorRamp2(
c(-2, -1, 0, 1, 2),
c("#2c7bb6", "#abd9e9", "#f7f7f7", "#fdae61", "#d7191c")
)
ht <- Heatmap(
mat_scaled,
col = col_fun,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_row_names = FALSE,
show_column_names = FALSE,
show_row_dend = FALSE,
row_title_rot = 0,
top_annotation = col_anno,
right_annotation = right_anno,
name = "Accessibility",
use_raster = TRUE
)
pdf(paste0(output, "/DIFF/01_diffPeak/All_Clusters_DiffPeak_Heatmap.pdf"), width = 12, height = 10)
draw(ht,padding = unit(c(10, 6, 6, 10), "mm"))
dev.off()
}
} else {
for (i in cluster) {
data.sig <- da_peaks[da_peaks$cluster == i & da_peaks$padj < pval_adj & da_peaks$logFC > logFC, ]
if (nrow(data.sig) > 0) {
sel_peaks <- c()
peak_group_anno <- c()
motif_text_map <- setNames(rep("", 2), c(case_group, control_group))
for (grp in c(case_group, control_group)) {
data.grp <- data.sig[data.sig$group == grp, ]
if (nrow(data.grp) > 0) {
data.grp <- data.grp[order(data.grp$logFC, decreasing = TRUE), ]
top_n <- min(nrow(data.grp), 500)
grp_peaks <- data.grp$feature[1:top_n]
sel_peaks <- c(sel_peaks, grp_peaks)
peak_group_anno <- c(peak_group_anno, rep(grp, top_n))
prefix <- paste0("c_", gsub(" ", "", i), "_", gsub(" ", "", grp))
motif_file <- paste0(output, "/DIFF/02_motif/", prefix, "_diffPeak_motifs.xls")
motif_text_map[grp] <- get_group_motif_text(motif_file, fe_cut = fe_cut, top_n = 20)
}
}
cells_in_cluster <- rownames(obj@meta.data)[obj@meta.data[[clusters_col]] == i]
sub_meta <- obj@meta.data[cells_in_cluster, , drop = FALSE]
cell_order <- order(factor(sub_meta[[group_comparison]], levels = c(case_group, control_group)))
cells_use <- cells_in_cluster[cell_order]
if (length(cells_use) > 0 && length(sel_peaks) > 0) {
active_groups <- c(case_group, control_group)[c(case_group, control_group) %in% unique(peak_group_anno)]
motif_text_anno <- motif_text_map[active_groups]
DefaultAssay(obj) <- "ATAC"
mat <- GetAssayData(obj, assay = "ATAC", slot = "data")[sel_peaks, cells_use, drop = FALSE]
mat <- as.matrix(mat)
group_vec <- as.character(sub_meta[cells_use, group_comparison, drop = TRUE])
pb <- make_pseudobulk(mat = mat, groups = group_vec, pb_size = pseudobulk_size)
mat <- pb$mat
col_df <- data.frame(Group = pb$group, row.names = colnames(mat), stringsAsFactors = FALSE)
mat_scaled <- t(scale(t(mat)))
mat_scaled[is.na(mat_scaled)] <- 0
mat_scaled[mat_scaled > 2] <- 2
mat_scaled[mat_scaled < -2] <- -2
cols_use <- ensure_color_map(col_df$Group, group_colors)
col_anno <- HeatmapAnnotation(
df = col_df,
col = list(Group = cols_use),
show_annotation_name = FALSE
)
mark <- build_mark(peak_group_anno, motif_text_anno)
right_anno <- rowAnnotation(
Motifs = anno_mark(
at = mark$at, # 仍然是你算好的目标行中心
labels = mark$labels, # 文本可长，可多行
which = "row",
side = "right",
labels_gp = gpar(fontsize = 8, lineheight = 1.2,fontface = "bold",col = "black"),
lines_gp = gpar(col = "black", lty = 2), # 线样式
link_width = unit(4, "mm"),
# just = c("right","top"),# 横向“引出段”，越大越明显折线感
padding = unit(0.5, "mm"), # 标签间距
extend = unit(c(1, 5), "mm") # 把标签整体往右推，制造更多折线空
),
width = unit(5, "cm")
)
col_fun <- circlize::colorRamp2(
c(-2, -1, 0, 1, 2),
c("#2c7bb6", "#abd9e9", "#f7f7f7", "#fdae61", "#d7191c")
)
ht <- Heatmap(
mat_scaled,
col = col_fun,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_row_names = FALSE,
show_column_names = FALSE,
show_row_dend = FALSE,
row_title_rot = 0,
top_annotation = col_anno,
right_annotation = right_anno,
name = "Accessibility",
column_title = paste0("Cluster: ", i),
use_raster = TRUE
)
pdf(paste0(output, "/DIFF/01_diffPeak/c_", gsub(" ", "", i), "_Group_DiffPeak_Heatmap.pdf"), width = 8, height = 7)
draw(ht,padding = unit(c(10, 6, 6, 10), "mm"))
dev.off()
}
}
}
}
}, error = function(e) {
message("Error in drawing heatmaps: ", e$message)
})R
options(repr.plot.width = 15, repr.plot.height = 7)
draw(ht,padding = unit(c(10, 6, 6, 10), "mm"))
💡 说明
这里的示例热图展示了 T cells 这个细胞群组间差异 Peak 的染色质可及性（Accessibility）分布，并在右侧通过连线标记了这些区域内显著富集的核心转录因子 Motif。
- 热图颜色 (Z-score)：颜色从蓝色到红色（或按设定的色板），反映了特定染色质区域在各细胞群或样本组间的相对开放程度。红色越深表示该区域越开放，蓝色越深表示该区域趋于闭合。
- 横轴 (列)：代表被拟批量（Pseudobulk）平滑化处理后的细胞群（或对比组）。顶部的颜色条（Column Annotation）用于区分不同的细胞类型或样本。
- 纵轴 (行)：代表鉴定出的显著差异 Peak（特异性开放区域）。行没有显示具体名称，但按照其归属的细胞群或分组进行了聚类和排序。
- 右侧标签 (Motifs)：提取了在对应差异 Peak 集合中富集倍数（Fold Enrichment）最高的 Top Motif 名称。通过观察这些标签，可以快速直观地将特异性开放区域与其潜在的上游驱动转录因子联系起来。
5.3 Motif 富集结果汇总
在前面的循环分析中，程序已经为每一个细胞群（或对比组）单独生成了差异 Peak 的 Motif 富集结果文件（如 c_Bcells_diffPeak_motifs.xls）。为了方便后续的全局预览、统计与下游的通路分析，我们需要将这些分散的文件合并为一个完整的总表。
这段代码的核心逻辑如下：
- 自动检索文件：通过正则表达式
motifs\\.xls$扫描02_motif结果目录下所有单独生成的富集结果表。 - 安全读取与合并：利用循环依次读取这些文件。在读取过程中，程序内置了安全机制，会自动跳过大小为 0 的空文件或行数为 0 的空数据框（这通常发生在某些细胞群没有鉴定出显著富集 Motif 的情况下）。
- 数据类型统一：在进行
rbind合并前，强制将所有列转换为字符型（character），避免因个别文件中某列（如cluster纯数字名称）类型推断不一致而导致的合并失败。 - 输出汇总结果：最终生成一个包含所有组别、所有细胞群显著富集 Motif 信息的
merged_data总表，便于研究人员进行全局检索与比较。
R
# 最简单的合并版本
data_dir <- paste0(output, "/DIFF/02_motif")
xls_files <- list.files(path = data_dir,
pattern = "motifs\\.xls$",
full.names = TRUE)
merged_data <- data.frame()
for (file in xls_files) {
# 跳过空文件
if (file.size(file) == 0) next
# 读取数据
df <- read.table(file, sep = "\t", header = TRUE,
stringsAsFactors = FALSE)
# 跳过空数据框
if (nrow(df) == 0) next
# 强制所有列为字符型
df[] <- lapply(df, as.character)
# 合并
merged_data <- rbind(merged_data, df)
}
# 处理cluster列
if (nrow(merged_data) > 0 && "cluster" %in% colnames(merged_data)) {
merged_data$cluster <- as.character(merged_data$cluster)
}
cat("合并了", nrow(merged_data), "行数据\n")
head(merged_data)output
合并了 1283 行数据
| motif | observed | background | percent.observed | percent.background | fold.enrichment | pvalue | motif.name | p.adjust | cluster | group | motif_tf | tf_in_rna | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | |
| 1 | MA1548.1 | 13 | 5506 | 41.9354838709677 | 13.765 | 3.04652988528643 | 0.000113293420229509 | PLAGL2 | 0.0229419175964756 | B cells | 25030508_pbmc_1_arc | PLAGL2 | TRUE |
| 2 | MA1648.1 | 454 | 7960 | 36.6720516962843 | 19.9 | 1.84281666815499 | 3.60213197709286e-44 | TCF12(var.2) | 3.24191877938358e-42 | B cells | XYRD_pbmc_2_arc | TCF12 | TRUE |
| 3 | MA1558.1 | 471 | 8486 | 38.0452342487884 | 21.215 | 1.79331766433129 | 5.39243273131294e-43 | SNAI1 | 3.9707913748759e-41 | B cells | XYRD_pbmc_2_arc | SNAI1 | TRUE |
| 4 | MA0830.1 | 438 | 7946 | 35.3796445880452 | 19.865 | 1.78100400644577 | 1.99335037013957e-38 | TCF4 | 1.15329557129504e-36 | B cells | XYRD_pbmc_2_arc | TCF4 | TRUE |
| 5 | MA0107.1 | 297 | 4474 | 23.9903069466882 | 11.185 | 2.14486427775487 | 2.60148879508546e-38 | RELA | 1.31700370251201e-36 | B cells | XYRD_pbmc_2_arc | RELA | TRUE |
| 6 | MA0105.1 | 277 | 4121 | 22.3747980613893 | 10.3025 | 2.17178335951365 | 2.07812263197394e-36 | NFKB1 | 8.41639665949446e-35 | B cells | XYRD_pbmc_2_arc | NFKB1 | TRUE |
💡 说明
该表展示了“差异开放Peak中哪些Motif被富集”，并提供统计显著性、富集倍数与RNA表达支持信息。每一行对应一个Motif候选，核心关注
p.adjust、fold.enrichment和tf_in_rna三个指标。
- motif：Motif ID，如
MA0494.1，通常来自JASPAR，用于唯一标识序列模式。- motif.name：可读的TF/Motif名称，如
Bach1::Mafk、RUNX1，双TF一般表示复合或家族相关模式。- observed：目标差异Peak集合中命中该Motif的Peak数量。
- background：背景Peak集合中命中该Motif的数量。
- percent.observed / percent.background：目标集与背景集中Motif命中百分比，用于直观看富集方向。
- fold.enrichment：富集倍数，即目标比例/背景的比例，>1表示富集，数值越大富集越强。
- pvalue：原始显著性P值。
- p.adjust：多重检验校正后的P值，通常以此作为显著性判断主标准，认为p.adjust < 0.05为显著富集。
- cluster / group：细胞群和比较组别。
- motif_tf：规范化后的TF名称，用于与RNA基因名匹配。
- tf_in_rna：是否有RNA表达支持；TRUE表示该Motif对应TF在目标细胞群中表达比例达到10%这个阈值。
5.4 TF 活性与 DiffPeak 富集 Motif 重叠分析（韦恩图）
要确定一个 TF 是否真正在当前细胞状态中发挥了关键调控作用，仅靠单一维度的分析往往是不够的。本步骤通过韦恩图，将前面分析的两个独立维度的证据进行交叉验证：
- 维度一：序列富集证据 —— 回答了“该 TF 的结合序列在差异开放的染色质区域中是否异常密集？”
- 维度二：全基因组活性变化证据 —— 回答了“在全基因组范围内，该 TF 潜在结合位点的整体开放程度是否发生了显著改变？”
通过对比这两个集合，我们可以将大量的候选 TF 进行过滤。韦恩图中间的交集部分，代表了那些既在局部差异区域被招募，又在全基因组层面表现出活跃状态的 TF 。
💡 解读结论
位于韦恩图交集区域的 Motif，具有最高级别的置信度，它们是核心主导 TF。在设计后续的分子生物学机制验证（如 ChIP-seq、CRISPR 敲除）或构建上游核心调控网络时，优先将这些交集 TF 作为第一优先级的候选靶点。
R
library(VennDiagram)
library(grid)
motif_dir <- file.path(output, "DIFF", "02_motif")
venn_dir <- file.path(motif_dir, "motif_overlap_venn")
dir.create(venn_dir, recursive = TRUE, showWarnings = FALSE)
all_stat <- list()
k <- 0
if (group_comparison == "") {
tf_file <- file.path(motif_dir, "All_Clusters_TF_Markers.xls")
if (!file.exists(tf_file)) {
warning(paste("跳过: 找不到文件", tf_file))
} else {
tf_df <- read.delim(
tf_file,
sep = "\t", check.names = FALSE, stringsAsFactors = FALSE
)
for (cl in as.character(cluster)) {
cl_s <- gsub(" ", "", cl)
motif_file <- file.path(motif_dir, paste0("c_", cl_s, "_diffPeak_motifs.xls"))
if (!file.exists(motif_file)) {
warning(paste("跳过: 找不到文件", motif_file))
next
}
motif_df <- read.delim(
motif_file,
sep = "\t", check.names = FALSE, stringsAsFactors = FALSE
)
tf_col <- intersect(c("gene", "feature", "motif", "motif_id", "motif.name", "name"), colnames(tf_df))[1]
motif_col <- intersect(c("motif", "motif_id", "motif.name", "feature", "gene", "name"), colnames(motif_df))[1]
tf_raw <- trimws(as.character(tf_df[tf_df$cluster == cl, tf_col, drop = TRUE]))
motif_raw <- trimws(as.character(motif_df[[motif_col]]))
tf_set <- unique(tf_raw[nzchar(tf_raw)])
motif_set <- unique(motif_raw[nzchar(motif_raw)])
inter <- intersect(tf_set, motif_set)
jac <- if (length(union(tf_set, motif_set)) == 0) NA_real_ else length(inter) / length(union(tf_set, motif_set))
out_prefix <- file.path(venn_dir, paste0("c_", cl_s, "_TFvsDiffPeak_motif_overlap_venn"))
title_txt <- paste0("Cluster: ", cl)
write.table(
data.frame(
label = title_txt,
TF_marker_n = length(tf_set),
DiffPeak_motif_n = length(motif_set),
overlap_n = length(inter),
jaccard = jac,
stringsAsFactors = FALSE
),
paste0(out_prefix, "_summary.xls"),
sep = "\t", quote = FALSE, row.names = FALSE
)
write.table(
data.frame(Motif = sort(inter), stringsAsFactors = FALSE),
paste0(out_prefix, "_overlap_motifs.xls"),
sep = "\t", quote = FALSE, row.names = FALSE
)
out_file <- paste0(out_prefix, ".pdf")
grDevices::pdf(out_file, width = 6, height = 6, onefile = TRUE)
grid::grid.newpage()
if (length(tf_set) == 0 && length(motif_set) == 0) {
grid::grid.text(paste0(title_txt, "\n无可用motif"), x = 0.5, y = 0.5)
} else {
venn_grob <- VennDiagram::draw.pairwise.venn(
area1 = length(tf_set),
area2 = length(motif_set),
cross.area = length(inter),
category = c("Activity markers", "Enriched motifs"),
fill = c("#60a5fa", "#f59e0b"),
inverted = length(tf_set) < length(motif_set),
alpha = c(0.5, 0.5),
cex = 1.2,
cat.cex = 1.0,
scaled = FALSE,
ind = FALSE,
cat.pos = c(-30, 30),
cat.dist = c(-0.05, -0.05)
)
grid::grid.draw(venn_grob)
grid::grid.text(title_txt, x = 0.5, y = 0.96, gp = grid::gpar(fontface = "bold"))
}
grDevices::dev.off()
k <- k + 1
all_stat[[k]] <- data.frame(mode = "cluster", cluster = cl, group = "", stringsAsFactors = FALSE)
}
}
} else {
for (cl in as.character(cluster)) {
cl_s <- gsub(" ", "", cl)
tf_file <- file.path(motif_dir, paste0("c_", cl_s, "_TF_Markers.xls"))
if (!file.exists(tf_file)) {
warning(paste("跳过: 找不到文件", tf_file))
next
}
tf_df <- read.delim(
tf_file,
sep = "\t", check.names = FALSE, stringsAsFactors = FALSE
)
for (g in as.character(c(case_group, control_group))) {
g_s <- gsub(" ", "", g)
motif_file <- file.path(motif_dir, paste0("c_", cl_s, "_", g_s, "_diffPeak_motifs.xls"))
if (!file.exists(motif_file)) {
warning(paste("跳过: 找不到文件", motif_file))
next
}
motif_df <- read.delim(
motif_file,
sep = "\t", check.names = FALSE, stringsAsFactors = FALSE
)
tf_col <- intersect(c("gene", "feature", "motif", "motif_id", "motif.name", "name"), colnames(tf_df))[1]
motif_col <- intersect(c("motif", "motif_id", "motif.name", "feature", "gene", "name"), colnames(motif_df))[1]
tf_raw <- trimws(as.character(tf_df[tf_df$group == g, tf_col, drop = TRUE]))
motif_raw <- trimws(as.character(motif_df[[motif_col]]))
tf_set <- unique(tf_raw[nzchar(tf_raw)])
motif_set <- unique(motif_raw[nzchar(motif_raw)])
inter <- intersect(tf_set, motif_set)
jac <- if (length(union(tf_set, motif_set)) == 0) NA_real_ else length(inter) / length(union(tf_set, motif_set))
out_prefix <- file.path(venn_dir, paste0("c_", cl_s, "_", g_s, "_TFvsDiffPeak_motif_overlap_venn"))
title_txt <- paste0("Cluster: ", cl, " | Group: ", g)
write.table(
data.frame(
label = title_txt,
TF_marker_n = length(tf_set),
DiffPeak_motif_n = length(motif_set),
overlap_n = length(inter),
jaccard = jac,
stringsAsFactors = FALSE
),
paste0(out_prefix, "_summary.xls"),
sep = "\t", quote = FALSE, row.names = FALSE
)
write.table(
data.frame(Motif = sort(inter), stringsAsFactors = FALSE),
paste0(out_prefix, "_overlap_motifs.xls"),
sep = "\t", quote = FALSE, row.names = FALSE
)
out_file <- paste0(out_prefix, ".pdf")
grDevices::pdf(out_file, width = 6, height = 6, onefile = TRUE)
grid::grid.newpage()
if (length(tf_set) == 0 && length(motif_set) == 0) {
grid::grid.text(paste0(title_txt, "\n无可用motif"), x = 0.5, y = 0.5)
} else {
venn_grob <- VennDiagram::draw.pairwise.venn(
area1 = length(tf_set),
area2 = length(motif_set),
cross.area = length(inter),
category = c("Activity markers", "Enriched motifs"),
fill = c("#60a5fa", "#f59e0b"),
inverted = length(tf_set) < length(motif_set),
alpha = c(0.5, 0.5),
cex = 1.2,
cat.cex = 1.0,
scaled = FALSE,
ind = FALSE,
cat.pos = c(-30, 30),
cat.dist = c(-0.05, -0.05)
)
grid::grid.draw(venn_grob)
grid::grid.text(title_txt, x = 0.5, y = 0.96, gp = grid::gpar(fontface = "bold"))
}
grDevices::dev.off()
k <- k + 1
all_stat[[k]] <- data.frame(mode = "group", cluster = cl, group = g, stringsAsFactors = FALSE)
}
}
}
if (length(all_stat) > 0) {
write.table(
do.call(rbind, all_stat),
file.path(venn_dir, "venn_generated_items.xls"),
sep = "\t", quote = FALSE, row.names = FALSE
)
message("完成：韦恩图与overlap结果输出到 ", venn_dir)
} else {
message("未生成任何韦恩图结果，请检查对应的输入文件是否存在。")
}R
options(repr.plot.width = 10, repr.plot.height = 7)
grid::grid.draw(venn_grob)
💡 说明
这里的韦恩图示例展示了两种不同计算维度的转录因子 (TF) 的重叠情况，旨在寻找高置信度的核心调控因子：
- Activity markers (蓝色区域)：基于 ChromVAR 计算的 TF 活性矩阵，通过差异分析 (
wilcoxauc) 直接鉴定出的在当前细胞群或实验组中活性显著升高的转录因子。- Enriched motifs (橙色区域)：基于当前细胞群或实验组的特异性开放区域 (差异 Peaks)，使用
FindMotifs进行序列富集分析得到的转录因子。- 交集区域 (Overlap)：同时满足以上两个条件的转录因子。这代表这些 TF 不仅其结合基序在特异性开放区域中被大量发现，且在整个基因组范围内的总体结合活性也是显著上调的。交集中的 TF 极有可能是驱动当前细胞状态或响应实验条件的核心调控因子。