Signac 单细胞 ATAC — Seurat 生态下的单细胞染色质分析工具包¶

一句话说明¶

Signac 是与 Seurat 深度集成的 R 单细胞 ATAC-seq 分析包，复用 Seurat 的所有可视化和多组学整合功能，特别适合需要同时分析 ATAC + RNA（多组学 CITE-seq/Multiome）的场景。

安装与配置¶

# 安装 Signac（最新稳定版，截止2025年为 1.14+）
install.packages("Signac")

# 安装核心依赖
BiocManager::install(c(
  "BSgenome.Hsapiens.UCSC.hg38",    # hg38 基因组序列（基序分析用）
  "EnsDb.Hsapiens.v86",              # Ensembl 基因注释（人类 GRCh38）
  "biovizBase",
  "TFBSTools"                        # 转录因子结合位点分析
))
install.packages("Seurat")

# 安装 MACS3（peak calling）
# pip install macs3   （在终端中安装）

# 验证安装
library(Signac)
library(Seurat)
packageVersion("Signac")

核心用法¶

library(Signac)                        # 加载 Signac
library(Seurat)                        # 加载 Seurat（必须配合使用）
library(GenomeInfoDb)                  # 染色体信息处理
library(EnsDb.Hsapiens.v86)           # 人类基因注释数据库

# ── 创建 Signac 对象 ──────────────────────────────────
# 从 Cellranger ATAC 输出读取数据
counts <- Read10X_h5('filtered_peak_bc_matrix.h5')   # 读取 Peak × Cell 矩阵

# 创建染色质分析对象（ChromatinAssay）
chrom_assay <- CreateChromatinAssay(
  counts = counts,                     # Peak 计数矩阵
  sep = c(":", "-"),                   # Peak 坐标分隔符（chr:start-end）
  genome = 'hg38',                     # 参考基因组
  fragments = 'atac_fragments.tsv.gz', # Cellranger 生成的片段文件
  min.cells = 10,                      # 每个 peak 至少在 10 个细胞中有信号
  min.features = 200                   # 每个细胞至少 200 个 peak 有信号
)

# 创建 Seurat 对象（包装 ChromatinAssay）
seurat_atac <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",                     # 检测层名称
  meta.data = metadata                 # 细胞元数据（含 TSS、片段数等）
)

参数详解¶

参数	函数	说明
`sep`	`CreateChromatinAssay`	Peak 坐标分隔符（默认 `c(":", "-")`）
`genome`	`CreateChromatinAssay`	参考基因组（`hg38`/`mm10`）
`min.cells`	`CreateChromatinAssay`	Peak 过滤：最少在 N 个细胞中有信号
`n.features`	`RunTFIDF`	TF-IDF 归一化用的特征数
`reduction`	`RunSVD`	奇异值分解（相当于 ATAC 的 PCA）
`dims`	`RunUMAP`	用于 UMAP 的 SVD 维度（跳过第1维）

实战案例¶

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)

# ── 完整 Signac 分析流程 ──────────────────────────────

# 假设已创建 seurat_atac（从上面继续）

# 1. 添加基因注释信息（用于后续 Gene Activity 分析）
annotations <- GetGRangesFromEnsDb(
  ensdb = EnsDb.Hsapiens.v86,          # 使用 Ensembl v86 注释
  standard.chromosomes = TRUE          # 只保留标准染色体
)
seqlevelsStyle(annotations) <- 'UCSC' # 转换为 UCSC 格式（chr1 而非 1）
Annotation(seurat_atac) <- annotations # 添加注释到对象

# 2. 质控指标计算
# TSS 富集分数（转录起始位点附近的信号富集，高值=好细胞）
seurat_atac <- TSSEnrichment(
  object = seurat_atac,
  fast = FALSE                         # fast=TRUE 更快但不生成热图
)

# 核小体信号强度（检测核小体 ~200bp 的周期性，高值=可能 doublet）
seurat_atac <- NucleosomeSignal(object = seurat_atac)

# 读取片段数（blacklist 区域的片段比例）
seurat_atac$pct_reads_in_peaks <-    # 落在 peak 内的片段比例
  seurat_atac$peak_region_fragments / seurat_atac$passed_filters * 100
seurat_atac$blacklist_ratio <-       # 落在 blacklist 区域的比例（越低越好）
  seurat_atac$blacklist_region_fragments / seurat_atac$peak_region_fragments

# 可视化质控指标
VlnPlot(
  object = seurat_atac,
  features = c('nCount_peaks', 'TSS.enrichment', 'nucleosome_signal',
               'pct_reads_in_peaks', 'blacklist_ratio'),
  pt.size = 0.1                        # 显示散点（设为0不显示）
)

