支持中心

简体中文

English

简体中文

English

Appearance

ATAC + RNA 多组学：Footprint分析 ​

作者: SeekGene
时长: 14 分钟
字数: 3.4k 字
更新: 2026-07-07
阅读: 0 次

1. 教程简介 ​

本教程基于 Seurat 和 Signac 等单细胞多组学分析工具，专门用于对 SeekArc 单细胞多组学数据（scRNA-seq + scATAC-seq） 进行转录因子足迹（Footprint）分析。

在单细胞 ATAC-seq 分析中，虽然 Motif 富集分析可以告诉我们哪些转录因子的结合位点在特定细胞群或差异开放区域中富集，但它无法证明转录因子是否实际结合在那里。

转录因子足迹（Footprint）分析通过在染色质切割轨迹层面客观评估 TF 潜在的物理结合行为。其基本原理是：当 TF 结合在 DNA 上时，会阻挡 Tn5 酶对该区域的切割，从而在富集峰的中心形成一个信号相对较低的保护区（即“足迹”），而两侧则出现信号高峰。

通过对比不同细胞群或实验组别中该保护区深度的变化，可以进一步验证候选转录因子在目标状态下是否发生了真实的结合与脱落事件。

R
suppressPackageStartupMessages(suppressWarnings({
    library(presto)
    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(foreach)
    library(doParallel)
    library(base64enc) 
}))

2. 输入文件准备 ​

2.1 输入文件要求 ​

开始分析前，需要准备以下几类核心输入文件：

  • input.rds：预处理完成的 Seurat 对象文件。该对象必须包含 RNA 和 ATAC 两个 Assay。
  • meta.tsv：细胞元数据文件（制表符分隔）。必须包含 barcode 列，以及用于定义细胞类型（如 Celltype）和样本分组（如 Sample）的对应列。
  • genome.fa：对应物种的参考基因组 FASTA 文件。用于提取 DNA 序列。
  • ATAC 片段文件及索引：构建 ATAC Assay 所依赖的 fragments.tsv.gz 及其 .tbi 索引文件。这是进行 Footprint 分析最核心的文件。
R
# --- 输入参数配置 ---此处参数按照自己需求填写参数

## fpath：fragment文件所在的目录
fpath = "/path/to/PBMC_demo"

## rds：输入的 Seurat 对象 RDS 文件路径，包含已预处理的单细胞多组学数据
rds = "/path/to/input.rds"

## meta：元数据文件路径（TSV格式），包含样本信息、细胞类型注释等元数据
meta = "/path/to/meta.tsv"

## species：物种信息
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"

## motif_names：需要进行 Footprint 分析的 Motif 名称列表，多个用逗号分隔。
tf_names = "GATA2,CEBPA,EBF1"
R
# 参数处理，无需改动
output <- "."
footprint_dir <- paste0(output, "/footprint")
dir.create(footprint_dir, recursive = TRUE, showWarnings = FALSE)
celltypes <- strsplit(celltypes,",")[[1]]
downsample_num <- as.numeric(downsample_num)
tf_list <- strsplit(tf_names, ",")[[1]]

3. 数据加载与预处理 ​

3.1 加载 Seurat 对象与元数据 ​

R
#读取数据，数据预处理
obj <- readRDS(rds)
if (meta != ""){
  meta <- read.table(meta,header=T,sep='	',check.names=F)
  rownames(meta) <- meta$barcode
  obj <- AddMetaData(obj, meta)
}
obj <- subset(obj, subset = !!sym(clusters_col) %in% celltypes)

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

R
#对细胞进行下采样
if (downsample == "TRUE"){
  cells_to_keep <- c()
  for (ctype in celltypes) {
    ctype_cells <- rownames(obj@meta.data[obj@meta.data[[clusters_col]] == ctype, ])
    if (length(ctype_cells) > downsample_num) {
      ctype_cells <- sample(ctype_cells, downsample_num)
    }
    cells_to_keep <- c(cells_to_keep, ctype_cells)
  }
  obj <- subset(obj, cells = cells_to_keep)
}

3.3 修改 fragment 文件路径 ​

在 Seurat 对象的 ATAC 模态中，Fragment 文件并不会把海量序列直接存进内存，而是仅仅保存了文件在磁盘上的路径（Paths）。因此，当数据在不同的服务器、容器之间传来传去时，原本保存在对象中的旧路径经常会因为目录结构改变或权限问题而失效，导致后续无法读取。这里的代码会遍历对象中的所有 Fragment 记录，提取出文件名，并将其与我们指定的最新有效目录（fpath）重新拼接，实现路径的实时更新。

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)
}

3.4 准备参考基因组 ​

Footprint 分析的一个核心难点是 Tn5 酶的切割偏好性校正（Tn5 Sequence Bias Correction）。Tn5 酶在切割 DNA 时并不是完全随机的，它对特定的六聚体（hexamer）序列有固有的切割偏好。如果不去除这种由 DNA 序列本身引起的偏好误差，我们观察到的信号低谷（“足迹”）可能仅仅是因为该区域的序列恰好不容易被 Tn5 切割，而并非真的有转录因子在此结合保护。

