177_scATAC-seq数据预处理¶

一句话概述¶

单细胞ATAC-seq（scATAC-seq）数据预处理包括比对、去重、质控、片段文件生成、peak calling和特征矩阵构建，是后续聚类和细胞类型注释的基础，主要使用CellRanger-ATAC、ArchR或Signac等工具完成。

核心知识点表格¶

知识点	说明
scATAC-seq原理	利用Tn5转座酶切割开放染色质区域，检测单细胞水平的染色质可及性
核心输出	fragments.tsv.gz片段文件
主要质控指标	TSS enrichment、片段数、核小体信号、FRiP
Peak calling	使用MACS2在聚合数据上调用peaks
特征矩阵	peak-by-cell稀疏二值/计数矩阵
降维方法	TF-IDF + LSI（Latent Semantic Indexing）
分析工具	CellRanger-ATAC、ArchR、Signac（Seurat扩展）
与scRNA-seq区别	信号稀疏、二值化特征、需要不同的标准化策略

步骤详解¶

第一步：CellRanger-ATAC上游处理¶

白话解释：CellRanger-ATAC是10x Genomics提供的官方流程，从原始FASTQ文件开始，完成比对、去重、细胞过滤和片段文件生成。就像scRNA-seq用CellRanger一样。

技术细节：CellRanger-ATAC内部使用BWA进行比对，自动过滤低质量reads和PCR重复，生成fragments.tsv.gz（每行记录一个片段的染色体、起止位置和所属barcode）。

# 下载参考基因组
cellranger-atac mkref \
    --config=refdata-config.yaml

# 或使用预构建的参考
# 下载 refdata-cellranger-arc-GRCh38-2020-A-2.0.0

# 运行CellRanger-ATAC
cellranger-atac count \
    --id=sample1 \
    --reference=/path/to/refdata-cellranger-atac-GRCh38-2020-A-2.0.0 \
    --fastqs=/path/to/fastqs/ \
    --sample=Sample1 \
    --localcores=16 \
    --localmem=64

# 主要输出文件
# outs/fragments.tsv.gz        # 片段文件（最重要）
# outs/fragments.tsv.gz.tbi    # 索引文件
# outs/filtered_peak_bc_matrix/ # peak-cell矩阵
# outs/singlecell.csv          # 单细胞指标
# outs/peaks.bed                # 调用的peaks
# outs/summary.json             # 汇总统计

第二步：使用Signac进行质控¶

白话解释：质控是决定哪些细胞数据可用的关键步骤。好的细胞应该有足够多的片段、高的TSS富集分数和清晰的核小体模式。

技术细节：关键QC指标：(1) TSS enrichment score - 转录起始位点的信号富集，反映信噪比；(2) nucleosome signal - 核小体周期性片段模式；(3) FRiP - peaks中的片段比例；(4) blacklist ratio - 黑名单区域的片段比例。

library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Hsapiens.v86)

# 加载数据
counts <- Read10X_h5("outs/filtered_peak_bc_matrix.h5")
metadata <- read.csv("outs/singlecell.csv", header = TRUE, row.names = 1)

# 创建ChromatinAssay
chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = "outs/fragments.tsv.gz",
  min.cells = 10,
  min.features = 200
)

# 创建Seurat对象
pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

# 添加基因注释
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
Annotation(pbmc) <- annotations

# 计算QC指标
pbmc <- NucleosomeSignal(pbmc)         # 核小体信号
pbmc <- TSSEnrichment(pbmc, fast = FALSE) # TSS富集

# 查看QC指标
head(pbmc@meta.data[, c("nCount_peaks", "TSS.enrichment",
                          "nucleosome_signal", "blacklist_ratio")])

第三步：质控过滤¶

# 可视化QC指标
library(ggplot2)
library(patchwork)

