单细胞ATAC-seq分析（scATAC-seq）¶

一句话概述¶

单细胞ATAC-seq通过检测单个细胞中的染色质开放区域揭示基因调控状态，核心分析工具ArchR和Signac可完成从fragment到peak calling、降维聚类、peak-gene linking和motif分析的完整流程。

核心知识点表格¶

知识点	说明
ATAC-seq原理	Tn5转座酶切割开放染色质，测序插入位点反映可及性
Fragment file	scATAC-seq核心输入文件，记录每个barcode的fragment坐标
ArchR	R语言的综合scATAC-seq分析框架，内存效率高
Signac	Seurat生态的scATAC-seq扩展包
Peak calling	从aggregate数据识别显著富集的开放区域
Gene Activity Score	通过启动子/基因体附近的可及性估计"基因表达"
Peak-Gene Linking	将远端peak与目标基因关联（类似enhancer-gene）
Motif enrichment	在差异peak中寻找富集的转录因子结合基序
LSI/TF-IDF	Latent Semantic Indexing，scATAC-seq标准降维方法
ChromVAR	估计每个细胞中TF motif的可及性偏差

各步骤详解¶

第一步：数据预处理与质控¶

白话解释： scATAC-seq数据通常经Cell Ranger ATAC处理后得到fragment文件。质控关注每个细胞的fragment数量、TSS富集度（好的ATAC数据应该在基因启动子区域富集）、以及核小体条带模式。

技术细节： - 输入：fragments.tsv.gz（Cell Ranger ATAC输出） - 质控指标：(1) nFragments > 1000-3000；(2) TSS enrichment > 4-8；(3) Nucleosome ratio < 4；(4) Blacklist ratio < 0.05 - Fragment size分布应显示核小体周期性模式（~200bp, ~400bp峰） - TSS enrichment反映信噪比

代码示例（ArchR）：

library(ArchR)

# 设置线程和基因组
addArchRThreads(threads = 16)
addArchRGenome("hg38")

# 创建Arrow文件（ArchR核心数据结构）
inputFiles <- c("sample1/outs/fragments.tsv.gz",
                "sample2/outs/fragments.tsv.gz")
names(inputFiles) <- c("Sample1", "Sample2")

ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  minTSS = 4,           # 最低TSS enrichment
  minFrags = 1000,      # 最低fragment数
  addTileMat = TRUE,    # 添加tile矩阵（500bp bins）
  addGeneScoreMat = TRUE # 添加gene activity score
)

# 质控可视化
# TSS enrichment vs nFragments
df <- getCellColData(ArchRProj, select = c("TSSEnrichment", "nFrags"))
ggplot(as.data.frame(df), aes(x = log10(nFrags), y = TSSEnrichment)) +
  geom_hex(bins = 100) +
  theme_minimal() +
  geom_hline(yintercept = 4, lty = "dashed") +
  geom_vline(xintercept = log10(1000), lty = "dashed")

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

代码示例（Signac）：

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

# 读取10x数据
counts <- Read10X_h5("sample1/outs/filtered_peak_bc_matrix.h5")
metadata <- read.csv("sample1/outs/singlecell.csv", row.names = 1)
fragpath <- "sample1/outs/fragments.tsv.gz"

# 获取基因注释
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"

# 创建ChromatinAssay
chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = fragpath,
  annotation = annotations
)

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

# 质控
obj <- NucleosomeSignal(obj)
obj <- TSSEnrichment(obj)

# 过滤
obj <- subset(obj,
  nCount_peaks > 3000 &
  nCount_peaks < 100000 &
  nucleosome_signal < 4 &
  TSS.enrichment > 2
)

第二步：降维与聚类¶

白话解释： scATAC-seq数据是极度稀疏的二值矩阵（每个位点在每个细胞中要么开放要么关闭）。标准方法是用TF-IDF加权后做SVD（类似PCA），这叫LSI（Latent Semantic Indexing）。然后就和scRNA-seq一样做UMAP和聚类。

技术细节： - TF-IDF：Term Frequency-Inverse Document Frequency，给稀有特征更高权重 - LSI components通常取2-30（去掉第1个component，它通常反映测序深度） - 聚类使用graph-based方法（与Seurat类似） - ArchR使用iterative LSI进一步优化