为了精确计算和校正这种序列偏好性，底层算法必须获取每一个 Peak 区域真实的 DNA 碱基序列，因此我们需要加载并对齐参考基因组。

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")
# 获取染色体长度
genome_seqlens <- setNames(as.numeric(width(genome)), names(genome))

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)

💡 代码处理细节：

Signac 官网教程要求安装物种对应的 BSgenome R 包，为了避免安装包失败或版本不一致的风险。本教程改用 Biostrings::readDNAStringSet() 直接读取自定义的 FASTA 文件。

特别强调：这里传入的 ref_genome（FASTA 文件）必须与你最初进行上游序列比对（Alignment，例如使用 seekarctools时）所使用的参考基因组 FASTA 文件完全一致，否则会导致坐标错位，使 Footprint 结果完全错误。

此外，代码还进行了严格的染色体名称清洗（去除 FASTA header 中冗余的描述符）并实现了与 ATAC 对象内部染色体的强制对齐。这种处理主动剔除了未在 ATAC 矩阵中出现的杂散染色体（Scaffolds/Contigs），大幅降低了内存消耗并彻底杜绝了后续由于染色体不匹配引发的越界报错。

4. Footprint 分析与可视化 ​

由于计算 Footprint 涉及到大量的底层序列遍历，过程较为耗时，建议开启并行计算。

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 全局变量大小限制
if (requireNamespace("memuse", quietly = TRUE)) {
    total_mem <- as.numeric(memuse::Sys.meminfo()$totalram)
    options(future.globals.maxSize = total_mem * 0.7)
} else {
    options(future.globals.maxSize = 50 * 1024^3)
}
output
Footprint 并行计算已启动，总核心数=32，使用 22 个工作线程

4.1 提取 Footprint 矩阵 ​

我们使用 Footprint() 函数提取特定 Motif 附近的 Tn5 插入信号矩阵。该过程会定位包含 Motif 的 Peak 区域，并在其上下游一定范围内（默认前后各 250 bp）精细统计切割事件的分布频率。

能回答的生物学问题 通过统计 Tn5 切割轨迹，我们可以直接且客观地观察到：转录因子在其 Motif 结合位点是否发生了真实的物理结合，并成功阻挡了 Tn5 酶的切割从而形成了一个低信号的“足迹保护区”。

R
# 直接提取匹配到的 Motif TF 名称，以便出图时显示 TF 名称
motif_obj <- Motifs(obj[["ATAC"]])
motif_names_all <- as.character(motif_obj@motif.names)
motif_ids_all <- names(motif_obj@motif.names)

motifs_use <- c()
for (tf in tf_list) {
    # 支持大小写不敏感匹配
    idx <- grep(paste0("^", tf, "$"), motif_names_all, ignore.case = TRUE)
    if (length(idx) == 0) {
        idx <- grep(tf, motif_names_all, ignore.case = TRUE)
    }
    if (length(idx) > 0) {
        # 直接使用匹配到的 TF 名称
        motifs_use <- c(motifs_use, motif_names_all[idx[1]])
        message(paste0("Found Motif name for ", tf, ": ", motif_names_all[idx[1]]))
    } else {
        warning(paste0("Warning: Motif for TF ", tf, " not found in the Seurat object."))
    }
}
motifs_use <- unique(motifs_use)
if (length(motifs_use) == 0) {
    stop("Error: None of the provided TF names were found in the Motif database.")
}

# 对齐染色体信息
atac_ranges <- granges(obj[["ATAC"]])
valid_chr <- intersect(as.character(seqlevels(atac_ranges)), names(genome_seqlens))
atac_ranges <- keepSeqlevels(atac_ranges, valid_chr, pruning.mode = "coarse")
seqlengths(atac_ranges)[valid_chr] <- genome_seqlens[valid_chr]     

# 过滤越界 Peak
atac_ranges <- trim(atac_ranges)
keep_idx <- start(atac_ranges) >= 1 & end(atac_ranges) <= seqlengths(atac_ranges)[as.character(seqnames(atac_ranges))]
keep_idx[is.na(keep_idx)] <- FALSE
obj[["ATAC"]]@ranges <- atac_ranges[keep_idx]
DefaultAssay(obj) <- "ATAC"

