scATAC-seq 差异可及性与功能富集分析教程

作者: SeekGene

时长: 32 分钟

字数: 7.5k 字

更新: 2026-02-27

阅读: 0 次

ATAC + RNA 双组学 Notebooks 分析指南差异富集分析

背景介绍

ATAC_DiffEnrich 是面向单细胞多组学数据中 scATAC-seq 的差异可及性与富集分析流程，旨在从不同细胞群或条件之间系统地鉴定差异可及染色质区域（DARs），解析其潜在调控因子（TF）及受影响的生物学过程/信号通路。该流程基于 Signac/Seurat 的 ATAC 工作流，结合 presto 的快速差异检验、JASPAR/TFBSTools 的 motif 资源、chromVAR 的 motif 活性评估，以及 clusterProfiler 的 GO/KEGG 富集分析，提供从“差异峰 → 调控因子 → 功能通路”的一体化解读。

核心原理

差异可及性检测（DARs）
- 在细胞群（或分组）间采用 Wilcoxon 秩和检验（presto::wilcoxauc）评估每个峰的可及性差异，并进行多重检验校正，获得显著 DAR 集合。
峰-基因连接（Peak-to-Gene）
- 以 ClosestFeature 将显著差异峰关联至邻近基因（如 TSS 附近），用于后续功能富集与通路解释。
转录因子推断（Motif/ChromVAR）
- 基于 JASPAR motif 集对差异峰做富集分析（FindMotifs），并用 chromVAR 量化全局 motif 偏好与细胞层面的 motif 活性变化。
功能与通路富集（GO/KEGG）
- 将峰邻近基因映射到人/鼠注释库，使用 clusterProfiler 进行 GO/KEGG 富集，刻画受影响的生物过程与信号通路。

目的与意义

细胞类型/状态特异调控图谱：识别在特定细胞群中开放的调控元件与核心 TF。
机制层解析：从差异峰—motif—chromVAR 活性—GO/KEGG 多层证据支持调控网络假设。
条件对比与干预线索：在病例-对照或处理-未处理分组中定位被扰动的调控通路与候选靶点。
结果可解释性提升：将无基因注释的峰信号转译为可读的基因与通路层结论。

本教程涵盖内容

差异 Peaks 分析：按细胞群或分组比较，输出显著 DAR 表格与可视化。
Top peaks 展示：对 LogFC 靠前的代表性峰绘制 CoveragePlot。
Motif 富集与活性：差异峰的 motif 富集、motif 序列图与 chromVAR UMAP 展示。
峰邻近基因：输出显著差异峰附近基因表。
GO/KEGG 富集：对峰邻近基因进行富集，生成柱状图/气泡图与汇总表。
交互式结果浏览：显著 DAR、motif、富集结果支持筛选与导出。

预期分析结果

显著差异 Peaks 列表：按细胞群汇总的显著 DAR 及其统计指标。
代表性信号图：Top 差异峰的覆盖度图，用于直观评估峰级别差异。
Motif 与 TF 线索：富集的 TF motif、对应序列图与细胞层 motif 活性变化。
功能通路概览：GO/KEGG 富集结果与可视化，指向受影响的关键生物学过程。
一站式输出目录：包含差异表、峰-基因表、motif/富集图与综合汇总文件，便于下游报告与复用。

参数设置

请先挂载将要分析的云平台相关流程数据，挂载方式请参照 jupyter 使用教程

具体参数填写：

rds：挂载的流程相关的 rds 数据，在/home/{UserName}/workspace/project/{UserName}/目录下，如下列项目数据/home/shumeng/workspace/data/AY1752167131383/
meta：与 rds 同目录下的 metadata 文件
outdir: 分析结果保存路径
clusters_col：细胞类型列名，选择要分析的细胞类型或者聚类对应的标签。假如想对注释好的细胞类型进行分析，则选择其对应的标签，例如 CellAnnotation，与<细胞类型>配合使用。
celltypes：细胞类型，多选，选择要分析的细胞类型或者聚类结果，如 Monocyte 和 Macrophage 等。
species：物种信息
group_comparison：分组比较，选择要进行比较分析的组别名，如要对 Ctrl 组和 STIM 组的基因活力进行比较，则选择 Ctrl 组和 STIM 组对应的标签名，与下面的<case_group：>和<control_group：>搭配适用
case_group：比较组，如上述的 STIM 组。
control_group：对照组，如上述的 Ctrl 组。

# Parameters
outdir = "/home/shumeng/workspace/project/shumeng/"
rds = "/home/shumeng/workspace/data/AY1752167131383/input.rds"
meta = "/home/shumeng/workspace/data/AY1752167131383/meta.tsv"
species = "mouse"
clusters_col = "wknn_res.0.5_d30_l2_50"
celltypes = "0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18"
group_comparison = ""
case_group = ""
control_group = ""

环境准备

R 包加载

请选择 common_r 这个环境进行该整合教程的学习

suppressPackageStartupMessages(suppressWarnings({
    library(Signac)
    library(Seurat)
    library(JASPAR2020)
    library(TFBSTools)
    library(BSgenome.Hsapiens.UCSC.hg38)
    library(BSgenome.Mmusculus.UCSC.mm10)
    library(patchwork)
    library(ggplot2)
    library(presto)
    library(clusterProfiler)
    library(stringr)
    library(tibble)
    library(dplyr)
    library(DOSE)
    library(org.Hs.eg.db)
    library(org.Mm.eg.db)
    library(enrichplot)
    library(foreach)
    library(doParallel)
    library(base64enc) 
    library(DT) 
    library(KEGG.db)
}))