p1 <- VlnPlot(pbmc, features = "nCount_peaks", pt.size = 0) + NoLegend()
p2 <- VlnPlot(pbmc, features = "TSS.enrichment", pt.size = 0) + NoLegend()
p3 <- VlnPlot(pbmc, features = "nucleosome_signal", pt.size = 0) + NoLegend()
p4 <- VlnPlot(pbmc, features = "blacklist_ratio", pt.size = 0) + NoLegend()
p1 | p2 | p3 | p4

# TSS富集图（应该呈现以TSS为中心的尖峰）
TSSPlot(pbmc, group.by = "high.tss") + NoLegend()

# 片段长度分布（应看到核小体周期性：~200bp, ~400bp）
FragmentHistogram(pbmc)

# 过滤细胞
pbmc <- subset(
  pbmc,
  subset = nCount_peaks > 3000 &      # 最小片段数
    nCount_peaks < 100000 &            # 最大片段数（去除doublet）
    TSS.enrichment > 3 &              # TSS富集
    nucleosome_signal < 4 &            # 核小体信号
    blacklist_ratio < 0.01             # 黑名单比例
)

cat("过滤后细胞数:", ncol(pbmc), "\n")

第四步：标准化与降维¶

白话解释：scATAC-seq不能像scRNA-seq那样用简单的log标准化，而是用TF-IDF（词频-逆文档频率），然后用LSI（类似PCA的方法）降维。

技术细节：TF-IDF借鉴自自然语言处理：TF是每个peak在每个细胞中的频率，IDF降低普遍存在的peak的权重。LSI然后对TF-IDF矩阵做SVD分解。注意第一个LSI分量通常与测序深度相关，应该被排除。

# 标准化：TF-IDF
pbmc <- RunTFIDF(pbmc)

# 选择top features
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q0")

# LSI降维
pbmc <- RunSVD(pbmc)

# 检查LSI分量与测序深度的相关性
DepthCor(pbmc)
# 第1个分量通常与深度高度相关，聚类时应排除