# 尝试直接计算，如果因为边界问题失败，则使用安全的降级策略
fp_obj <- NULL
use_fallback <- FALSE
tryCatch({
    fp_obj <- Footprint(
        object = obj,
        motif.name = motifs_use,
        genome = genome,
        assay = "ATAC"
    )
}, error = function(e1) {
    message("直接 Footprint 计算失败，尝试使用安全边界过滤... 错误原因: ", e1$message)
    up <- 250
    down <- 250
    motif_obj <- Motifs(obj[["ATAC"]])
    region_list <- list()
    key_vec <- c()
    for (m in motifs_use) {
        # 提取该 Motif 对应的所有结合位点
        regs <- Signac:::GetFootprintRegions(motif.obj = motif_obj, motif.name = m)
        # 过滤掉不存在于参考基因组中的染色体
        regs <- keepSeqlevels(regs, intersect(seqlevels(regs), names(genome_seqlens)), pruning.mode = "coarse")
        # 过滤掉越过染色体边界的位点
        chr_len <- genome_seqlens[match(as.character(seqnames(regs)), names(genome_seqlens))]
        ok <- !is.na(chr_len) & start(regs) > up & end(regs) <= (chr_len - down)
        regs <- regs[ok]
        if (length(regs) > 0) {
            region_list[[length(region_list) + 1]] <- regs
            key_vec <- c(key_vec, m)
        }
    }
    if (length(region_list) == 0) {
        stop("Fallback failed: No valid motif regions remain after boundary filtering.")
    }
    fp_obj <<- Footprint(
        object = obj,
        regions = region_list,
        key = key_vec,
        genome = genome,
        assay = "ATAC",
        upstream = up,
        downstream = down,
        compute.expected = FALSE # 如果报错，可能需要关闭期望值计算
    )
    # 更新实际使用的 motif 列表
    motifs_use <<- key_vec
    use_fallback <<- TRUE
})

if (is.null(fp_obj)) {
    stop("Footprint 计算最终失败。")
}

💡 进阶说明：为什么本教程的代码比官方单行 Footprint() 调用复杂得多？

理想情况下，直接运行 Footprint(obj, motif.name = "GATA2") 即可完成计算。但在真实的单细胞分析场景中，简单调用极易因底层数据格式或边界细节引发报错崩溃。为保证分析流程的极高健壮性（Robustness），我们在脚本中内置了以下几层自动容错机制：

  1. ** TF 名称智能匹配**：用户通常习惯输入易读的 TF 简称（如 "GATA2"），而底层数据库往往有特定的复杂命名（如 "GATA2::CEBPA" 等）。代码内置了大小写不敏感与模糊匹配逻辑，自动帮你桥接易读名称与底层 Motif ID。
  2. ** 染色体信息严格对齐**：ATAC 的 Peak 数据常夹杂未组装的散碎染色体片段（Scaffolds/Contigs），若参考基因组未包含这些片段，提取序列时会引发致命错误。代码会自动拦截并清洗这些不匹配的片段。
  3. ** 越界峰（Out-of-bounds Peaks）安全裁剪**：Footprint 计算需向 Motif 侧翼延伸 250 bp。若某结合位点紧贴染色体物理边缘，强行延伸将越界导致函数崩溃。代码会提前精准识别并剔除这些高危边缘位点。
  4. ** 异常降级计算策略（Fallback）**：Signac 原生的 Footprint() 极度脆弱，诸如 Tn5 背景期望偏倚计算失败等微小异常都会导致整个流程中断。我们通过 tryCatch 构建了安全网：当常规计算失败时，立即无缝切换至“手动提取+严格边界过滤+关闭期望值计算”的备用方案，最大程度确保你能顺利拿到最终的分析图表。

4.2 可视化 Footprint 结果 ​

获取到 Footprint 矩阵后，我们可以使用 PlotFootprint() 函数进行可视化。我们可以比较不同细胞类型，或者比较同一细胞类型下不同分组的足迹深度差异。

R
if (use_fallback) {
    p_fp_cluster <- PlotFootprint(
        object = fp_obj,
        features = motifs_use,
        group.by = clusters_col,
        assay = "ATAC",
        show.expected = FALSE,
        normalization = "subtract"
    )
} else {
    p_fp_cluster <- PlotFootprint(
        object = fp_obj,
        features = motifs_use,
        group.by = clusters_col,
        assay = "ATAC"
    )
}
R
options(repr.plot.width = 10, repr.plot.height = 10)

p_fp_cluster + patchwork::plot_layout(ncol = 1)
output
Warning message:
“Removed 4545 rows containing missing values or values outside the scale range
(\`geom_label_repel()\`).”
Warning message:
“Removed 4599 rows containing missing values or values outside the scale range
(\`geom_label_repel()\`).”
Warning message:
“Removed 4599 rows containing missing values or values outside the scale range
(\`geom_label_repel()\`).”

说明：

  1. X轴 表示距离 Motif 中心的距离（bp）。
  2. Y轴 表示归一化后的 Tn5 切割频率。
  3. 如果看到中心区域（0 bp 处）有一个明显的“凹陷”（即切割频率降低），而两侧有“高峰”，这说明该转录因子在这些细胞中发生了真实的结合。
  4. 如果某一组（如 Case 组）的凹陷比另一组（如 Control 组）更深，说明该转录因子在 Case 组中的结合活性更强。
0 条评论·0 条回复