单细胞 ATAC-seq 分析¶

一句话概述¶

scATAC-seq 在单细胞水平检测染色质可及性，揭示细胞类型特异的基因调控景观，主流分析工具为 ArchR 和 Signac。

核心知识点表格¶

知识点	说明
scATAC-seq	单细胞染色质可及性测序，10x Genomics 最常用
Fragment file	记录每个细胞每个 DNA 片段的坐标和条形码
Peak-Cell 矩阵	类似 scRNA-seq 的 count 矩阵，行=峰，列=细胞
TF-IDF	词频-逆文档频率标准化，scATAC-seq 的核心标准化方法
LSI	Latent Semantic Indexing，基于 SVD 的降维方法
Gene Activity Score	用启动子附近的 ATAC 信号估算基因"活性"
chromVAR	计算转录因子 motif 的可及性变异
ArchR	R 语言全功能 scATAC-seq 分析包（Nature Genetics 2021）
Signac	Seurat 生态的 scATAC-seq 扩展包
SnapATAC2	Rust+Python 高性能 scATAC-seq 分析工具

白话解释原理¶

想象一下： 在一个器官里有很多不同类型的细胞。每种细胞打开了不同的"基因开关"。

scATAC-seq 做的事情： 1. 把细胞一个一个分开（10x Genomics 微流控） 2. 给每个细胞的 DNA 贴上独特条形码 3. 用 Tn5 转座酶探测每个细胞中哪些位置的染色质是打开的 4. 测序后，根据条形码把数据分回到每个细胞

核心挑战： scATAC-seq 数据极其稀疏——每个细胞只能捕获几千个片段（而基因组有几百万个可能的开放位点），就像一张几乎全黑的照片，需要特殊方法来"看清"。

各步骤详解¶

第一步：CellRanger ATAC 预处理¶

# 10x Genomics CellRanger ATAC 处理原始数据
cellranger-atac count \
    --id=sample1 \                     # 样本ID
    --reference=/data/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 \  # 参考基因组
    --fastqs=/data/fastq/ \            # FASTQ目录
    --sample=sample1 \                 # FASTQ中的样本名
    --localcores=16 \                  # 线程数
    --localmem=64                      # 内存(GB)

# CellRanger ATAC 输出关键文件:
# outs/
# ├── filtered_peak_bc_matrix/     # 峰-细胞矩阵（只含有效细胞）
# │   ├── matrix.mtx               # 稀疏矩阵
# │   ├── barcodes.tsv              # 细胞条形码
# │   └── peaks.bed                 # 峰坐标
# ├── fragments.tsv.gz              # Fragment 文件（核心！）
# ├── fragments.tsv.gz.tbi          # Fragment 索引
# ├── singlecell.csv                # 单细胞QC指标
# └── summary.csv                   # 总体统计摘要

第二步：使用 Signac 分析（Seurat 生态）¶

# 安装
# BiocManager::install(c("Signac", "Seurat", "EnsDb.Hsapiens.v86",
#                         "BSgenome.Hsapiens.UCSC.hg38"))

library(Signac)                # 加载Signac
library(Seurat)                # 加载Seurat
library(EnsDb.Hsapiens.v86)    # 加载基因注释
library(ggplot2)               # 加载绑图包

# 1. 读取数据
counts <- Read10X_h5("filtered_peak_bc_matrix.h5")  # 读取峰-细胞矩阵

# 创建 ChromatinAssay
chrom_assay <- CreateChromatinAssay(
    counts = counts,                    # 计数矩阵
    sep = c(":", "-"),                   # 峰坐标分隔符
    fragments = "fragments.tsv.gz",      # Fragment文件路径
    min.cells = 10,                      # 至少10个细胞有信号的峰
    min.features = 200                   # 至少200个峰有信号的细胞
)

# 创建 Seurat 对象
pbmc <- CreateSeuratObject(
    counts = chrom_assay,               # 染色质assay
    assay = "peaks"                     # assay名称
)

# 添加基因注释
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)  # 获取注释
seqlevelsStyle(annotations) <- "UCSC"                            # 转为UCSC格式
Annotation(pbmc) <- annotations                                  # 添加到对象

# 2. 质量控制
# 计算QC指标
pbmc <- NucleosomeSignal(object = pbmc)    # 核小体信号（单核小体/无核小体比）
pbmc <- TSSEnrichment(object = pbmc)        # TSS富集分数

# 添加黑名单比例
pbmc$blacklist_ratio <- FractionCountsInRegion(
    object = pbmc,
    assay = "peaks",
    regions = blacklist_hg38               # ENCODE黑名单区域
)

# 过滤低质量细胞
pbmc <- subset(
    x = pbmc,
    subset = nCount_peaks > 1000 &         # 至少1000个片段
        nCount_peaks < 100000 &            # 最多10万个片段
        nucleosome_signal < 2 &            # 核小体信号 < 2
        TSS.enrichment > 2 &              # TSS富集 > 2
        blacklist_ratio < 0.025            # 黑名单比例 < 2.5%
)
cat("过滤后细胞数:", ncol(pbmc), "\n")     # 打印细胞数