代码示例（ArchR）：

# Iterative LSI降维
proj <- addIterativeLSI(
  ArchRProj = proj,
  useMatrix = "TileMatrix",
  name = "IterativeLSI",
  iterations = 2,
  clusterParams = list(resolution = c(0.2), sampleCells = 10000, n.start = 10),
  varFeatures = 25000,
  dimsToUse = 1:30
)

# 添加Harmony batch correction（如有多样本）
proj <- addHarmony(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  name = "Harmony",
  groupBy = "Sample"
)

# 聚类
proj <- addClusters(
  input = proj,
  reducedDims = "IterativeLSI",
  method = "Seurat",
  name = "Clusters",
  resolution = 0.8
)

# UMAP
proj <- addUMAP(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  name = "UMAP",
  nNeighbors = 30,
  minDist = 0.5
)

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

代码示例（Signac）：

# 归一化 + LSI
obj <- RunTFIDF(obj)
obj <- FindTopFeatures(obj, min.cutoff = "q0")  # 使用所有features
obj <- RunSVD(obj)

# 检查component与深度的相关性（去掉高度相关的）
DepthCor(obj)  # 通常去掉component 1

# 聚类
obj <- RunUMAP(obj, reduction = "lsi", dims = 2:30)
obj <- FindNeighbors(obj, reduction = "lsi", dims = 2:30)
obj <- FindClusters(obj, algorithm = 3, resolution = 0.5)

# 可视化
DimPlot(obj, label = TRUE) + NoLegend()

第三步：Peak Calling与细胞类型注释¶

白话解释： 之前用的是Cell Ranger给的peak或tile矩阵。现在可以按cluster分组做更精确的peak calling——每个细胞类型有自己特异的开放区域。细胞类型注释可以通过Gene Activity Score（类似RNA表达）或与配对scRNA-seq数据整合来实现。

技术细节： - ArchR使用MACS2对每个cluster的pseudo-bulk数据做peak calling - Signac也调用MACS2 - Gene Activity Score：累计基因启动子±2kb和gene body的可及性信号 - 与scRNA-seq整合：label transfer（Seurat/ArchR均支持）

代码示例（ArchR）：

# 基于聚类的peak calling
proj <- addGroupCoverages(ArchRProj = proj, groupBy = "Clusters")
proj <- addReproduciblePeakSet(
  ArchRProj = proj,
  groupBy = "Clusters",
  pathToMacs2 = findMacs2()  # 需要安装MACS2
)
proj <- addPeakMatrix(proj)

# Gene Activity Score可视化（辅助细胞类型注释）
markerGenes <- c("CD34", "CD14", "CD3D", "MS4A1", "GATA1")
plotEmbedding(proj,
  colorBy = "GeneScoreMatrix",
  name = markerGenes,
  embedding = "UMAP",
  quantCut = c(0.01, 0.95))

# 与scRNA-seq整合做label transfer
# 假设已有seRNA（Seurat对象）
proj <- addGeneIntegrationMatrix(
  ArchRProj = proj,
  useMatrix = "GeneScoreMatrix",
  matrixName = "GeneIntegrationMatrix",
  reducedDims = "IterativeLSI",
  seRNA = seRNA,                    # scRNA-seq Seurat对象
  addToArrow = FALSE,
  groupRNA = "celltype",            # scRNA-seq的细胞类型列
  nameCell = "predictedCell",
  nameGroup = "predictedGroup",
  nameScore = "predictedScore"
)

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

白话解释： 找出不同细胞类型/条件之间开放程度显著不同的peak区域，类似于scRNA-seq中的差异表达基因分析。

代码示例（ArchR）：

# Marker peaks（每个cluster vs 所有其他cluster）
markerPeaks <- getMarkerFeatures(
  ArchRProj = proj,
  useMatrix = "PeakMatrix",
  groupBy = "Clusters",
  bias = c("TSSEnrichment", "log10(nFrags)"),  # 校正技术偏差
  testMethod = "wilcoxon"
)