# UMAP可视化
pbmc <- RunUMAP(pbmc, reduction = "lsi", dims = 2:30)  # 排除第1个分量
pbmc <- FindNeighbors(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindClusters(pbmc, algorithm = 3, resolution = 0.5)

DimPlot(pbmc, label = TRUE) + NoLegend()

第五步：使用ArchR进行预处理（替代方案）¶

白话解释：ArchR是另一个流行的scATAC-seq分析框架，使用Arrow文件存储数据，支持更大规模的数据集和更多的功能。

library(ArchR)

# 设置参考基因组
addArchRGenome("hg38")

# 创建Arrow文件
ArrowFiles <- createArrowFiles(
  inputFiles = "outs/fragments.tsv.gz",
  sampleNames = "sample1",
  minTSS = 4,            # 最低TSS enrichment
  minFrags = 1000,        # 最低片段数
  addTileMat = TRUE,      # 添加tile矩阵
  addGeneScoreMat = TRUE  # 添加基因评分矩阵
)

# 创建ArchR项目
proj <- ArchRProject(
  ArrowFiles = ArrowFiles,
  outputDirectory = "ArchR_output",
  copyArrows = TRUE
)

# Doublet检测
proj <- addDoubletScores(proj)
proj <- filterDoublets(proj)

# LSI降维
proj <- addIterativeLSI(
  proj,
  useMatrix = "TileMatrix",
  name = "IterativeLSI",
  iterations = 2,
  clusterParams = list(resolution = c(0.2), sampleCells = 10000),
  varFeatures = 25000
)

# 聚类
proj <- addClusters(proj, reducedDims = "IterativeLSI", resolution = 0.8)

# UMAP
proj <- addUMAP(proj, reducedDims = "IterativeLSI")

# 可视化
plotEmbedding(proj, colorBy = "cellColData", name = "Clusters")

实战命令速查¶

# CellRanger-ATAC快速运行
cellranger-atac count --id=sample --reference=refdata --fastqs=fastqs/ --localcores=16

# 检查fragments文件
zcat fragments.tsv.gz | head -5
# chr1    10073   10133   ACGTACGTACGTACGT-1    1

# 对fragments文件建索引（如缺失）
tabix -p bed fragments.tsv.gz

# 统计片段数
zcat fragments.tsv.gz | wc -l

# 片段数每个barcode的分布
zcat fragments.tsv.gz | cut -f4 | sort | uniq -c | sort -rn | head

# Signac快速流程
pbmc <- CreateSeuratObject(counts = CreateChromatinAssay(counts, fragments = "fragments.tsv.gz"))
pbmc <- NucleosomeSignal(pbmc)
pbmc <- TSSEnrichment(pbmc)
pbmc <- subset(pbmc, nCount_peaks > 3000 & TSS.enrichment > 3)
pbmc <- RunTFIDF(pbmc) %>% FindTopFeatures() %>% RunSVD() %>%
  RunUMAP(reduction = "lsi", dims = 2:30) %>%
  FindNeighbors(reduction = "lsi", dims = 2:30) %>%
  FindClusters(resolution = 0.5)

面试常问点¶

Q1: scATAC-seq和scRNA-seq的数据特点有什么不同？ A: scATAC-seq数据极度稀疏（每个peak在每个细胞中通常只有0或1个片段），特征空间巨大（数十万peaks），二值化特征。scRNA-seq虽然也稀疏，但基因数量有限（~2-3万），counts值可以>1。这导致两者的标准化、降维和聚类方法完全不同。

Q2: 为什么scATAC-seq用TF-IDF+LSI而不是PCA？ A: scATAC-seq的peak矩阵近似二值化，不符合PCA假设的正态分布。TF-IDF来自文本分析，适合处理稀疏的出现/缺失数据。LSI是对TF-IDF矩阵做SVD，等价于文本分析中的潜在语义分析，能有效提取主要变异模式。

Q3: TSS enrichment score反映了什么？ A: TSS enrichment score衡量转录起始位点附近的片段富集程度。高质量的ATAC-seq数据应该在TSS处有显著的信号富集（因为启动子区域通常是开放的）。低TSS enrichment表明信噪比差，可能是死细胞或技术问题。

Q4: 第一个LSI分量为什么要排除？ A: 第一个LSI分量通常与测序深度（total fragment count）高度相关，反映的是技术变异而非生物学差异。包含它会导致聚类主要由测序深度驱动。可以用DepthCor函数检查每个分量与深度的相关性。

Q5: 如何处理scATAC-seq中的doublet？ A: ArchR内置了doublet检测功能（addDoubletScores），基于模拟合成doublet的方法。Signac中可以使用AMULET工具。doublet在scATAC中比scRNA更难检测，因为数据更稀疏。

Q6: Gene activity score是什么？ A: 由于scATAC-seq不直接测量基因表达，可以通过计算每个基因body和启动子区域的ATAC片段数来估计"基因活性"。这个score可用于初步细胞类型注释，但准确性不如直接的scRNA-seq。

易错点¶

fragments文件缺少索引：Signac需要.tbi索引文件，缺失会报错。用tabix重建索引
参考基因组版本不一致：注释文件的基因组版本必须与比对使用的参考一致
过滤阈值过严或过松：TSS enrichment和片段数阈值需根据具体数据调整，不能照搬默认值
LSI dims选择不当：dims参数应排除与深度相关的分量，具体要看DepthCor结果
Peak calling策略：在所有细胞上联合call peaks还是按cluster分别call会影响结果

补充知识¶

scATAC-seq质控标准参考值¶

指标	优秀	可接受	差
TSS enrichment	>10	4-10	<4
片段数/细胞	>10,000	3,000-10,000	<3,000
FRiP	>0.3	0.15-0.3	<0.15
Blacklist ratio	<0.005	0.005-0.02	>0.02

Signac vs ArchR对比¶

特性	Signac	ArchR
框架	Seurat扩展	独立R包
数据存储	内存	Arrow文件（磁盘）
大数据支持	一般	好（适合>10万细胞）
学习曲线	低（熟悉Seurat即可）	中等
多样本整合	需手动处理	内置支持