# 3. 标准化和降维
# TF-IDF 标准化
pbmc <- RunTFIDF(pbmc)                      # Term Frequency - Inverse Document Frequency

# 特征选择（选择高变异峰）
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q0")  # 使用所有峰

# LSI 降维（类似 PCA，但适合稀疏二值数据）
pbmc <- RunSVD(pbmc)                         # SVD降维
# 注意：第1个LSI成分通常与测序深度相关，需要排除
DepthCor(pbmc)                               # 检查每个LSI成分与深度的相关性

# 4. 聚类
pbmc <- RunUMAP(
    object = pbmc,
    reduction = "lsi",                       # 使用LSI降维结果
    dims = 2:30                              # 排除第1个成分
)
pbmc <- FindNeighbors(
    object = pbmc,
    reduction = "lsi",
    dims = 2:30
)
pbmc <- FindClusters(
    object = pbmc,
    algorithm = 3,                           # SLM算法
    resolution = 0.5                         # 分辨率
)

# 可视化
DimPlot(pbmc, label = TRUE) + NoLegend()     # UMAP图

第三步：Gene Activity 和细胞类型注释¶

# 1. 计算 Gene Activity Score
# 用启动子上游2kb + 基因体的ATAC信号估算基因"活性"
gene_activities <- GeneActivity(pbmc)       # 计算基因活性

# 添加到Seurat对象
pbmc[["RNA"]] <- CreateAssayObject(counts = gene_activities)  # 作为新assay
pbmc <- NormalizeData(                       # 标准化
    object = pbmc,
    assay = "RNA",
    normalization.method = "LogNormalize"
)

# 2. 与 scRNA-seq 参考数据整合注释
# 加载 scRNA-seq 参考
rna_ref <- readRDS("pbmc_rna_reference.rds")  # scRNA-seq参考数据

# 找到锚点
transfer_anchors <- FindTransferAnchors(
    reference = rna_ref,                      # RNA参考
    query = pbmc,                              # ATAC查询
    reduction = "cca",                         # CCA方法
    dims = 1:30
)

# 转移标签
predicted_labels <- TransferData(
    anchorset = transfer_anchors,             # 锚点
    refdata = rna_ref$celltype,               # 参考的细胞类型
    weight.reduction = pbmc[["lsi"]],         # 权重
    dims = 2:30
)

pbmc <- AddMetaData(object = pbmc,            # 添加预测标签
                     metadata = predicted_labels)
DimPlot(pbmc, group.by = "predicted.id",       # 按预测细胞类型着色
        label = TRUE) + NoLegend()

第四步：Motif 分析和差异可及性¶

# 1. Motif 富集分析（chromVAR）
library(chromVAR)         # 加载chromVAR
library(JASPAR2020)       # JASPAR motif数据库
library(TFBSTools)        # TF结合位点工具

# 获取motif矩阵
pfm <- getMatrixSet(
    x = JASPAR2020,                           # JASPAR数据库
    opts = list(collection = "CORE",           # 核心集合
                tax_group = "vertebrates",     # 脊椎动物
                all_versions = FALSE)          # 只要最新版
)

# 添加 motif 信息
pbmc <- AddMotifs(
    object = pbmc,
    genome = BSgenome.Hsapiens.UCSC.hg38,     # 基因组序列
    pfm = pfm                                  # motif矩阵
)

# 运行 chromVAR
pbmc <- RunChromVAR(
    object = pbmc,
    genome = BSgenome.Hsapiens.UCSC.hg38      # 基因组
)
# 现在可以用 FeaturePlot 查看特定TF的motif可及性
FeaturePlot(pbmc, features = "MA0139.1",       # CTCF motif
            min.cutoff = "q10", max.cutoff = "q90")

# 2. 差异可及性分析
# 找差异可及峰（类似 FindMarkers）
da_peaks <- FindMarkers(
    object = pbmc,
    ident.1 = "CD4 T",                        # 细胞类型1
    ident.2 = "CD8 T",                        # 细胞类型2
    test.use = "LR",                           # logistic回归
    latent.vars = "nCount_peaks",              # 控制测序深度
    min.pct = 0.05                             # 至少5%的细胞有信号
)
# 在差异峰中做motif富集
enriched_motifs <- FindMotifs(
    object = pbmc,
    features = rownames(da_peaks[da_peaks$p_val_adj < 0.05 &
                                 da_peaks$avg_log2FC > 0.5, ])
)
head(enriched_motifs, 10)                      # 查看top 10 motif

第五步：使用 ArchR 分析（替代方案）¶

# ArchR 的优势: 更快、内存更低、功能更全面、处理大数据集
# install.packages("devtools")
# devtools::install_github("GreenleafLab/ArchR")

library(ArchR)  # 加载ArchR