# 获取marker peak列表
markerList <- getMarkers(markerPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")

# 热图展示
heatmapPeaks <- plotMarkerHeatmap(
  seMarker = markerPeaks,
  cutOff = "FDR <= 0.01 & Log2FC >= 1",
  transpose = TRUE
)
draw(heatmapPeaks)

# 浏览器轨迹图（类似IGV）
plotBrowserTrack(
  ArchRProj = proj,
  groupBy = "Clusters",
  geneSymbol = "CD14",
  upstream = 50000,
  downstream = 50000
)

第五步：Peak-Gene Linking¶

白话解释： Peak-gene linking试图找出远端的开放区域（可能是enhancer）和它们调控的目标基因之间的关系。原理是：如果一个peak的可及性和某基因的表达（或gene activity score）在不同细胞中高度相关，那它们可能有调控关系。

技术细节： - 基于peak可及性和gene activity/expression的相关性 - 通常限制在一定距离内（如±250kb） - 需要足够的细胞数量保证统计力 - 结果类似于eQTL但在单细胞水平 - 可结合Hi-C数据验证

代码示例（ArchR）：

# Peak-to-gene linking
proj <- addPeak2GeneLinks(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  useMatrix = "GeneIntegrationMatrix"  # 如果有配对RNA-seq
)

# 获取peak-gene links
p2g <- getPeak2GeneLinks(
  ArchRProj = proj,
  corCutOff = 0.45,
  resolution = 1,
  returnLoops = FALSE
)

# 可视化特定基因的linked peaks
plotBrowserTrack(
  ArchRProj = proj,
  groupBy = "Clusters",
  geneSymbol = "MYC",
  upstream = 200000,
  downstream = 200000,
  loops = getPeak2GeneLinks(proj, corCutOff = 0.45, resolution = 100, returnLoops = TRUE)
)

# 基因组轨迹上画loop连线
plotPeak2GeneHeatmap(ArchRProj = proj, groupBy = "Clusters")

代码示例（Signac）：

# Signac中的peak-gene linking
# 需要有配对的RNA assay
DefaultAssay(obj) <- "peaks"

# 链接peaks to genes
obj <- LinkPeaks(
  object = obj,
  peak.assay = "peaks",
  expression.assay = "RNA",  # 需要配对RNA数据
  genes.use = c("MYC", "CD14", "MS4A1")
)

# 可视化
CoveragePlot(
  object = obj,
  region = "MYC",
  features = "MYC",
  expression.assay = "RNA",
  extend.upstream = 200000,
  extend.downstream = 50000
)

第六步：Motif分析与转录因子活性¶

白话解释： 在差异开放的peak区域中搜索转录因子(TF)结合基序(motif)的富集情况，推断哪些TF在特定细胞类型中活跃。ChromVAR更进一步，为每个细胞计算每个motif的可及性偏差分数。

技术细节： - Motif数据库：JASPAR、CIS-BP、ENCODE - Motif enrichment：超几何检验或Fisher's exact test - ChromVAR：计算每个细胞对每个motif的deviation z-score - 可识别"pioneer factors"（在关闭染色质中也结合的TF）

代码示例（ArchR）：

# 添加motif注释
proj <- addMotifAnnotations(ArchRProj = proj, motifSet = "cisbp", name = "Motif")

# Motif enrichment in marker peaks
enrichMotifs <- peakAnnoEnrichment(
  seMarker = markerPeaks,
  ArchRProj = proj,
  peakAnnotation = "Motif",
  cutOff = "FDR <= 0.01 & Log2FC >= 1"
)

# 热图展示
heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 7, transpose = TRUE)
draw(heatmapEM)

# ChromVAR：单细胞TF活性
proj <- addBgdPeaks(proj)
proj <- addDeviationsMatrix(
  ArchRProj = proj,
  peakAnnotation = "Motif",
  force = TRUE
)

# 可视化特定TF的deviation score
plotEmbedding(
  ArchRProj = proj,
  colorBy = "MotifMatrix",
  name = c("z:GATA1_383", "z:SPI1_322", "z:PAX5_657"),
  embedding = "UMAP",
  imputeWeights = getImputeWeights(proj)
)

# TF footprinting
motifPositions <- getPositions(proj)
proj <- addFootprints(ArchRProj = proj, positions = motifPositions, groupBy = "Clusters")

plotFootprints(
  seFoot = getFootprints(proj),
  ArchRProj = proj,
  normMethod = "Subtract",
  plotName = "Footprints-Subtract",
  addDOC = FALSE,
  smoothWindow = 5
)