#创建必要的输出目录
output <- paste0(outdir,'/output')
celltypes <- strsplit(celltypes,",")[[1]]
# pct <- as.numeric(pct)
# logfc <- as.numeric(logfc)
# pval_adj <- as.numeric(pval_adj)
# top_n <- as.numeric(top_n)

cl <- makeCluster(min(detectCores(), 25))
registerDoParallel(cl)

options(future.globals.maxSize = 10000 * 1024^2)

dir.create(paste0(output,'/DIFF/01_diffPeak'), recursive = TRUE)
dir.create(paste0(output,'/DIFF/02_topPeak'), recursive = TRUE)
dir.create(paste0(output,'/DIFF/03_motif'), recursive = TRUE)
dir.create(paste0(output,'/DIFF/04_closestFeature'), recursive = TRUE)
dir.create(paste0(output,'/clusterProfiler/GO'), recursive = TRUE)
dir.create(paste0(output,'/clusterProfiler/KEGG'), recursive = TRUE)

#数据读取
obj <- readRDS(rds)
if (meta != ""){
  meta <- read.table(meta,header=T,sep='\t',check.names=F)
  rownames(meta) <- meta$barcodes
  obj <- AddMetaData(obj, meta)
}
obj <- subset(obj, subset = !!sym(clusters_col) == celltypes)

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(dirname(dirname(outdir)), '/', filename)
}

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

Motif 资源与 chromVAR 初始化

Motif 分析用于从差异可及性区域（DARs）推断潜在调控因子（TF）及其活性变化。本流程在显著差异峰上进行 motif 富集（FindMotifs），并结合 chromVAR 计算的 motif 偏离分数（deviation）进行细胞层面的活性可视化，从而在“峰级差异 → TF 线索 → 细胞活性”三层面形成证据闭环。

按物种加载 JASPAR2020 的 PFM 矩阵，将 motif 映射到峰集，并使用对应参考基因组运行 chromVAR；同时设置后续富集与注释所需的物种与数据库标识。

功能概述
- 获取人/鼠的 JASPAR PFM 集：getMatrixSet(JASPAR2020, opts = list(species = 9606/10090))
- 将 motif 写入对象：AddMotifs(object = obj, genome = BSgenome.* , pfm = pfm)
- 计算 motif 偏离分数：RunChromVAR(object = obj, genome = BSgenome.*)
- 统一下游物种标识：
  - human → species = "hsa", db = "org.Hs.eg.db", spe = "human", 基因组 hg38
  - mouse → species = "mmu", db = "org.Mm.eg.db", spe = "mouse", 基因组 mm10
使用要点
- 参考基因组需与物种一致（human→BSgenome.Hsapiens.UCSC.hg38，mouse→BSgenome.Mmusculus.UCSC.mm10）

if (species == "human"){
    pfm <- getMatrixSet(
    x = JASPAR2020,
    opts = list(species = 9606, all_versions = FALSE)
    )
    obj <- AddMotifs(
    object = obj,
    genome = BSgenome.Hsapiens.UCSC.hg38,
    # genome = genome,
    pfm = pfm
    )
    obj <- RunChromVAR(
    object = obj,
    genome = BSgenome.Hsapiens.UCSC.hg38
    # genome = genome,
    )
    species = "hsa"
    db = 'org.Hs.eg.db'
    spe = "human"
}else if (species == "mouse"){
    pfm <- getMatrixSet(
    x = JASPAR2020,
    opts = list(species = 10090, all_versions = FALSE)
    )
    obj <- AddMotifs(
    object = obj,
    genome = BSgenome.Mmusculus.UCSC.mm10,
    pfm = pfm
    )
    obj <- RunChromVAR(
    object = obj,
    genome = BSgenome.Mmusculus.UCSC.mm10
    )

    species = "mmu"
    db = 'org.Mm.eg.db'
    spe = "mouse"
}

细胞群间差异 Peaks 分析

差异 peaks 分析

为在不同细胞群（或分组条件）间识别差异可及性区域（DARs），本流程采用 Wilcoxon 秩和检验进行差异可及性（DA）评估，并使用 Benjamini–Hochberg（BH）进行多重检验校正。presto::wilcoxauc 面向 Seurat 对象的高效实现，可在大规模 scATAC 数据上显著提速且保持统计稳健性。