# 1. 创建 Arrow 文件（ArchR的核心数据格式，基于HDF5）
ArrowFiles <- createArrowFiles(
    inputFiles = "fragments.tsv.gz",       # Fragment文件
    sampleNames = "sample1",               # 样本名
    minTSS = 4,                            # 最低TSS富集分数
    minFrags = 1000,                       # 最少片段数
    addTileMat = TRUE,                     # 添加500bp Tile矩阵
    addGeneScoreMat = TRUE                 # 添加Gene Score矩阵
)

# 2. 创建 ArchRProject
proj <- ArchRProject(
    ArrowFiles = ArrowFiles,               # Arrow文件
    outputDirectory = "ArchR_output",       # 输出目录
    copyArrows = TRUE                      # 复制Arrow文件
)

# 3. 自动分析流程
proj <- addDoubletScores(proj)              # 双细胞检测
proj <- filterDoublets(proj)                # 过滤双细胞
proj <- addIterativeLSI(proj)               # 迭代LSI降维
proj <- addClusters(proj, resolution = 0.8) # 聚类
proj <- addUMAP(proj)                       # UMAP可视化
proj <- addGroupCoverages(proj)             # 计算组别覆盖度
proj <- addReproduciblePeakSet(proj)        # 调用可重复峰
proj <- addMotifAnnotations(proj)           # 添加motif注释

# 4. 可视化
plotEmbedding(proj, colorBy = "Clusters")   # UMAP按cluster着色

常见报错与解决¶

报错	原因	解决方案
`Low median fragments per cell (<3000)`	文库质量不佳	检查实验步骤，考虑降低过滤阈值
`First LSI component correlated with depth`	正常现象	在下游分析中排除第1个成分（dims = 2:30）
`CreateChromatinAssay: fragments file error`	fragments文件未建索引	用 tabix 建索引: `tabix -p bed fragments.tsv.gz`
`Memory error with ArchR`	数据量大	ArchR 设计上支持大数据，检查 addTileMat 参数
`Poor clustering resolution`	数据稀疏或细胞数太少	调整 resolution 参数，或使用更多峰
`GeneActivity scores not matching RNA`	基因活性只是粗略估计	通过与 scRNA-seq 整合来准确注释

速查表¶

┌──────────────────────────────────────────────────────────┐
│              scATAC-seq 分析速查                          │
├──────────────────────────────────────────────────────────┤
│ 工具对比:                                                 │
│   Signac — Seurat生态, R语言, 灵活                       │
│   ArchR — 全功能, 处理大数据, HDF5存储                    │
│   SnapATAC2 — Rust+Python, 最快                          │
│                                                          │
│ 分析流程:                                                 │
│   CellRanger ATAC → 读取 → QC → TF-IDF →                │
│   LSI → 聚类 → UMAP → Gene Activity →                   │
│   细胞注释 → Motif → 差异可及性                          │
│                                                          │
│ QC 标准:                                                 │
│   片段数: >1000 per cell                                 │
│   TSS 富集: >2 (最低), >4 (推荐)                        │
│   核小体信号: <2                                         │
│   黑名单比例: <2.5%                                      │
│                                                          │
│ 关键步骤:                                                 │
│   标准化: TF-IDF (非 LogNormalize!)                      │
│   降维: LSI/SVD (非 PCA!)                                │
│   第1个LSI成分通常需要排除                                │
│                                                          │
│ 细胞注释: Label Transfer 从 scRNA-seq                    │
│ Motif 分析: chromVAR + JASPAR                            │
└──────────────────────────────────────────────────────────┘

面试高频问题¶

scATAC-seq 数据为什么比 scRNA-seq 更稀疏？怎么处理？ 答：scRNA-seq 中一个基因有多个 mRNA 拷贝可被捕获，而 scATAC-seq 中一个位点最多只有 2 个 DNA 拷贝（二倍体），实际捕获率更低，导致每个细胞只有几千个非零值。处理策略：(1) 使用 TF-IDF 标准化而非 log-normalize；(2) 降维用 LSI/SVD 而非 PCA；(3) 使用更大的特征集（如 500bp tiles 而非 peaks）；(4) 利用 Gene Activity Score 聚合到基因水平降低稀疏度。
TF-IDF 标准化是什么？为什么 scATAC-seq 要用它？ 答：TF-IDF 来自自然语言处理，分两步：(1) TF（term frequency）——对每个细胞内的信号除以总片段数，校正测序深度；(2) IDF（inverse document frequency）——对每个峰乘以 log(总细胞数/有此峰的细胞数)，给稀有峰更高权重（因为稀有的开放位点更能区分细胞类型）。这比简单的 log-normalize 更适合二值化的稀疏数据。
如何将 scATAC-seq 数据与 scRNA-seq 整合？ 答：主要方法有两种。(1) Label transfer（Signac/Seurat）：先计算 Gene Activity Score 作为 scATAC 的"伪表达"，再用 CCA 或 RPCA 找到两个数据集之间的锚点，将 scRNA 的细胞标签转移到 scATAC。(2) WNN（Weighted Nearest Neighbor）：如果是 10x Multiome 的配对数据（同一个细胞同时有 RNA 和 ATAC），可以用 WNN 学习每个细胞中两种数据的相对权重，构建加权最近邻图。ArchR 也有内置的整合方法。