实战命令（可复制）¶

# ===== 环境安装 =====
# R包
# install.packages("devtools")
# devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())
# BiocManager::install(c("Signac", "Seurat", "chromVAR", "motifmatchr"))

# MACS2/MACS3（peak calling必需）
# 推荐安装MACS3（支持scATAC-seq的FRAG格式和barcode过滤，且内置HMMRATAC用于ATAC-seq）
pip install macs3
# 或
conda install -c bioconda macs3
# MACS2仍可用但已停止新功能开发，仅修复安装问题
# pip install macs2

# ===== Cell Ranger ATAC（单独ATAC-seq实验）=====
cellranger-atac count \
  --id=sample1 \
  --reference=/ref/refdata-cellranger-atac-GRCh38-2020-A-2.0.0 \
  --fastqs=/data/fastqs/sample1/ \
  --localcores=16 \
  --localmem=64

# ===== Cell Ranger ARC（Multiome ATAC+GEX联合实验）=====
# 注意：Cell Ranger ARC 用于 10x Multiome 试剂盒（同时测 ATAC + RNA），
# Cell Ranger ATAC 用于单独的 scATAC-seq 试剂盒，两者使用不同的参考基因组包
# cellranger-arc count \
#   --id=sample1 \
#   --reference=/ref/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 \
#   --libraries=libraries.csv \
#   --localcores=16 \
#   --localmem=64

# ===== ArchR完整流程 =====
# 见上文R代码各步骤

# ===== Signac完整流程 =====
# 见上文R代码各步骤

面试常问点¶

Q1: scATAC-seq和bulk ATAC-seq的主要区别？¶

A: (1) 分辨率：scATAC-seq在单细胞水平，能区分不同细胞类型的调控状态；bulk是所有细胞的平均信号。(2) 数据特点：scATAC-seq极度稀疏（每个细胞仅覆盖基因组的2-5%），需要特殊的降维和归一化方法。(3) Peak calling：scATAC-seq通常先聚类再对pseudo-bulk做peak calling。(4) 分析框架不同（LSI vs 标准peak分析）。

Q2: 为什么scATAC-seq用LSI而不是PCA？¶

A: scATAC-seq数据是极度稀疏的近似二值矩阵，不满足PCA对连续高斯分布的假设。LSI来自自然语言处理（文档-词频分析），通过TF-IDF加权强调信息性特征（稀有但有区分度的peak），然后做SVD降维。这更适合稀疏计数数据的结构。

Q3: ArchR和Signac如何选择？¶

A: ArchR：(1) 内存效率更高（使用Arrow文件格式，不需全量加载）；(2) 功能更全面（内置peak-gene linking、footprinting、轨迹分析、doublet检测）；(3) 适合大数据集（>100k细胞）。Signac：(1) 与Seurat生态无缝整合；(2) 如果同时分析scRNA-seq+scATAC-seq（multiome），Signac更方便；(3) 学习曲线较低。此外还有SnapATAC2（Python/Rust），2024年独立基准测试（Genome Biology）显示它在库大小偏差校正和聚类上表现优于ArchR和Signac，适合Python生态用户。

Q4: Gene Activity Score准确吗？¶

A: Gene Activity Score是近似估计，有局限性：(1) 假设启动子/gene body可及性与表达正相关，但这并非总是成立；(2) 对远端增强子调控的基因估计不准；(3) 分辨率有限。10x Multiome（同时测RNA+ATAC）可以提供真实的配对数据。建议将Gene Activity Score主要用于细胞类型注释的辅助手段，不作为定量表达估计。

Q5: Peak-Gene Linking的原理和局限？¶

A: 原理：在跨细胞层面计算peak可及性与基因表达/score的相关性（Pearson/Spearman），距离限制在±250-500kb。局限：(1) 相关性不等于因果（可能是间接关联）；(2) 受细胞数量影响大；(3) 无法检测跨染色体调控；(4) 需要Hi-C或CRISPR扰动实验验证。

Q6: ChromVAR的原理和用途？¶

A: ChromVAR为每个细胞、每个motif计算一个deviation z-score：观察到该motif相关peak的可及性是否显著偏离基于GC含量等预期的背景。用途：(1) 识别不同细胞类型中活跃的TF；(2) 聚类和可视化辅助；(3) 比较条件间TF活性变化。优点是不需要peak calling，直接在fragment数据上计算。

