ATAC + RNA 多组学：差异可及性分析
时长: 49 分钟
字数: 11.2k 字
更新: 2026-07-07
阅读: 0 次
1. 教程简介
本教程基于 Seurat + Signac + clusterProfiler 等单细胞多组学分析工具，用于对 SeekArc 单细胞多组学数据（RNA + ATAC）进行 差异 peak 分析、peak 邻近基因注释及下游功能解读。
该教程面向已完成降维聚类的 SeekArc 多组学数据，通过识别细胞群间或实验条件间的差异开放染色质区域（differential peaks），结合 ClosestFeature() 进行邻近基因注释，将表观基因组层面的染色质开放性变化映射到潜在的靶基因上，实现以下核心分析目标：
- 鉴定细胞类型特异性开放染色质区域：识别不同细胞簇特有的开放 peak，辅助细胞身份定义；
- 挖掘条件特异性调控元件：比较不同处理组或疾病组的差异 peak，揭示环境刺激或病理状态下的染色质重塑事件；
- 差异 peak 的靶基因预测：利用最近邻法则将差异 peak 关联到潜在受调控基因，为后续 RNA 层面的差异表达分析提供表观基因组维度的补充证据；
- 富集分析与生物学通路挖掘：将注释到的邻近基因列表输入
clusterProfiler，进行 GO/KEGG 等功能富集分析，揭示差异开放区域参与的生物学过程。
注意
💡 文档中展示的图表仅为部分代表性的可视化结果。如需查看所有细胞群或对比组的完整分析结果（包含详细的数据表格与所有高清绘图），请前往
../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.db) #人类基因注释数据库
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 多组学对象进行高级差异分析（DiffEnrich）。开始分析前，需要准备以下几类核心输入文件：
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/fpath"
## rds：输入的 Seurat 对象 RDS 文件路径，包含已预处理的单细胞多组学数据
rds = "/path/to/input.rds"
## meta：元数据文件路径（TSV格式），包含样本信息、细胞类型注释等元数据
meta = "/path/to/meta.tsv"
## species：物种信息，"others"表示非标准模式物种（如人/小鼠以外的物种）
species = "human"
## 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"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_topPeak'), recursive = TRUE)
dir.create(paste0(output,'/DIFF/03_closestFeature'), recursive = TRUE)2.2 函数定义
在正式开展分析之前，需要预先定义几个核心的辅助函数，以便在后续的分析循环中进行高效的重复调用。脚本主要定义了以下三个功能函数：
make_pseudobulk函数 目的是将稀疏的单细胞矩阵按指定的分组和大小（pb_size）进行拟批量（pseudobulk）聚合，用于降低单细胞数据的稀疏性，从而提高后续相关性计算及热图可视化的稳健性。go_enrich函数 该函数主要利用clusterProfiler包的enrichGO功能对输入的靶基因（Symbol 格式）进行 GO（Gene Ontology）富集分析。分析涵盖了生物学过程（BP）、细胞组分（CC）和分子功能（MF）三个层面（ont='ALL'）。除了输出完整的富集结果表格（.xls），函数还会自动生成展示 Top 20 显著富集通路的柱状图（barplot）和气泡图（dotplot），并在绘图时自动对过长的通路名称进行换行处理，以保证排版美观。kegg_enrich函数 该函数负责执行 KEGG 通路富集分析。由于 KEGG 数据库要求使用 ENTREZID 作为输入，函数内部会先进行基因 ID 的映射处理，并调用在线 KEGG 数据库（use_internal_data = FALSE）进行富集计算。为了提高结果的易读性，函数内置了结果清洗逻辑：剥离通路描述中冗余的物种后缀信息，并将富集结果中的 ENTREZID 重新转换为易读的 Symbol 基因名，最终导出详尽的结果表格及对应的柱状图和气泡图。
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
#定义一个GO富集分析的函数
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 = 8, height = 8)
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 = 8, height = 8)
print(dotplot(go,showCategory=20)+ scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
dev.off()
}R
#定义一个KEGG富集分析的函数
kegg_enrich <- function(eg, species, outdir, prefix) {
genenames <- eg$SYMBOL
names(genenames) <- eg$ENTREZID
genelist <- eg$ENTREZID
# 使用 KEGG 在线数据库（支持所有物种，包括 bta）
kegg <- enrichKEGG(
gene = genelist,
organism = species, # bta
keyType = "kegg",
pAdjustMethod = "BH",
pvalueCutoff = 1,
qvalueCutoff = 1,
use_internal_data = FALSE
)
# 解析基因名
gene_list <- strsplit(kegg$geneID, split = "/")
geneName <- unlist(lapply(gene_list, function(x) paste(genenames[x], collapse = "/")))
# 清理 KEGG 描述
kegg@result$Description <- sub("\\s+-\\s+[^\\(]+\\s*\\([^)]+\\)\\s*$", "",
as.character(kegg@result$Description))
# 输出表格
KEGGenrich <- data.frame(
cluster = prefix,
kegg[, 1:6],
kegg[, 10:12],
geneName,
kegg[, 14, drop = 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) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
dev.off()
pdf(paste0(outdir, "/", prefix, "_KEGG_bar.pdf"), width = 8, height = 8)
print(barplot(kegg, showCategory = 20) +
scale_y_discrete(labels = function(x) str_wrap(x, width = 50)))
dev.off()
}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)
DefaultAssay(obj) <- "ATAC"
Idents(obj) <- clusters_col💡 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 Path)
在单细胞对象转移或跨服务器读取时，ATAC 模态下的片段文件路径经常会失效。程序会遍历 Seurat 对象中 ATAC Assay 的所有 Fragment 记录，提取文件名，并将其路径重定向到当前的输出目录（fpath），确保后续计算 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. 开始差异 Peak 分析
4.1 细胞群间/组间差异 Peak 鉴定
presto 包的 wilcoxauc 函数做差异peak分析相较于Seurat的FindAllMarkers和FindMarkers而言，更快速和高效，因此本教程利用 Wilcoxon 秩和检验对 ATAC Assay 的开放矩阵进行差异分析。
根据分析目的不同，一般可以进行两种差异分析模式。
- 单聚类 Marker 鉴定模式（Cluster vs All）
- 触发条件：当
group_comparison参数为空（""）时。 - 分析逻辑：程序会将数据集中每一个细胞类型（或聚类）与剩余所有细胞进行对比。这种模式主要用于寻找定义某种特定细胞类型（如 B cells, T cells）的专属开放区域（Cluster-specific Marker Peaks）。
- 触发条件：当
- 细胞类型内组间差异模式（Case vs Control）
- 触发条件：当
group_comparison参数不为空（如设置为"Sample"）时。 - 分析逻辑：程序会首先提取指定对比的两个分组（
case_group和control_group），然后遍历每一个细胞类型，在同一种细胞类型内部比较疾病组（实验组）与对照组的差异。 - 应用场景：这种模式常用于探索特定细胞亚群在疾病发生、药物处理或不同发育阶段下的特异性表观遗传改变。程序会在计算完成后，将所有细胞类型内部的组间差异结果合并为一个总表（
all_diffPeak.xls）并自动导出。
- 触发条件：当
本教程将通过判断 group_comparison 参数是否为空，来选择执行哪种分析模式
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')
}R
head(da_peaks[da_peaks$padj < pval_adj & da_peaks$logFC > logFC, ])| feature | group | avgExpr | logFC | statistic | auc | pval | padj | pct_in | pct_out | cluster | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | |
| T cells.130 | chr1-1692209-1693568 | 25030508_pbmc_1_arc | 0.7942903 | 0.3801981 | 777228.5 | 0.6371901 | 2.159255e-29 | 6.461503e-26 | 50.927835 | 31.8489066 | T cells |
| T cells.1059 | chr1-12510125-12511334 | 25030508_pbmc_1_arc | 0.4425744 | 0.2932666 | 658140.5 | 0.5395589 | 2.695462e-12 | 6.892835e-10 | 11.958763 | 4.2544732 | T cells |
| T cells.2161 | chr1-25774578-25775743 | 25030508_pbmc_1_arc | 0.4612178 | 0.2908842 | 689600.5 | 0.5653506 | 4.535117e-18 | 3.051895e-15 | 20.618557 | 8.2703777 | T cells |
| T cells.4675 | chr1-59814184-59815421 | 25030508_pbmc_1_arc | 0.4529309 | 0.3608838 | 687916.5 | 0.5639700 | 3.876126e-28 | 9.912029e-25 | 16.288660 | 3.8170974 | T cells |
| T cells.4936 | chr1-64876900-64877760 | 25030508_pbmc_1_arc | 0.3995816 | 0.2821482 | 649553.0 | 0.5325187 | 1.865536e-11 | 3.864215e-09 | 9.278351 | 2.9423459 | T cells |
| T cells.5445 | chr1-83508232-83508942 | 25030508_pbmc_1_arc | 0.2853563 | 0.2793229 | 640625.0 | 0.5251993 | 4.262219e-26 | 8.347758e-23 | 5.154639 | 0.1192843 | T cells |
💡 解读指南
上表交互式展示了差异peak的结果表，每一行代表一个差异可及性 Peak，用于回答“哪个组在该 Peak 上更开放、显著性如何、效应量多大”。
- feature：Peak 坐标（如
chr1-4846826-4847624），即开放性具有差异的peak。- group：组别信息，表示该 Peak 在该组更开放，该列只出现在组间差异分析中。
- cluster：细胞群名称，该peak具体在哪个细胞群内差异开放。
- avgExpr：该 Peak 在目标组中的平均可及性水平，是ATAC 信号均值，越大通常表示整体更开放。
- logFC：差异倍数效应量，对差异倍数做log处理后的数值，正值越大，差异倍数越大。
- statistic：Wilcoxon 检验统计量，用于排序与显著性计算。
- auc：区分度指标，0.5 近似无差异，越偏离 0.5 说明区分能力越强。
- pval：原始 P 值，表示差异显著性。
- padj：多重检验校正后的 P 值，通常用它作为显著性判断主指标。
- pct_in：目标组中该 Peak 可检测到信号的细胞百分比。
- pct_out：对照组/其余组中该 Peak 可检测到信号的细胞百分比；与
pct_in对比可直观看到组间覆盖差异。
R
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')4.2 CoveragePlot 可视化差异 Peak
单纯的统计学表格（如 p-value 和 logFC）虽然能找出显著变化的区间，但无法直观反映真实的染色质开放状态。为了验证统计结果的可靠性并直观展示表观遗传差异，我们需要回到原始的测序数据层面进行检查。
下面的代码模块遍历了每一个细胞群（Cluster），针对其显著差异 Peak 进行深度可视化：
- 提取 Top 候选区：在通过严格的显著性阈值过滤后，程序会按效应量（
logFC）进行降序排列，自动提取变化最剧烈的 Top 差异 Peak。 - 绘制染色质覆盖图：调用 Signac 的
CoveragePlot函数，以提取出的 Top Peak 坐标为中心，向上下游各延伸 5000 bp。该图会提取并叠加所选细胞群在实验组与对照组（或其余细胞群）的 Tn5 酶切插入片段信号，生成直观的基因组 Track 图。
通过 CoveragePlot，我们可以清晰地看到不同比较组别在特定基因组区间上“峰型（Peak）”的高低变化，这为我们的差异统计结果提供了最直接的生物学视觉证据。
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)R
# 遍历所有的 Cluster，提取 Top 差异 Peak 进行 CoveragePlot 可视化
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) {
data <- data[valid_rows, ]
}
# 根据 group_comparison 决定分析模式
if (group_comparison == "") {
analysis_groups <- c("cluster_marker")
} else {
analysis_groups <- c(case_group, control_group)
}
for (grp in analysis_groups) {
tryCatch({
# 根据不同模式过滤显著差异 Peak
if (grp == "cluster_marker") {
data.sig <- data[data$padj < pval_adj & data$logFC > logFC, ]
prefix <- paste0('c_', gsub(" ", "", i))
} else {
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) next
# 提取 logFC 最大的 Top 1 差异 Peak（可根据需要修改提取数量）
peaks.top <- head(data.sig[order(data.sig$logFC, decreasing = TRUE), ])[['feature']]
if (length(peaks.top) > 0) {
print(paste0("Drawing CoveragePlot for Top Peak in ", prefix, ": ", peaks.top[1]))
# 绘制 CoveragePlot
if (group_comparison == "") {
p_cov <- CoveragePlot(
object = obj,
region = peaks.top[1], # 这里只取 Top 1 画图，避免图太大
extend.upstream = 5000,
extend.downstream = 5000,
nrow = 1
)
} else {
p_cov <- CoveragePlot(
object = obj,
region = peaks.top[1],
extend.upstream = 5000,
extend.downstream = 5000,
nrow = 1,
group.by = group_comparison,
idents = c(case_group, control_group)
)
}
# 打印到屏幕（Notebook 专用）
print(p_cov)
# 保存 PDF 到本地
ggsave(paste0(output, '/DIFF/02_topPeak/', prefix, '_diffPeak_top.pdf'),
p_cov, width = 24, height = 3)
}
}, error = function(e) {
message(sprintf("Error drawing CoveragePlot for cluster %s, group %s: %s", i, grp, e$message))
})
}
}, error = function(e) {
message(sprintf("Error processing cluster %s: %s", i, e$message))
})
}R
options(repr.plot.width = 10, repr.plot.height = 4)
p_cov
💡 说明
该图示例展示其中一组Top差异Peak在不同组或cluster中的ATAC覆盖信号差异。每个小面板对应一个候选关键Peak附近区域，通过比较上下多组信号峰高和峰形，验证统计差异是否在原始信号层面可见。
- 上/下轨道：分别是两组或多个细胞簇在同一区域的标准化ATAC覆盖信号；峰越高代表该区域可及性越强。
- Normalized signal：纵轴为标准化覆盖强度，可在同一面板内比较两组高低。
- chr position (bp)：横轴为基因组坐标，显示该Peak及其上下游邻域。
- Genes 轨道：显示邻近基因结构与方向，便于判断差异Peak可能影响的功能基因。
- Peak 标记：灰色短条为调用到的开放区间，重点看Top差异Peak所在位置是否与组间信号变化一致。
- 多面板并列：每个面板一个Top peak位点，便于快速筛出“统计显著且图像上也分离明显”的高可信候选区域。
4.3 差异 Peak 全局热图展示
利用 ComplexHeatmap 绘制热图，全景展示特异性开放区域的信号变化：
数据平滑与标准化 提取各组 Top 500 差异 Peak 进行伪批量化（Pseudobulk）合并，并通过 Z-score 标准化增强对比度，避免单细胞层面的数据过于稀疏。
双重展示模式 根据分析策略，自动生成“所有细胞群 Marker 大热图”（区分不同细胞群）或“特定细胞群组间差异热图”（对比实验组与对照组）。
Top Peak 精准标注 动态提取各组排名前 5 的最显著差异 Peak 坐标，利用
anno_mark功能在热图右侧进行精准的连线标注，帮助快速定位变化最剧烈的核心基因组区域。
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
}
if (group_comparison == "") {
all_top_peaks <- c()
peak_cluster_anno <- c()
# 新增：用于收集要标注的 Peak 及其所在分组
mark_peaks_list <- list()
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))
# 提取当前 Cluster 前 5 个最显著的 Peak 用于标注
top5_peaks <- head(sel_peaks, 5)
mark_peaks_list[[as.character(i)]] <- top5_peaks
}
}
if (length(all_top_peaks) > 0) {
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
)
# 计算被选中标注的 Peak 在最终矩阵中的确切行号
target_peaks <- unlist(mark_peaks_list)
target_idx <- match(target_peaks, rownames(mat_scaled))
# 去除可能出现的 NA（虽然理论上不会）
valid_idx <- !is.na(target_idx)
mark_at <- target_idx[valid_idx]
mark_labels <- target_peaks[valid_idx]
right_anno <- rowAnnotation(
Peaks = anno_mark(
at = mark_at, # 具体到每一行的索引
labels = mark_labels, # 对应的真实 Peak 坐标
which = "row",
side = "right",
labels_gp = gpar(fontsize = 8, fontface = "bold", col = "black"),
lines_gp = gpar(col = "black", lty = 1),
link_width = unit(4, "mm"),
padding = unit(0.5, "mm"),
extend = unit(c(1, 1), "mm")
)
)
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()
# 新增：用于收集当前 cluster 内部两组的 Top 5 Peak
mark_peaks_list <- list()
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))
#提取该组的前 5 个最显著 Peak
mark_peaks_list[[grp]] <- head(grp_peaks, 5)
}
}
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) {
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
)
#：计算被选中标注的 Peak 在当前热图矩阵中的行号
target_peaks <- unlist(mark_peaks_list)
target_idx <- match(target_peaks, rownames(mat_scaled))
valid_idx <- !is.na(target_idx)
mark_at <- target_idx[valid_idx]
mark_labels <- target_peaks[valid_idx]
right_anno <- rowAnnotation(
Peaks = anno_mark(
at = mark_at,
labels = mark_labels,
which = "row",
side = "right",
# 1. 稍微缩小字体，让出更多空间 (从 8 改为 6 或 7)
labels_gp = gpar(fontsize = 6, fontface = "bold", col = "black"),
lines_gp = gpar(col = "black", lty = 2, lwd = 0.5),
# 2. 拉长连线：让文字离热图稍微远一点，视觉上更开阔
link_width = unit(8, "mm"),
# 3. 增大文本间距：强制相邻的文本上下散开
padding = unit(1.5, "mm"),
# 4. 扩大延伸范围：允许连线向热图的上下方寻找更多空白区域来放置文本
extend = unit(c(15, 15), "mm")
)
)
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"))
5. 差异 Peak 过滤与邻近基因注释
5.1 差异 Peak 邻近基因注释
在获取各细胞群或实验组间的初始差异 Peak 后，需要进一步对结果进行严格过滤与生物学注释，从而将抽象的表观遗传信号转化为具体的基因调控线索。针对筛选出的高置信度差异 Peaks，我们可以利用 Signac 的 ClosestFeature 函数，将其物理坐标精确映射到参考基因组上距离最近的基因。这一步将非编码区的表观遗传开放变化，直接关联到潜在被调控的靶基因上。注释结果（包含基因名称、物理距离等信息）输出至 03_closestFeature 目录，为后续的 GO/KEGG 功能通路富集分析奠定基础。
R
for (i in cluster) {
tryCatch({
DefaultAssay(obj) <- "ATAC"
data <- da_peaks[da_peaks$cluster == i, ]
# 1. 使用正则表达式检查 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("正在分析：", grp))
tryCatch({
# 2. 差异 Peak 的显著性阈值过滤
if (grp == "cluster_marker") {
# For single cluster, filter by padj and 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))
}
# 如果没有符合条件的显著 Peak，跳过该组
if (nrow(data.sig) == 0) {
message(sprintf("No significant peaks found for %s", prefix))
next
}
print(dim(data.sig))
# 导出过滤后的显著差异 Peak 表格
write.table(data.sig, paste0(output, '/DIFF/01_diffPeak/', prefix, '_diffPeak_significant.xls'),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = '\t')
# 3. 差异 Peak 邻近靶基因注释
peaks <- data.sig$feature
genes <- ClosestFeature(obj, regions = peaks)
# 补充元数据信息
genes$cluster <- as.character(i)
if (grp != "cluster_marker") {
genes$group <- as.character(grp)
}
# 导出就近基因注释表格
write.table(genes, paste0(output, '/DIFF/03_closestFeature/', prefix, '_diffPeak_sig_genes.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(i, " 分析完成"))
}, error = function(e) {
message(sprintf("Error in cluster %s: %s", i, e$message))
})
}R
data_dir <- paste0(output,"/DIFF/03_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)
head(merged_data)| tx_id | gene_name | gene_id | gene_biotype | type | closest_region | query_region | distance | cluster | group | |
|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <int> | <chr> | <chr> | |
| 1 | ENST00000423796 | AC114498.1 | ENSG00000235146 | lncRNA | exon | chr1-594235-594768 | chr1-630357-631124 | 35588 | B cells | 25030508_pbmc_1_arc |
| 2 | ENST00000481276 | SLC35E2B | ENSG00000189339 | protein_coding | exon | chr1-1692449-1692709 | chr1-1692209-1693568 | 0 | B cells | 25030508_pbmc_1_arc |
| 3 | ENST00000650835 | LINC01964 | ENSG00000260840 | lncRNA | exon | chr2-85065895-85067347 | chr2-85089328-85090330 | 21980 | B cells | 25030508_pbmc_1_arc |
| 4 | ENST00000666960 | AL078621.4 | ENSG00000287165 | lncRNA | exon | chr2-113582795-113583152 | chr2-113603235-113604680 | 20082 | B cells | 25030508_pbmc_1_arc |
| 5 | ENST00000440223 | UBE2F | ENSG00000184182 | protein_coding | exon | chr2-237966827-237967132 | chr2-237943150-237943999 | 22827 | B cells | 25030508_pbmc_1_arc |
| 6 | ENST00000453168 | STT3B | ENSG00000163527 | protein_coding | exon | chr3-31532638-31533312 | chr3-31226203-31227780 | 304857 | B cells | 25030508_pbmc_1_arc |
💡 说明
上表展示了差异开放 Peak 映射到基因组上最近基因（Closest Feature）的结果。每一行代表一个差异 Peak 及其对应的潜在靶基因，用于回答“这个异常开放的染色质区域最可能调控哪个基因”。
- tx_id：转录本 ID（如
ENST00000423796），代表与该差异 Peak 物理距离最近的具体转录本。- gene_name：基因符号（如
SLC35E2B），即上述转录本所属的常用基因名称，便于直观了解基因功能。- gene_id：基因的 Ensembl ID（如
ENSG00000235146），用于唯一标识该基因。- gene_biotype：基因的生物学类型，例如
protein_coding（蛋白编码基因）或lncRNA（长链非编码 RNA），帮助判断调控产物的性质。- type：Peak 落入的基因组元件类型，如
exon（外显子）、intron（内含子）、promoter（启动子）等。- closest_region：基因组上参考注释元件的具体坐标区间。
- query_region：我们输入的、发生显著差异开放的 Peak 的坐标区间（如
chr1-630357-631124）。- distance：该差异 Peak 距离最近基因元件的物理距离（单位为碱基对 bp）。如果距离为 0，说明该 Peak 直接与靶基因存在重叠（Overlap）。
- cluster：细胞群名称，标明该 Peak 是在哪个特定的细胞群（如
B cells）中鉴定出来的。- group：组别信息，标明该 Peak 是在哪个比较组（如
25030508_pbmc_1_arc）中呈现特异性高开放状态。
5.2 邻近基因功能注释
在完成了从坐标到基因的映射后，可进一步将这些零散的靶基因放入宏观的生物学网络中进行考察。通过执行 GO 与 KEGG 富集分析，理解这些发生染色质开放状态改变的区域，最终究竟驱动了细胞怎样的表型变化。
该部分分析结果文件说明： 所有的富集分析结果均保存在 03_closestFeature 输出目录中。针对每一个细胞群或对比组（如 prefix），程序会自动生成以下结果文件：
- 表格数据：包含完整的富集统计信息及对应靶基因列表（如
[prefix]_GOenrich.xls和[prefix]_KEGGenrich.xls）。 - 可视化图表：分别展示 Top 20 显著富集通路的柱状图（
*_GO_bar.pdf、*_KEGG_bar.pdf）与气泡图（*_GO_dot.pdf、*_KEGG_dot.pdf），便于直接用于后续的报告和展示。
R
if (species == "human") {
kegg_code <- "hsa" # 新建专门表示 KEGG 缩写的变量，不覆盖原 species
db <- "org.Hs.eg.db"
spe <- "human" # 如果你后续非要用 spe，可以保留，但其实直接用原变量 species 就行
} else if (species == "mouse") {
kegg_code <- "mmu"
db <- "org.Mm.eg.db"
spe <- "mouse"
} else {
stop(paste0("不支持该物种: '", species,
"'\n请确保已安装对应物种的基因注释数据库。\n",
"可用选项: 'human' 或 'mouse'"))
}R
files <- list.files(paste0(output,"/DIFF/03_closestFeature"),pattern="diffPeak_sig_genes.xls")
for (f in files){
prefix <- gsub('_diffPeak_sig_genes.xls','',f)
f <- paste0(output,"/DIFF/03_closestFeature/",f)
print(f)
print(prefix)
clus_1<-read.table(f,header=T,sep='\t')
clus<-clus_1$gene_name
eg <- NULL
tryCatch({
eg <- bitr(clus, fromType="SYMBOL", toType=c("ENSEMBL","ENTREZID"), OrgDb=db)
print(head(eg))
}, error = function(e) {
message(paste("处理文件", f, "时出错:", e$message))
})
if (is.null(eg) || nrow(eg) == 0) {
message(paste0("跳过富集分析（无有效映射基因）: ", prefix))
next
}
tryCatch({
go_enrich(eg,db,paste0(output,"/03_closestFeature/clusterProfiler/GO"),prefix)
}, error = function(e) {
message(paste0("GO富集分析错误: ", prefix, " - ", e$message))
})
tryCatch({
kegg_enrich(eg, species, paste0(output,"/03_closestFeature/clusterProfiler/KEGG"), prefix)
}, error = function(e) {
message(paste0("KEGG富集分析错误: ", prefix, " - ", e$message))
})
}R
# 合并 GO 富集分析结果
go_dir <- paste0(output,"/03_closestFeature/clusterProfiler/GO")
merged_data_go <- data.frame()
if(dir.exists(go_dir)) {
xls_files <- list.files(path = go_dir,
pattern = "\\.xls$",
full.names = TRUE)
for (file in xls_files) {
df <- read.table(file, sep="\t", header=TRUE, fill=TRUE, quote="", check.names=FALSE)
merged_data_go <- bind_rows(merged_data_go, df)
}
if (nrow(merged_data_go) > 0) {
write.table(merged_data_go,
file = file.path(go_dir, "all_GOenrich.xls"),
row.names = FALSE,
sep="\t",
quote=FALSE)
}
}R
head(merged_data_go)| cluster | ONTOLOGY | ID | Description | GeneRatio | BgRatio | RichFactor | FoldEnrichment | zScore | pvalue | p.adjust | qvalue | geneID | Count | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <int> | |
| 1 | c_Bcells | BP | GO:0030098 | lymphocyte differentiation | 101/1678 | 432/18805 | 0.2337963 | 2.620107 | 10.662915 | 7.040828e-20 | 4.066782e-16 | 3.474463e-16 | PRKCZ/PLA2G2D/ARID1A/PIK3R3/SEMA4A/NTRK1/CD1D/SLAMF6/LY9/PTPRC/HLX/SP3/ITGA4/DOCK10/INPP5D/HDAC4/PLCL2/TGFBR2/CMTM7/CCR9/TLR9/FOXP1/RHOH/AP3B1/FNIP1/IL4/CD74/IRF4/CD83/SOX4/HLA-DRB1/HLA-DOA/RUNX2/PRDM1/FOXO3/TRAF3IP2/VNN1/MYB/ARID1B/CCR6/LFNG/CARD11/ACTB/HDAC9/IKZF1/CDK6/CHD7/TPD52/SMARCA2/IFNA2/PAX5/SHB/SYK/TNFSF8/FUT7/IL2RA/VSIR/ZMIZ1/HHEX/BLNK/RAG2/SPI1/PTPRJ/MS4A1/MEN1/SYVN1/POU2AF1/RNF41/STAT6/USP44/DTX1/ZFP36L1/JAG2/DLL4/CD19/PLCG2/IRF8/ZFPM1/IKZF3/RARA/CCR7/SMARCE1/CD79B/C17orf99/PTPN2/TCF3/AP3D1/CLEC4G/SMARCA4/LYL1/BRD4/IL12RB1/NFKBID/CD79A/CLPTM1/LILRB2/RUNX1/SMARCB1/PATZ1/NFAM1/SASH3 | 101 |
| 2 | c_Bcells | BP | GO:0042113 | B cell activation | 73/1678 | 287/18805 | 0.2543554 | 2.850509 | 9.888159 | 8.608820e-17 | 2.486227e-13 | 2.124113e-13 | LAPTM5/PIK3R3/VAV3/NTRK1/FCRL1/PTPRC/MSH6/SP3/ITGA4/SLC39A10/DOCK10/INPP5D/HDAC4/PLCL2/CMTM7/TLR9/PRKCD/CD38/BANK1/GAPT/MEF2C/FNIP1/IL4/MZB1/CD74/CDKN1A/TNFRSF21/AKIRIN2/TRAF3IP2/CCR6/LFNG/CARD11/HDAC9/LYN/TPD52/IFNA2/PAX5/SHB/SYK/HHEX/BLNK/CD81/SWAP70/RAG2/SPI1/PTPRJ/MS4A1/SYVN1/TBC1D10C/POU2AF1/STAT6/SLC15A4/ZFP36L1/PRKCB/CD19/NOD2/PLCG2/IRF8/IKZF3/CD79B/C17orf99/PTPN2/TCF3/LYL1/CD22/TYROBP/CD79A/ERCC1/SHLD1/NFATC2/TNFRSF13C/NFAM1/SASH3 | 73 |
| 3 | c_Bcells | BP | GO:0050853 | B cell receptor signaling pathway | 33/1678 | 78/18805 | 0.4230769 | 4.741336 | 10.363831 | 3.495526e-15 | 6.730053e-12 | 5.749834e-12 | VAV3/PTPRC/FCMR/IGKC/SLC39A10/PLCL2/FOXP1/KLHL6/CD38/BANK1/MEF2C/SH2B2/BLK/LYN/CD72/SYK/BLNK/CD81/MS4A1/PRKCH/IGHA2/IGHA1/IGHD/PRKCB/CD19/PLCG2/CD79B/CD22/CD79A/NFATC2/MAPK1/IGLC7/NFAM1 | 33 |
| 4 | c_Bcells | BP | GO:0051251 | positive regulation of lymphocyte activation | 74/1678 | 337/18805 | 0.2195846 | 2.460839 | 8.470086 | 2.109515e-13 | 2.769138e-10 | 2.365818e-10 | PRKCZ/ARID1A/VAV3/CD1D/PTPRC/HLX/NCK2/SLC39A10/IGFBP2/INPP5D/TGFBR2/TLR9/HES1/CD38/RHOH/PPP3CA/AP3B1/MEF2C/IL4/CD74/CD83/SOX4/HLA-DRB1/HLA-DQA1/HLA-DOB/HLA-DMB/HLA-DOA/CDKN1A/AKIRIN2/CD24/FOXO3/FYN/VNN1/ARID1B/CARD11/ACTB/LYN/DOCK8/SMARCA2/SHB/SYK/IL2RA/MAP3K8/VSIR/ZMIZ1/IGF2/CD81/SPI1/CCDC88B/STAT6/AKT1/CSK/NOD2/RARA/CCR7/SMARCE1/TCF3/AP3D1/TMIGD2/CD209/SMARCA4/BRD4/FCHO1/IL12RB1/NFKBID/TYROBP/LILRB2/SHLD1/NFATC2/RUNX1/SMARCB1/TNFRSF13C/EFNB1/SASH3 | 74 |
| 5 | c_Bcells | BP | GO:0030217 | T cell differentiation | 72/1678 | 324/18805 | 0.2222222 | 2.490399 | 8.470180 | 2.397107e-13 | 2.769138e-10 | 2.365818e-10 | PRKCZ/PLA2G2D/ARID1A/PIK3R3/SEMA4A/CD1D/SLAMF6/LY9/PTPRC/HLX/SP3/TGFBR2/CCR9/FOXP1/RHOH/AP3B1/IL4/CD74/IRF4/CD83/SOX4/HLA-DRB1/HLA-DOA/RUNX2/PRDM1/FOXO3/TRAF3IP2/VNN1/MYB/ARID1B/CCR6/LFNG/CARD11/ACTB/CDK6/CHD7/SMARCA2/IFNA2/SHB/SYK/TNFSF8/FUT7/IL2RA/VSIR/ZMIZ1/RAG2/SPI1/MEN1/STAT6/USP44/DTX1/ZFP36L1/JAG2/DLL4/ZFPM1/IKZF3/RARA/CCR7/SMARCE1/PTPN2/AP3D1/CLEC4G/SMARCA4/BRD4/IL12RB1/NFKBID/CLPTM1/LILRB2/RUNX1/SMARCB1/PATZ1/SASH3 | 72 |
| 6 | c_Bcells | BP | GO:1903037 | regulation of leukocyte cell-cell adhesion | 81/1678 | 395/18805 | 0.2050633 | 2.298102 | 8.161330 | 7.459736e-13 | 7.181239e-10 | 6.135306e-10 | PRKCZ/PLA2G2D/ARID1A/PTAFR/LAPTM5/CD1D/PTPRC/HLX/NCK2/ITGA4/IGFBP2/TGFBR2/CBLB/HES1/RHOH/PPP3CA/AP3B1/IL4/CD74/DUSP22/CD83/SOX4/RIPOR2/HLA-DRB1/HLA-DQA1/HLA-DOB/HLA-DMB/HLA-DOA/TNFRSF21/CD24/FOXO3/FYN/ARG1/VNN1/ARID1B/MAD1L1/CARD11/ACTB/LYN/DOCK8/SMARCA2/IFNA2/SHB/GCNT1/SYK/FUT7/IL2RA/MAP3K8/VSIR/ZMIZ1/IGF2/CD81/RAG2/CCDC88B/LRRC32/ETS1/WNK1/DTX1/AKT1/CSK/NOD2/RARA/CCR7/SMARCE1/PTPN2/AP3D1/TMIGD2/CLEC4G/CD209/SMARCA4/BRD4/FCHO1/IL12RB1/NFKBID/LILRB2/RUNX1/SMARCB1/ADORA2A/TNFRSF13C/EFNB1/SASH3 | 81 |
R
kegg_dir <- paste0(output,"/03_closestFeature/clusterProfiler/KEGG")
merged_data_kegg <- data.frame()
if(dir.exists(kegg_dir)) {
xls_files <- list.files(path = kegg_dir,
pattern = "\\.xls$",
full.names = TRUE)
for (file in xls_files) {
df <- read.table(file, sep="\t", header=TRUE, fill=TRUE, quote="", check.names=FALSE)
merged_data_kegg <- bind_rows(merged_data_kegg, df)
}
if (nrow(merged_data_kegg) > 0) {
write.table(merged_data_kegg,
file = file.path(kegg_dir, "all_KEGGenrich.xls"),
row.names = FALSE,
sep="\t",
quote=FALSE)
}
}R
head(merged_data_kegg)| cluster | category | subcategory | ID | Description | GeneRatio | BgRatio | pvalue | p.adjust | qvalue | geneName | Count | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <chr> | <int> | |
| 1 | c_Erythroblast | Cellular Processes | Transport and catabolism | hsa04144 | Endocytosis | 127/3122 | 250/9382 | 5.188384e-09 | 1.821123e-06 | 1.070445e-06 | ACAP3/LDLRAP1/EPS15/DNAJC6/SH3GLB1/DNM3/ASAP2/RAB10/EHD3/RAB11FIP5/PSD4/ACTR3/CXCR4/STAM2/WIPF1/CXCR1/AGAP1/CAV3/IQSEC1/RAB5A/TGFBR2/PDCD6IP/RHOA/ARF4/CHMP2B/CBLB/SNX4/RAB7A/PRKCI/PLD1/ACAP2/TFRC/FGFR3/GRK4/ARAP2/PDGFRA/FGFR4/RUFY1/HLA-G/HLA-E/HSPA1A/HSPA1B/SNX3/IGF2R/CYTH3/WIPF3/SMURF1/CAV2/CAPZA2/AGAP3/VPS37A/RAB11FIP1/ARFGEF1/CHMP4C/ASAP1/CHMP5/PIP5K1B/DNM1/SH3GLB2/IL2RA/STAM/KIF5B/PARD3/AGAP5/ZFYVE27/GRK5/FGFR2/AP2A2/VPS37C/EHD1/ARRB1/HSPA8/VPS26B/IQSEC3/CAPZA3/RNF41/KIF5A/MDM2/EEA1/WASHC4/GIT2/ARPC3/VPS37B/GRK1/SNX6/ARF6/NEDD4/SNX1/RAB11A/SMAD3/SH3GL3/IGF1R/RAB11FIP3/ARRB2/PLD2/EPN2/AP2B1/WIPF2/EPN3/CLTC/SMURF2/CYTH1/CHMP6/RAB31/SMAD2/NEDD4L/VPS4B/SH3GL1/DNM2/LDLR/EPS15L1/AP2S1/EHD2/CYTH2/AP2A1/EPN1/SNX5/CHMP4B/ITCH/SRC/ARFGEF2/CLTCL1/GRK3/IL2RB/ARFGAP3/SH3KBP1/IQSEC2 | 127 |
| 2 | c_Erythroblast | Environmental Information Processing | Signal transduction | hsa04310 | Wnt signaling pathway | 94/3122 | 178/9382 | 4.748674e-08 | 8.333922e-06 | 4.898632e-06 | DVL1/ROR1/PRKACB/VANGL1/VANGL2/LGR6/WNT9A/ROCK2/PPP3R1/FZD7/FZD5/WNT6/WNT7A/CELSR3/RHOA/GSK3B/RUVBL1/RYK/TBL1XR1/SENP2/MAPK10/PPP3CA/LEF1/SFRP2/CTNND2/APC/MCC/TCF7/CSNK1A1/CAMK2A/FBXW11/CCND3/ANKRD6/MAP3K7/RAC1/SFRP4/FZD1/CUL1/FZD3/FZD6/MYC/PRKACG/ROR2/INVS/DKK1/PPP3CB/FRAT1/FRAT2/SFRP5/BTRC/TCF7L2/CTBP2/CSNK2A3/LGR4/FOSL1/LRP5/WNT11/WNT5B/CCND2/WIF1/LGR5/CSNK1A1L/DAAM1/CCDC88C/SMAD3/TLE3/AXIN1/CREBBP/PRKCB/SIAH1/NFATC3/TLE7/NLK/FZD2/WNT9B/RNF43/PRKCA/RAC3/APCDD1/NFATC1/APC2/TLE2/CSNK2A1/PLCB1/PLCB4/NFATC2/ZNRF3/CSNK1E/RBX1/CELSR1/TBL1X/PRICKLE3/GPC4/TBL1Y | 94 |
| 3 | c_Erythroblast | Cellular Processes | Cellular community - eukaryotes | hsa04510 | Focal adhesion | 103/3122 | 203/9382 | 1.575654e-07 | 1.843515e-05 | 1.083607e-05 | PIK3R3/VAV3/RAP1A/TNR/LAMC2/PPP1R12B/LAMB3/CAPN2/ROCK2/SOS1/ITGA4/FN1/COL4A3/CAV3/RAF1/ITGA9/LAMB2/RHOA/FLNB/GSK3B/MYLK/ITGB5/COL6A5/COL6A6/PIK3CB/PIK3CA/PDGFRA/MAPK10/IBSP/PDGFC/ITGA2/PIK3R1/THBS4/PDGFRB/MYLK4/TNXB/CCND3/VEGFA/FYN/LAMA2/PDGFA/RAC1/ITGB8/HGF/RELN/CAV2/BRAF/PTK2/PIP5K1B/SHC3/LAMC3/VAV2/ITGA8/ITGB1/VCL/PTEN/BAD/PAK1/BIRC3/BIRC2/CCND2/VWF/EMP1/ITGB7/RAP1B/PPP1CC/PXN/FLT1/COL4A2/ARHGAP5/SOS2/ACTN1/PAK6/SHC4/TLN2/MAP2K1/IGF1R/PDPK1/EMP2/PRKCB/MYLK3/BCAR1/CRK/ERBB2/ITGA2B/ITGB3/PRKCA/RAC3/MYL12A/LAMA3/BCL2/VAV1/COMP/ACTN4/PAK4/MYLK2/SRC/LAMA5/COL9A3/MAPK1/PARVB/ELK1/COL4A5 | 103 |
| 4 | c_Erythroblast | Environmental Information Processing | Signal transduction | hsa04151 | PI3K-Akt signaling pathway | 165/3122 | 362/9382 | 4.772663e-07 | 3.865820e-05 | 2.272307e-05 | GNB1/CASP9/EPHA2/PIK3R3/JAK1/GNG12/PKN2/CSF1/IL6R/EFNA4/NTRK1/FASLG/TNR/LAMC2/LAMB3/PPP2R5A/SOS1/BCL2L11/ATF2/ITGA4/CREB1/FN1/IRS1/COL4A3/EIF4E2/RAF1/ITGA9/LAMB2/GSK3B/ITGB5/COL6A5/COL6A6/PPP2R3A/PIK3CB/PIK3CA/GNB4/FGFR3/PPP2R2C/PDGFRA/KIT/EREG/AREG/FGF5/IBSP/EIF4E/NFKB1/FGF2/PDGFC/GHR/ITGA2/PIK3R1/THBS4/IL4/CSF1R/PDGFRB/FGF18/FGFR4/TNXB/ATF6B/CCND3/VEGFA/FOXO3/LAMA2/SGK1/MYB/PDGFA/RAC1/ITGB8/CREB5/YWHAG/MAGI2/HGF/CDK6/RELN/PIK3CG/CREB3L2/NOS3/RHEB/ANGPT2/PPP2R2A/EIF4EBP1/SGK3/ANGPT1/MYC/PTK2/JAK2/IFNB1/TEK/SYK/LAMC3/TSC1/RXRA/IL2RA/ITGA8/ITGB1/RET/DDIT4/PTEN/PIK3AP1/FGF8/FGFR2/CREB3L1/CHRM1/BAD/FGF19/CCND2/NTF3/VWF/NR4A1/ITGB7/MDM2/HSP90B1/FGF9/FLT1/COL4A2/PCK2/SOS2/GNG2/PPP2R5E/PPP2R5C/HSP90AA1/FGF7/MAP2K1/IGF1R/PDPK1/IL4R/CD19/RBL2/YWHAE/ERBB2/BRCA1/ITGA2B/ITGB3/GNGT2/PRKCA/RPTOR/LAMA3/PHLPP1/BCL2/STK11/MAP2K2/CREB3L3/INSR/EPOR/PKN1/JAK3/COMP/GNG8/FGF21/GYS1/FLT3LG/ANGPT4/BCL2L1/SGK2/PCK1/LAMA5/COL9A3/IFNAR2/IFNAR1/MAPK1/OSM/IL2RB/IL3RA/FGF16/COL4A5 | 165 |
| 5 | c_Erythroblast | NA | NA | hsa04517 | IgSF CAM signaling | 140/3122 | 299/9382 | 5.506865e-07 | 3.865820e-05 | 2.272307e-05 | PTPRF/PIK3R3/SSX2IP/VAV3/RAP1A/CD2/CADM3/CD84/LY9/CD244/F11R/SH2D1B/NFASC/CNTN2/SRGAP2/PLXNA2/SPTBN1/DOK1/NCK2/ACTR3/ITGA4/CD28/NRP2/INPP5D/RHOA/ROBO1/GSK3B/MYLK/PLXNA1/NCK1/PIK3CB/PRKCI/NLGN1/PIK3CA/IL1RAP/MAPK10/FGB/RAPGEF2/TRIO/PIK3R1/SLIT3/UNC5A/MYLK4/TUBB2A/MAPK14/FYN/EPB41L2/EZR/AFDN/RAC1/MAGI2/PLXNA4/CNTNAP2/DOK2/LYN/MTSS1/PTK2/PTPRD/TJP2/SPTAN1/ABL1/VAV2/TUBB8/PRKCQ/ITGB1/PARD3/ANK3/UNC5B/VCL/SORBS1/SLIT1/SPTBN2/LRFN4/PPFIA1/PAK1/RDX/NECTIN1/KIRREL3/ERC1/TUBA1B/GRIP1/RAP1B/ARPC3/PXN/FOXO1/LMO7/ARF6/ACTN1/PAK6/MAP2K1/NTRK3/CASKIN1/PDPK1/MYH11/PRKCB/ITGAL/ITGAM/MYLK3/CDH1/MTSS2/BCAR1/PLCG2/TUBB3/MYH10/NTN1/PMP22/MPP3/ITGB3/ICAM2/PRKCA/BAIAP2/MYL12A/TUBB6/CABLES1/CD226/MBP/MAP2K2/VAV1/CLEC4M/MAG/ACTN4/PAK4/SPTBN4/CADM4/PVR/NECTIN2/PPFIA3/MYH14/MYLK2/SRC/TUBB1/CXADR/ITGB2/MAPK1/MYH9/MAPK11/CASK/MSN/SH2D1A/PLXNA3 | 140 |
| 6 | c_Erythroblast | Environmental Information Processing | Signal transduction | hsa04015 | Rap1 signaling pathway | 104/3122 | 212/9382 | 1.118055e-06 | 5.800961e-05 | 3.409772e-05 | EPHA2/RAP1GAP/PIK3R3/VAV3/CSF1/RAP1A/MAGI3/EFNA4/SIPA1L2/ADCY3/RALB/FARP2/RAF1/RHOA/GNAI2/ADCY5/PIK3CB/P2RY1/PRKCI/PIK3CA/FGFR3/PDGFRA/KIT/FGF5/FGF2/PDGFC/RAPGEF2/FYB1/PIK3R1/CSF1R/PDGFRB/FGF18/FGFR4/RGS14/MAPK14/VEGFA/AFDN/PDGFA/RAC1/RAPGEF5/RALA/ADCY1/MAGI2/GNAI1/HGF/BRAF/ANGPT2/ANGPT1/TEK/GNAQ/RALGDS/VAV2/GRIN1/CALML5/APBB1IP/ITGB1/PARD3/PLCE1/FGF8/FGFR2/CTNND1/FGF19/DRD2/RAP1B/FGF9/FLT1/PRKD1/EVL/FGF7/TLN2/MAP2K1/IGF1R/PRKCB/LAT/ITGAL/ITGAM/ADCY7/GNAO1/CDH1/BCAR1/CRK/ADORA2B/MAP2K3/ITGA2B/ITGB3/SKAP1/PRKCA/RAC3/MAP2K2/VAV1/INSR/SIPA1L3/PRKD2/FGF21/ANGPT4/PLCB1/PLCB4/ID1/SRC/TIAM1/ITGB2/MAPK1/MAPK11/FGF16 | 104 |