# 3. 过滤低质量细胞
seurat_atac <- subset(
  x = seurat_atac,
  subset = nCount_peaks > 3000 &       # 至少 3000 个 peak 片段
    nCount_peaks < 100000 &             # 不超过 10 万（过滤 doublet）
    pct_reads_in_peaks > 15 &           # 至少 15% 落在 peak 内
    blacklist_ratio < 0.05 &            # blacklist 比例 < 5%
    nucleosome_signal < 4 &             # 核小体信号 < 4
    TSS.enrichment > 3                  # TSS 富集 > 3
)
print(paste("过滤后细胞数：", ncol(seurat_atac)))

# 4. 归一化和降维（TF-IDF + SVD）
seurat_atac <- RunTFIDF(seurat_atac)   # TF-IDF 归一化（适合稀疏 ATAC 数据）
seurat_atac <- FindTopFeatures(
  seurat_atac,
  min.cutoff = 'q5'                    # 选 top 95% 的 peak（过滤低信号）
)
seurat_atac <- RunSVD(seurat_atac)     # SVD（相当于 PCA）

# 检查第一个维度是否与测序深度相关（如相关则跳过）
DepthCor(seurat_atac)                  # 如果 dim1 与深度相关，UMAP 跳过 dim1

# 5. 聚类和 UMAP
seurat_atac <- RunUMAP(
  object = seurat_atac,
  reduction = 'lsi',
  dims = 2:30                          # 跳过第1维（通常代表测序深度噪声）
)
seurat_atac <- FindNeighbors(
  object = seurat_atac,
  reduction = 'lsi',
  dims = 2:30
)
seurat_atac <- FindClusters(
  object = seurat_atac,
  verbose = FALSE,
  algorithm = 3,                       # Leiden 算法
  resolution = 0.5
)
DimPlot(object = seurat_atac, label = TRUE)

# 6. 基因活跃度估计（将 ATAC 信号转换为类似 RNA 表达的矩阵）
gene_activities <- GeneActivity(seurat_atac)   # 以 TSS 为中心的 ATAC 信号

# 将基因活跃度矩阵添加到对象
seurat_atac[["RNA"]] <- CreateAssayObject(counts = gene_activities)
seurat_atac <- NormalizeData(
  object = seurat_atac,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(seurat_atac$nCount_RNA)
)

# 7. Peak Calling（重新用 MACS3 call peaks，比 Cellranger 更准）
DefaultAssay(seurat_atac) <- "peaks"
seurat_atac <- CallPeaks(
  seurat_atac,
  macs2.path = '/path/to/macs3',        # MACS3 路径
  group.by = "seurat_clusters"          # 按 cluster 分组 call peaks
)

# 8. 差异可及性分析
da_peaks <- FindMarkers(
  object = seurat_atac,
  ident.1 = '0',                        # 比较 cluster 0
  ident.2 = '1',                        # 与 cluster 1
  only.pos = TRUE,                      # 只返回正向差异（更可及）
  test.use = 'LR',                      # Logistic Regression（推荐 ATAC 用）
  latent.vars = 'nCount_peaks'          # 校正测序深度的影响
)

# 显示差异 peak 对应的基因
ClosestFeature(
  object = seurat_atac,
  regions = rownames(da_peaks)[1:10]    # 查找 top10 差异 peak 最近的基因
)

常见报错与解决¶

报错	原因	解决方法
`Fragment file not found`	片段文件路径错误	检查路径，文件需要有 .tbi 索引
`BSgenome not installed`	基因组包缺失	`BiocManager::install("BSgenome.Hsapiens.UCSC.hg38")`
`NucleosomeSignal` 很慢	片段文件很大	增加 `fast=TRUE` 参数
`SVD dim1` 与深度高度相关	正常现象	UMAP 使用 `dims=2:30` 跳过 dim1
`CallPeaks` 找不到 macs	MACS3 路径错误	`which macs3` 获取正确路径

速查表¶

# 安装
install.packages("Signac")
BiocManager::install("EnsDb.Hsapiens.v86")

# 核心流程（TF-IDF 降维）
seurat_atac <- RunTFIDF(seurat_atac)
seurat_atac <- FindTopFeatures(seurat_atac, min.cutoff='q5')
seurat_atac <- RunSVD(seurat_atac)
seurat_atac <- RunUMAP(seurat_atac, reduction='lsi', dims=2:30)
seurat_atac <- FindNeighbors(seurat_atac, reduction='lsi', dims=2:30)
seurat_atac <- FindClusters(seurat_atac, resolution=0.5)

# 与 ArchR 的对比
# Signac  → Seurat 生态，适合多组学整合
# ArchR   → 独立高性能框架，适合纯 ATAC 大数据分析