Q7: 如何处理scATAC-seq的batch effect？¶

A: (1) ArchR内置Harmony整合；(2) Signac可结合Seurat的IntegrateData或Harmony；(3) 在LSI空间中应用Harmony/LIGER等方法；(4) 确保不同batch使用相同的feature set（如统一的peak set）。

易错点¶

1. LSI中未去除第一个component¶

问题： LSI的第一个component通常高度相关于测序深度（技术因素），保留它会让聚类结果被深度主导。 正确做法： 使用dims 2:30（ArchR中自动处理，Signac中需手动指定）。验证：用DepthCor()检查。

2. Peak calling时未分组¶

问题： 对所有细胞的aggregate数据做单一peak calling，可能被优势细胞类型主导，丢失少数群体的peaks。 正确做法： 按cluster/celltype分组做peak calling（如ArchR的addReproduciblePeakSet），再合并去冗余得到union peak set。

3. 质控阈值过于宽松¶

问题： 混入了大量低质量细胞（低fragment数、低TSS enrichment），影响下游分析。 正确做法： 严格质控：nFrags>3000，TSS enrichment>6-8，nucleosome_signal<4，blacklist_fraction<0.01。宁可丢失一些真细胞也不要混入噪声。

4. 忽略fragment文件索引¶

问题： Signac需要fragment file的tabix索引（.tbi），缺失会报错或极慢。 正确做法： tabix -p bed fragments.tsv.gz 生成索引文件。

5. Motif分析解读过度¶

问题： Motif enrichment只说明序列中存在该TF的结合基序，不代表TF真的结合或有功能。 正确做法： 结合TF表达数据（如配对RNA-seq中TF是否表达）和footprinting分析（真正结合会留下保护印记）来验证。

6. 将scATAC聚类直接等同于scRNA聚类¶

问题： 两种数据的聚类结果可能不完全一致，因为染色质状态和转录不完全同步。 正确做法： 如有multiome数据用WNN (Weighted Nearest Neighbor)整合；否则用label transfer但保持对差异的关注。

补充知识¶

10x Multiome（联合scRNA+scATAC）¶

# Seurat/Signac处理Multiome数据
# 同一个细胞同时有RNA和ATAC数据
library(Seurat)
library(Signac)

# 读取multiome数据
inputdata.10x <- Read10X_h5("filtered_feature_bc_matrix.h5")
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# 创建multimodal对象
obj <- CreateSeuratObject(counts = rna_counts)
obj[["ATAC"]] <- CreateChromatinAssay(counts = atac_counts, ...)

# WNN整合
obj <- FindMultiModalNeighbors(obj, reduction.list = list("pca", "lsi"), dims.list = list(1:30, 2:30))
obj <- RunUMAP(obj, nn.name = "weighted.nn", reduction.name = "wnn.umap")
obj <- FindClusters(obj, graph.name = "wsnn", resolution = 0.8)

单细胞ATAC-seq分析（scATAC-seq）¶

一句话概述¶

核心知识点表格¶

各步骤详解¶

第一步：数据预处理与质控¶

第二步：降维与聚类¶

第三步：Peak Calling与细胞类型注释¶

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

第五步：Peak-Gene Linking¶

第六步：Motif分析与转录因子活性¶

实战命令（可复制）¶

面试常问点¶

Q1: scATAC-seq和bulk ATAC-seq的主要区别？¶

Q2: 为什么scATAC-seq用LSI而不是PCA？¶

Q3: ArchR和Signac如何选择？¶

Q4: Gene Activity Score准确吗？¶

Q5: Peak-Gene Linking的原理和局限？¶

Q6: ChromVAR的原理和用途？¶

Q7: 如何处理scATAC-seq的batch effect？¶

易错点¶

1. LSI中未去除第一个component¶

2. Peak calling时未分组¶

3. 质控阈值过于宽松¶

4. 忽略fragment文件索引¶

5. Motif分析解读过度¶

6. 将scATAC聚类直接等同于scRNA聚类¶

补充知识¶

10x Multiome（联合scRNA+scATAC）¶

最新进展¶

📚 相关文章推荐