if (group_comparison == "") {
    da_peaks <- wilcoxauc(obj, clusters_col, assay='data', seurat_assay='ATAC')
    colnames(da_peaks)[2] <- 'cluster'
    # write.table(da_peaks, 
    #            file=paste0(output,'/DIFF/all_diffPeak.xls'),
    #            quote=F, row.names=F, col.names=T, sep='\t')
} else {
    all_da_peaks <- list()
    obj <- subset(obj, subset = !!sym(group_comparison) == c(case_group, control_group))
    
    for (cell_type in unique(obj@meta.data[[clusters_col]])) {
        tryCatch({
            sub_obj <- subset(obj, subset = !!sym(clusters_col) == cell_type)
            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)
            all_da_peaks[[as.character(cell_type)]] <- da_peaks
            # write.table(da_peaks,
            #            file = paste0(output, '/DIFF/by_celltype/', cell_type, '_diffPeak.xls'),
            #            quote = F, row.names = F, col.names = T, sep = '\t')
        }, error = function(e) {
        message(sprintf("Error in cluster %s: %s", cell_type, 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')
}

差异 peaks 结果保存与可视化

下面基于差异可及性结果，对每个细胞群（cluster）进行显著峰筛选、峰-基因映射、motif 富集与活性展示，以及代表性 Top peaks 的覆盖度图绘制，直观呈现峰级差异与潜在转录调控线索。

核心流程
- 集群排序：使用通用排序函数 sort_clusters，优先按数值顺序排列（兼容混合字符编号）。
- 显著峰筛选：对每个 cluster 选取 padj < 0.05 的差异 peaks（feature 列），并导出结果表：
  - output/DIFF/01_diffPeak/c_cluster_diffPeak_significant.xls
- 峰-基因映射：ClosestFeature(obj, regions = peaks) 将显著 peaks 映射至邻近基因：
  - output/DIFF/04_closestFeature/c_cluster_diffPeak_sig_genes.xls
Motif 富集与展示：
- FindMotifs(object = obj, features = peaks) 输出富集表（含 p.adjust 与 fold.enrichment）：
  - output/DIFF/03_motif/c_cluster_diffPeak_motifs.xls
- 依据 fold.enrichment 取 Top 6 motif，MotifPlot 绘制序列图：
  - ..._diffPeak_motifs_top.pdf/png
Top peaks 覆盖度图（CoveragePlot）：
- 按 logFC 降序选取 Top 6 差异 peaks，窗口上下游各扩展 5 kb：
  - 无分组对比：整体绘制；
  - 有分组（group_comparison 非空）：按 case_group vs control_group 分面绘制；
- 输出：output/DIFF/02_topPeak/c_cluster_diffPeak_top.pdf/png
- chromVAR 活性 UMAP：
  - 切换至 chromvar assay，使用 FeaturePlot 展示 Top motifs 的偏离分数（ncol=6/3）：
    - output/DIFF/03_motif/c_cluster_diffPeak_motifs_top_umap.pdf/png

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, ]
        # write.table(data,paste0(output,'/DIFF/c_',gsub(" ","",i),'_diffPeak.xls'),quote=F,row.names=F,col.names=T,sep='\t')
        # data.sig <- data[data$p_val_adj < 0.05, ]
        data.sig <- data[data$padj < 0.05, ]
        write.table(data.sig,paste0(output,'/DIFF/01_diffPeak/c_',gsub(" ","",i),'_diffPeak_significant.xls'),quote=F,row.names=F,col.names=T,sep='\t')
        # peaks <- rownames(data.sig)
        peaks <- data.sig$feature
        
        genes <- ClosestFeature(obj, regions = peaks)
        genes$cluster <- as.character(i)
        write.table(genes,paste0(output,'/DIFF/04_closestFeature/c_',gsub(" ","",i),'_diffPeak_sig_genes.xls'),quote=F,row.names=F,col.names=T,sep='\t')
    
        enriched.motifs <- FindMotifs(object = obj, features = peaks)
        enriched.motifs$cluster <- as.character(i)
        write.table(enriched.motifs,paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs.xls'),quote=F,row.names=F,col.names=T,sep='\t')
    
        motif.top <- head(enriched.motifs[order(enriched.motifs$fold.enrichment, decreasing = TRUE), ])
        p1 <- MotifPlot(object = obj, motifs = rownames(motif.top), nrow = 1)
        p11 <- MotifPlot(object = obj, motifs = rownames(motif.top), nrow = 2)
    
        ggsave(paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top.pdf'), p11, width = 12, height = 6)
        # ggsave(paste0(output,'/DIFF/c_',gsub(" ","",i),'_diffPeak_motifs_top.png'), p1, width = 8, height = 4, dpi = 100, bg = "white")
        ggsave(paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top.png'), p1, width = 24, height = 3, dpi = 100, bg = "white")
    
        # peaks.top <- rownames(head(data.sig[order(data.sig$avg_log2FC, decreasing = TRUE), ]))
        peaks.top <- head(data.sig[order(data.sig$logFC, decreasing = TRUE), ])[['feature']]
        if (group_comparison == "") {
          p2 <- CoveragePlot(object = obj, region = peaks.top, extend.upstream = 5000, extend.downstream = 5000, nrow = 1)
          # p2 <- CoveragePlot(object = obj, region = peaks.top, extend.upstream = 5000, extend.downstream = 5000, nrow = 3)
          p22 <- CoveragePlot(object = obj, region = peaks.top, extend.upstream = 5000, extend.downstream = 5000, nrow = 2)
        } else {
          p2 <- CoveragePlot(object = obj, region = peaks.top, extend.upstream = 5000, extend.downstream = 5000, nrow = 1, group.by = group_comparison, idents = c(case_group, control_group))
          p22 <- CoveragePlot(object = obj, region = peaks.top, extend.upstream = 5000, extend.downstream = 5000, nrow = 2, group.by = group_comparison, idents = c(case_group, control_group))
        }
        ggsave(paste0(output,'/DIFF/02_topPeak/c_',gsub(" ","",i),'_diffPeak_top.pdf'), p22, width = 12, height = 6)
        # ggsave(paste0(output,'/DIFF/c_',gsub(" ","",i),'_diffPeak_top.png'), p2, width = 8, height = 4, dpi = 100, bg = "white")
        ggsave(paste0(output,'/DIFF/02_topPeak/c_',gsub(" ","",i),'_diffPeak_top.png'), p2, width = 24, height = 3, dpi = 100, bg = "white")
        
        DefaultAssay(obj) <- 'chromvar'
        p3 <- FeaturePlot(object = obj, features = rownames(motif.top), ncol = 6)
        # p3 <- FeaturePlot(object = obj, features = rownames(motif.top), ncol = 2)
        p33 <- FeaturePlot(object = obj, features = rownames(motif.top), ncol = 3)
        
        ggsave(paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top_umap.pdf'), p33, width = 12, height = 6)
        # ggsave(paste0(output,'/DIFF/c_',gsub(" ","",i),'_diffPeak_motifs_top.png'), p1, width = 8, height = 4, dpi = 100, bg = "white")
        ggsave(paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top_umap.png'), p3, width = 24, height = 3, dpi = 100, bg = "white")
    }, error = function(e) {
        message(sprintf("Error in cluster %s: %s", i, e$message))
    })
} 
saveRDS(obj,paste0(output, '/output.rds'))
write.table(da_peaks[da_peaks$padj < 0.05, ],paste0(output,'/DIFF/01_diffPeak/all_diffPeak_significant.xls'),quote=F,row.names=F,col.names=T,sep='\t')

细胞群间显著差异 Peaks 列表

计算得到的各细胞群间校正 p 值<0.05 的差异 Peaks 表格。

指标解读（差异 Peaks 表）

feature: 峰位点坐标，格式为 chr:start-end，用于唯一标识每个 ATAC 峰。
cluster: 目标细胞群（或标签）名称，表示该行统计是以此群体为“in-group”与其余细胞为“out-group”进行比较。
avgExpr: 目标细胞群内该峰的平均可及性信号（来源于 ATAC data 矩阵的归一化值，如 TF-IDF/归一化计数），数值越大表示该群体整体更开放。
logFC: 目标群体相对其余细胞的对数倍数变化，>0 表示在该 cluster 中更开放，数值越大差异越明显。
statistic: Wilcoxon 秩和检验的统计量（秩和），用于计算 p 值。
auc: AUC（Area Under Curve）衡量两组可分性，0.5 表示无差异；通常 0.6≈弱、0.7≈中等、≥0.8≈较强分离。
pval: 未校正 p 值，由 Wilcoxon 检验得到。
padj: 多重检验校正后的 p 值（Benjamini–Hochberg），用于显著性判断（常用阈值 < 0.05）。
pct_in: 目标细胞群中该峰为非零/可及的细胞占比（%），反映在群体内的“普遍性”。
pct_out: 其余细胞中该峰为非零/可及的细胞占比（%），用于与 pct_in 对比评估特异性。

options(warn = -1)          
options(message = FALSE)  

# da_peaks <- da_peaks %>% mutate(across(where(is.numeric), ~round(., 2)))

datatable(da_peaks[da_peaks$padj < 0.05, ], rownames = FALSE, filter = 'top', options = list(pageLength = 5, dom = 'Blfrtip', buttons = c('csv', 'excel', 'pdf')), extensions = 'Buttons') %>%
    formatSignif(  
        columns = names(da_peaks)[sapply(da_peaks, is.numeric)],
        digit = 3, zero.print = NA
    )

top peaks 展示

使用 CoveragePlot 对显著差异峰（padj < 0.05）中按 logFC 降序的 Top 6 进行可视化。每个面板展示该峰附近 ±5 kb 的归一化覆盖度曲线（Normalized signal），下方同步显示基因结构与预测 peaks，便于联读。

图片解读
- 上半部分彩色轨迹代表不同 cluster（或细胞类型）的覆盖度，轨迹越高表示该区域越开放。
- 横轴为基因组坐标；下方蓝色箭头为基因坐标（如 Ppp1r16a、Cpt、Mfsd3），灰色条为 peaks 位置。
- 同一峰在多条轨迹中的高度差异，揭示其跨 cluster 的开放度差异；目标组/cluster 曲线整体更高提示更强的染色质开放与潜在调控活性。

# 2. Peaks图表

peaks_links <- data.frame(
    Cluster = cluster,
    Preview = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/DIFF/02_topPeak/c_',gsub(" ","",i),'_diffPeak_top.png')
        file_path <- paste0('../output/DIFF/02_topPeak/c_', gsub(" ","", i), '_diffPeak_top.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    })
)

options(warn = -1)          
options(message = FALSE) 

# Peaks表格
datatable(peaks_links, filter = 'top', 
          escape = FALSE,
          options = list(
              pageLength = 1,
              dom = 'tp'
          ),
          caption = "Peak Coverage Plots",
          rownames = FALSE)

Top peaks 的 Motif 展示

前面已经对每个 cluster 的显著差异峰集合（padj < 0.05）上通过执行 FindMotifs，评估 JASPAR Motif 的富集，识别可能驱动该峰集开放的转录因子（TF），下面分别用并对 Top Motif 进行可视化与活性展示。

Motif 列表

差异 peak 富集 Motif 列表指标解读

motif：Motif 的唯一标识符（如 JASPAR ID）。
observed：在差异 peaks 区域中实际检测到该 motif 的次数。
background：在背景区域（非差异 peaks）中检测到该 motif 的次数。
percent.observed：该 motif 在差异 peaks 区域中的出现比例（%）。
percent.background：该 motif 在背景区域中的出现比例（%）。
fold.enrichment：富集倍数，表示 motif 在差异 peaks 区域的富集程度（percent.observed / percent.background），数值越大富集越显著。
pvalue：富集检验的原始 p 值，衡量 motif 富集的统计显著性。
motif.name：Motif 对应的转录因子名称或注释。
p.adjust：多重检验校正后的 p 值（如 FDR），常用阈值 < 0.05 判定显著富集。
cluster：对应的细胞群/cluster，表明该 motif 富集于哪个 cluster 的差异 peaks。

data_dir <- paste0(output,"/DIFF/03_motif")
xls_files <- list.files(path = data_dir, 
                       pattern = "motifs\\.xls$", 
                       full.names = TRUE)

merged_data <- data.frame()
for (file in xls_files) {
    df <- read.table(file, sep="\t", header=TRUE)
    df$cluster <- as.character(df$cluster)
    merged_data <- bind_rows(merged_data, df)
}
merged_data$cluster <- as.character(merged_data$cluster)

options(warn = -1)          
options(message = FALSE)  

# merged_data_motif <- merged_data %>% mutate(across(where(is.numeric), ~round(., 2)))
merged_data_motif <- merged_data 
datatable(merged_data_motif[merged_data_motif$p.adjust < 0.05, ], rownames = FALSE, filter = 'top', options = list(pageLength = 5, dom = 'Blfrtip', buttons = c('csv', 'excel', 'pdf')), extensions = 'Buttons') %>%
    formatSignif(  
        columns = names(merged_data_motif)[sapply(merged_data_motif, is.numeric)],
        digit = 3, zero.print = NA
    )

Motif 可视化（chromVAR UMAP）

按 cluster 展示 Top Motif 的序列 logo ，直观反映其碱基组成和保守性。
按 cluster 展示 Top Motif 的细胞层活性分布预览。每个小图为一个 Motif 的 chromVAR 偏离分数（deviation z-score）在 UMAP 空间的投影。

# 1. Motifs图表

motifs_links1 <- data.frame(
    Cluster = cluster,
    Preview = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top.png')
        file_path <- paste0('../output/DIFF/03_motif/c_', gsub(" ","", i), '_diffPeak_motifs_top.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    })
)

motifs_links2 <- data.frame(
    Cluster = cluster,
    Preview = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/DIFF/03_motif/c_',gsub(" ","",i),'_diffPeak_motifs_top_umap.png')
        file_path <- paste0('../output/DIFF/03_motif/c_', gsub(" ","", i), '_diffPeak_motifs_top_umap.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    })
)

options(warn = -1)          
options(message = FALSE)  

# Motifs表格
datatable(motifs_links1, filter = 'top', 
          escape = FALSE,
          options = list(
              pageLength = 1,
              dom = 'tp'
          ),
          caption = "Motifs Top Plots",
          rownames = FALSE)

datatable(motifs_links2, filter = 'top', 
          escape = FALSE,
          options = list(
              pageLength = 1,
              dom = 'tp'
          ),
          rownames = FALSE)

差异 peaks 附近的基因

单独的峰位点坐标难以直接解释其功能。我们使用 ClosestFeature 将显著差异 peaks（padj < 0.05）映射到最近的基因，从而推断可能受这些调控元件影响的候选基因，并为后续 GO/KEGG 富集分析提供输入。

指标解读（ClosestFeature 结果表）

tx_id: 最近转录本的 Ensembl 转录本 ID。
gene_name: 最近基因的基因符号（SYMBOL）。
gene_id: 最近基因的 Ensembl 基因 ID。
gene_biotype: 基因类型（如 protein_coding、lincRNA 等）。
type: 距离最近的功能注释类型（如 exon、UTR、intron、promoter 等）。
closest_region: 最近功能注释片段的基因组坐标。
query_region: 输入的差异 peak 坐标。
distance: 两者之间的距离（bp）。0 表示重叠；数值越小，调控关联可能性越高；常将 ≤2000 bp 视为启动子近端候选。
cluster: 该差异 peak 所属的细胞群标签。

data_dir <- paste0(output,"/DIFF/04_closestFeature")
xls_files <- list.files(path = data_dir, 
                       pattern = "diffPeak_sig_genes\\.xls$", 
                       full.names = TRUE)

merged_data <- data.frame()
for (file in xls_files) {
    df <- read.table(file, sep="\t", header=TRUE)
    df$cluster <- as.character(df$cluster)
    merged_data <- bind_rows(merged_data, df)
}
merged_data$cluster <- as.character(merged_data$cluster)

options(warn = -1)          
options(message = FALSE) 

# merged_data_genes <- merged_data %>% mutate(across(where(is.numeric), ~round(., 2)))
merged_data_genes <- merged_data 
datatable(merged_data_genes, rownames = FALSE, filter = 'top', options = list(pageLength = 5, dom = 'Blfrtip', buttons = c('csv', 'excel', 'pdf')), extensions = 'Buttons') %>%
    formatSignif(  
        columns = names(merged_data_genes)[sapply(merged_data_genes, is.numeric)],
        digit = 3, zero.print = NA
    )

差异 peaks 附近基因的功能富集

将 ClosestFeature 得到的邻近基因作为输入，使用 clusterProfiler 对每个 cluster 分别进行 GO 与 KEGG 富集，系统刻画差异调控区域可能影响的生物学过程与信号通路。

分析流程
- 基因准备：按 cluster 提取 gene_name，用 bitr 映射到 ENSEMBL/ENTREZID。
GO 富集：enrichGO(ont = "ALL", pAdjustMethod = "BH")，输出显著条目与可视化。
KEGG 富集：enrichKEGG(organism = hsa/mmu)，并可视化。

go_enrich <- function(eg,db,outdir,prefix) {
    genelist <- eg$SYMBOL
    go <- enrichGO(genelist, OrgDb=db, ont='ALL',pAdjustMethod = 'BH',qvalueCutoff = 1,pvalueCutoff = 1,keyType = 'SYMBOL')
    go1 <- data.frame(cluster=prefix, go)
    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)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
    pdf(paste0(outdir,'/',prefix,'_GO_dot.pdf',sep=''),width = 10, height = 8)
    print(dotplot(go,showCategory=20)+ scale_y_discrete(labels = function(x) 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)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
    png(paste0(outdir,'/',prefix,'_GO_dot.png',sep=''),width = 780, height = 580)
    print(dotplot(go,showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
}


kegg_enrich <- function(eg,species,outdir,prefix) {
    genenames<-eg$SYMBOL
    names(genenames)<-eg$ENTREZID
    genelist <- eg$ENTREZID
    download.KEGG.Path <- function (species) {
    keggpathid2extid.df <- clusterProfiler:::kegg_link(species, "pathway")
    if (is.null(keggpathid2extid.df)) {
      message <- paste("Failed to download KEGG data.", "Wrong 'species' or the network is unreachable.", 
                     "The 'species' should be one of organisms listed in", 
                     "'https://www.genome.jp/kegg/catalog/org_list.html'")
      stop(message)
    }
    keggpathid2extid.df[, 1] %<>% gsub("[^:]+:", "", .)
    keggpathid2extid.df[, 2] %<>% gsub("[^:]+:", "", .)
    keggpathid2name.df <- clusterProfiler:::kegg_list("pathway")
    keggpathid2name.df[, 1] %<>% gsub("map", species, .)
    keggpathid2extid.df <- keggpathid2extid.df[keggpathid2extid.df[, 1] %in% keggpathid2name.df[, 1], ]
    return(list(KEGGPATHID2EXTID = keggpathid2extid.df, KEGGPATHID2NAME = keggpathid2name.df))
  }
    
    assignInNamespace('download.KEGG.Path',download.KEGG.Path,ns = 'clusterProfiler')
    head(genelist)
    kegg <- enrichKEGG(genelist, organism = species, keyType = 'kegg', pAdjustMethod = 'BH',pvalueCutoff = 1,qvalueCutoff = 1,use_internal_data = TRUE)
    gene_list <- strsplit(kegg$geneID,split='/')
    geneName<-unlist(lapply(gene_list,FUN=function(x){paste(genenames[x],collapse='/')}))
   
    KEGGenrich <- data.frame(cluster=prefix, kegg[,1:8],geneName,kegg[,9,drop=F])
    write.table(KEGGenrich,file=paste(outdir,'/',prefix,'_KEGGenrich.xls',sep=''),sep='\t',quote=F,row.names=F)
    pdf(paste0(outdir,'/',prefix,'_KEGG_dot.pdf',sep=''),width = 8, height = 8)
    print(dotplot(kegg, showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
    pdf(paste0(outdir,'/',prefix,'_KEGG_bar.pdf',sep=''),width = 8, height = 8)
    print(barplot(kegg, showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
    png(paste0(outdir,'/',prefix,'_KEGG_dot.png',sep=''),width = 780, height = 580)
    print(dotplot(kegg, showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
    png(paste0(outdir,'/',prefix,'_KEGG_bar.png',sep=''),width = 780, height = 580)
    print(barplot(kegg, showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
    dev.off()
}

files <- list.files(paste0(output,"/DIFF/04_closestFeature"),pattern="diffPeak_sig_genes.xls")

# cores <- detectCores() - 1  
# cl <- makeCluster(cores)
# registerDoParallel(cl)

# foreach(f=files,.packages = c("clusterProfiler", "reshape2", "ggplot2", "stringr", "tibble", "dplyr", "DOSE", "enrichplot")) %dopar% {
for (f in files){
  prefix <- gsub('_diffPeak_sig_genes.xls','',f)
  f <- paste0(outdir,"/output/DIFF/04_closestFeature/",f)
  print(f)
  print(prefix)
  clus_1<-read.table(f,header=T,sep='\t')
  clus<-clus_1$gene_name

  tryCatch({
    eg <- bitr(clus, fromType="SYMBOL", toType=c("ENSEMBL","ENTREZID"), OrgDb=db)
    print(head(eg))
  }, error = function(e) {
    message(paste("处理文件", f, "时出错:", e$message))
  })
  
  tryCatch({
    go_enrich(eg,db,paste0(outdir,"/output/clusterProfiler/GO"),prefix)
  }, error = function(e) {
    message(paste0("GO富集分析错误: ", prefix, " - ", e$message))
  })
  tryCatch({
    kegg_enrich(eg, species, paste0(outdir,"/output/clusterProfiler/KEGG"), prefix)
  }, error = function(e) {
    message(paste0("KEGG富集分析错误: ", prefix, " - ", e$message))
  })

}
# stopCluster(cl)

# 合并 GO 富集分析结果
go_dir <- paste0(output,"/clusterProfiler/GO")
if(dir.exists(go_dir)) {  # 检查目录是否存在
    xls_files <- list.files(path = go_dir, 
                           pattern = "\\.xls$", 
                           full.names = TRUE)
    
    merged_data <- data.frame()
    for (file in xls_files) {
        df <- read.table(file, sep="\t", header=TRUE, fill=TRUE, quote="", check.names=FALSE)
        merged_data <- bind_rows(merged_data, df)
    }
    
    # 直接使用完整路径写入文件
    write.table(merged_data, 
                file = file.path(go_dir, "all_GOenrich.xls"), 
                row.names = FALSE,
                sep="\t",
                quote=FALSE)
}

GO 富集结果

表格汇总各 cluster 的 GO 富集条目与显著性统计，用于刻画差异 peaks 邻近基因的功能指向。

指标解读
- cluster：结果所属的细胞群标签。
- ONTOLOGY：GO 域，BP（生物过程）/ CC（细胞组分）/ MF（分子功能）。
- ID：GO 术语编号（如 GO:0007264）。
- Description：GO 术语名称。
- GeneRatio：命中该术语的输入基因数/输入基因总数（越大说明在本次输入中命中比例越高）。
- BgRatio：该术语在背景基因集中的命中数/背景总数（用于与 GeneRatio 对比评估富集强度）。
- pvalue：富集检验原始 p 值。
- p.adjust：多重检验校正后的 p 值（BH/FDR），常用阈值 < 0.05 判定显著。
- qvalue：q 值估计（FDR 另一种表示），可用阈值 < 0.2 辅助筛选。
- geneID：命中该术语的基因列表，常用 “/” 分隔。

options(warn = -1)          
options(message = FALSE) 

# merged_data_go <- merged_data %>% mutate(across(where(is.numeric), ~round(., 2)))
merged_data_go <- merged_data 
if(nrow(merged_data_go[merged_data_go$p.adjust < 0.05, ]) > 0) {
datatable(merged_data_go[merged_data_go$p.adjust < 0.05, ], rownames = FALSE, filter = 'top', caption = "GO Enrichment",
          options = list(pageLength = 5, dom = 'Blfrtip', buttons = c('csv', 'excel', 'pdf')), 
          extensions = 'Buttons') %>%
    formatSignif(  
        columns = names(merged_data_go)[sapply(merged_data_go, is.numeric)],
        digit = 3, zero.print = NA
    )
}

GO 富集结果可视化

Dot（左图）横轴 GeneRatio：命中该术语的输入基因数/输入总数，越大代表相对更集中。纵轴 GO 术语（BP/CC/MF）。圆点大小 Count：命中基因的绝对数量，越大越多。颜色 p.adjust：颜色越偏红表示校正后显著性越高（p.adjust 越小）。解读：优先关注“红且大、靠右”的点（显著性高、命中多、比例大）。 Bar（右图）横轴 Count：命中基因数；条形越长命中越多。纵轴 GO 术语。颜色 p.adjust：颜色越红表示显著性越高。解读：从上到下对比条形长度与颜色，锁定“长且红”的条目。

# 3. 富集图表

go_links <- data.frame(
    Cluster = cluster,
    GO_dot = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/clusterProfiler/GO/c_', gsub(" ","", i), '_GO_dot.png')
        file_path <- paste0('../output/clusterProfiler/GO/c_', gsub(" ","", i), '_GO_dot.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    }),
    GO_bar = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/clusterProfiler/GO/c_', gsub(" ","", i), '_GO_bar.png')
        file_path <- paste0('../output/clusterProfiler/GO/c_', gsub(" ","", i), '_GO_bar.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    })
)

options(warn = -1)          
options(message = FALSE)  

# Enrich表格
datatable(go_links, filter = 'top', 
          escape = FALSE,
          options = list(
              pageLength = 1,
              dom = 'tp'
          ),
          caption = "Enrich Plots",
          rownames = FALSE)

- 柱状图：x 轴为注释到 GO 编号上的基因数，y 轴为 GO term 描述，颜色代表校正 p 值，颜色越红 p 值越小，颜色越蓝 p 值越大。
- 气泡图：x 轴为 GeneRatio，y 轴为 GO term 描述，颜色代表校正 p 值，颜色越红 p 值越小，颜色越蓝 p 值越大，点的大小表示富集在通路内的差异基因个数。

# 合并 KEGG 富集分析结果
kegg_dir <- paste0(output,"/clusterProfiler/KEGG")
if(dir.exists(kegg_dir)) {  # 检查目录是否存在
    xls_files <- list.files(path = kegg_dir, 
                           pattern = "\\.xls$", 
                           full.names = TRUE)
    
    merged_data <- data.frame()
    for (file in xls_files) {
        df <- read.table(file, sep="\t", header=TRUE, fill=TRUE, quote="", check.names=FALSE)
        merged_data <- bind_rows(merged_data, df)
    }
    
    # 直接使用完整路径写入文件
    write.table(merged_data, 
                file = file.path(kegg_dir, "all_KEGGenrich.xls"), 
                row.names = FALSE,
                sep="\t",
                quote=FALSE)
}

KEGG 富集结果

汇总各 cluster 的 KEGG 通路富集条目，用于刻画差异 peaks 邻近基因可能参与的信号通路与代谢路径（基于 enrichKEGG，物种代码 human→hsa，mouse→mmu）。

指标解读
- cluster：所属细胞群标签。
- category：KEGG 一级分类（如 Metabolism、Cellular Processes）。
- subcategory：KEGG 二级分类（如 Carbohydrate metabolism、Signal transduction）。
- ID：通路编号（如 mmu04015）。
- Description：通路名称。
- GeneRatio：输入基因中命中该通路的比例（命中数/输入总数），越大说明在当前输入中更集中。
- BgRatio：背景基因集中命中该通路的比例（命中数/背景总数），用于与 GeneRatio 对比评估富集强度。
- pvalue：富集检验原始 p 值。
- p.adjust：多重检验校正后的 p 值（FDR/BH），常用阈值 < 0.05 判定显著。
- geneID/geneName：命中该通路的基因列表，用 “/” 分隔。

options(warn = -1)          
options(message = FALSE) 

# merged_data_kegg <- merged_data %>% mutate(across(where(is.numeric), ~round(., 2)))
merged_data_kegg <- merged_data 
if(nrow(merged_data_kegg[merged_data_kegg$p.adjust < 0.05, ]) > 0) {
datatable(merged_data_kegg[merged_data_kegg$p.adjust < 0.05, ], rownames = FALSE, filter = 'top', caption = "KEGG Enrichment",
          options = list(pageLength = 5, dom = 'Blfrtip', buttons = c('csv', 'excel', 'pdf')), 
          extensions = 'Buttons') %>%
    formatSignif(  
        columns = names(merged_data_kegg)[sapply(merged_data_kegg, is.numeric)],
        digit = 3, zero.print = NA
    )
}

KEGG 富集结果可视化

Dot（左图）
- 横轴 GeneRatio：命中该通路的输入基因数/输入总数，越大越集中。
- 纵轴通路名称（含物种代码，如 mmu04015）。
- 点大小 Count：命中基因数，越大命中越多。
- 颜色 p.adjust：颜色越偏红表示显著性越高（校正后 p 值更小）。
- 解读：优先关注“红且大、靠右”的点（显著、命中多、比例高）。
Bar（右图）
- 横轴 Count：命中基因数；条形越长表示命中越多。
- 纵轴通路名称。
- 颜色 p.adjust：越红显著性越高。
- 解读：从上到下比较条形长度与颜色，锁定“长且红”的通路。

# 3. 富集图表

kegg_links <- data.frame(
    Cluster = cluster,
    KEGG_dot = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/clusterProfiler/KEGG/c_', gsub(" ","", i), '_KEGG_dot.png')
        file_path <- paste0('../output/clusterProfiler/KEGG/c_', gsub(" ","", i), '_KEGG_dot.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    }),
    KEGG_bar = sapply(cluster, function(i) {
        file_path_ab <- paste0(output,'/clusterProfiler/KEGG/c_', gsub(" ","", i), '_KEGG_bar.png')
        file_path <- paste0('../output/clusterProfiler/KEGG/c_', gsub(" ","", i), '_KEGG_bar.png')
        if(file.exists(file_path_ab)) {
            sprintf('<a href="%s"><img src="/placeholder.svg" height="200px"></a>', file_path, file_path)
        } else {
            "图片不存在"
        }
    })
)

options(warn = -1)          
options(message = FALSE)  

# Enrich表格
datatable(kegg_links, filter = 'top', 
          escape = FALSE,
          options = list(
              pageLength = 1,
              dom = 'tp'
          ),
          caption = "Enrich Plots",
          rownames = FALSE)

- 柱状图：x 轴为注释到 KEGG 通路上的基因数，y 轴为 KEGG 通路描述，颜色代表校正 p 值，颜色越红 p 值越小，颜色越蓝 p 值越大。
- 气泡图：x 轴为 GeneRatio，y 轴为 KEGG 通路描述，颜色代表校正 p 值，颜色越红 p 值越小，颜色越蓝 p 值越大，点的大小表示富集在通路内的差异基因个数。

结果文件

DIFF/*diffPeak.xls 差异 Peaks 总表
DIFF/*diffPeak_significant.xls 显著差异 Peaks 表格
DIFF/*diffPeak_top.pdf/png top6 的显著差异 Peaks 染色质可及性信号图
DIFF/*diffPeak_motifs.xls 显著差异 Peaks 的 motif 表格
DIFF/*diffPeak_motifs_top.pdf/png top6 的显著差异 Peaks 的 motif 图
DIFF/*diffPeak_sig_genes.xls 显著差异 peaks 附近的基因
clusterProfiler/*/all_*enrich.xls 富集结果总表
clusterProfiler/*/*enrich.xls 各细胞群富集结果
clusterProfiler/*/*bar.pdf/png 各细胞群富集柱状图
clusterProfiler/*/*dot.pdf/png 各细胞群富集气泡图

结果目录：目录下包括此项分析所涉及的所有图片以及表格等结果文件

文献案例解析

文献一：
文献《Single-cell multiomics analysis reveals regulatory programs in clear cell renal cell carcinoma》进行差异可及染色质区域（DAR）相关分析，发现肿瘤细胞的 DAR 在代谢相关的生物过程中显着富集。

参考资料

[1] Stuart, Tim, Avi Srivastava, Caleb Lareau, and Rahul Satija. "Multimodal single-cell chromatin analysis with Signac." BioRxiv (2020): 2020-11.

[2] Korsunsky I, A. Nathan N Millard, and S. Raychaudhuri. "presto: Fast Functions for Differential Expression using Wilcox and AUC." R package. https://rdrr.io/github/immunogenomics/presto (2019).

附录

.csv/.xls/.txt ：结果数据表格文件，文件以逗号或制表符（Tab）分隔。Linux/Mac 用户使用 less 或 more 命令查看；Windows 用户使用高级文本编辑器 Notepad++ 等查看，也可以用 Microsoft Excel 打开。
.png：结果图像文件，位图，无损压缩。
.pdf：结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用 Adobe Illustrator 进行图片编辑，用于文章发表等。

scATAC-seq 差异可及性与功能富集分析教程 ​

背景介绍 ​

核心原理 ​

目的与意义 ​

本教程涵盖内容 ​

预期分析结果 ​

参数设置 ​

环境准备 ​

Motif 资源与 chromVAR 初始化 ​

细胞群间差异 Peaks 分析 ​

差异 peaks 分析 ​

差异 peaks 结果保存与可视化 ​

细胞群间显著差异 Peaks 列表 ​

top peaks 展示 ​

Top peaks 的 Motif 展示 ​

Motif 列表 ​

Motif 可视化（chromVAR UMAP） ​

差异 peaks 附近的基因 ​

差异 peaks 附近基因的功能富集 ​

KEGG 富集结果可视化 ​

结果文件 ​

文献案例解析 ​

参考资料 ​

附录 ​