在本教程中，我们将分析由10x Genomics 官网提供的人类外周血单核细胞（PBMCs）的单细胞ATAC-seq数据集。本指南中使用的以下文件均可通过10x Genomics网站获取：
原始数据：The [Raw data]
元数据：The [Metadata]
片段文件：The [fragments file]
片段文件索引：The fragments file [index]
下载方式：
wget https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz.tbi
- 加载必要的数据包
library(Signac)
library(Seurat)
library(GenomicRanges)
## Warning: package 'GenomicRanges' was built under R version 4.5.2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(patchwork)
library(GenomeInfoDb)
预处理流程
在预处理染色质数据时，Signac 使用了两个相关输入文件的信息，这两个文件都可以通过 CellRanger 创建：
峰/细胞矩阵。这类似于用于分析单细胞 RNA-seq 的基因表达计数矩阵。但不同的是，矩阵的每一行代表基因组中的一个区域（即一个峰），该区域被预测为开放染色质区域。矩阵中的每个数值代表每个单一条形码（即一个细胞）在每个峰内定位的 Tn5 整合位点的数量。您可以在 10X 网站上找到更多详细信息。
片段文件。该文件列出了所有单细胞中全部独特片段的信息。它是一个大得多的文件，处理速度较慢，并存储在磁盘上（而不是内存中）。不过，保留此文件的优势在于，它包含了与每个单细胞相关的所有片段，而不仅仅是那些映射到峰区域的片段。关于片段文件的更多信息，可以在 10x Genomics 网站或 sinto 网站上找到。
我们首先使用 cellranger-atac 生成的峰/细胞矩阵和细胞元数据创建一个 Seurat 对象，并将磁盘上片段文件的路径存储在 Seurat 对象中。
counts <- Read10X_h5(filename = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv",
header = TRUE,
row.names = 1
)
chrom_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
fragments = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz",
min.cells = 10,
min.features = 200
)
pbmc <- CreateSeuratObject(
counts = chrom_assay,
assay = "peaks",
meta.data = metadata
)
如果没有 H5 文件怎么办？
ATAC-seq 数据使用一种自定义的测定对象——ChromatinAssay 来存储。这使得一些专门用于分析基因组单细胞测定（如 scATAC-seq）的功能得以实现。通过打印该测定对象，我们可以看到 ChromatinAssay 中可以包含的一些额外信息，包括基序信息、基因注释和基因组信息。
pbmc[['peaks']]
例如，我们可以对将 ChromatinAssay 设为活跃测定对象的 Seurat 对象（或直接对 ChromatinAssay 对象）调用 granges 函数，以查看与该对象中每个特征相关联的基因组范围。有关 ChromatinAssay 类的更多信息，请参阅对象交互指南。
granges(pbmc)
然后，我们移除那些对应于染色体骨架（例如 KI270713.1）或其他非标准序列的特征，只保留（22 + 2）条标准染色体。
peaks.keep <- seqnames(granges(pbmc)) %in% standardChromosomes(granges(pbmc))
pbmc <- pbmc[as.vector(peaks.keep), ]
我们还可以向 pbmc 对象添加人类基因组的基因注释信息。这样，下游的函数就可以直接从该对象中获取基因注释信息。
每个基因组组装版本都会发布多个补丁版本。在处理已比对的数据（例如我们将要使用的 10x Genomics 文件）时，建议使用与进行比对时相同的组装补丁版本对应的注释信息。
从数据集摘要中可以看出，10x Genomics 用于进行比对的参考包是“GRCh38-2020-A”，这对应的是 Ensembl v98 补丁版本。有关各个 Ensembl 版本的更多信息，您可以参考此网站。
library(AnnotationHub)
ah <- AnnotationHub()
# Search for the Ensembl 98 EnsDb for Homo sapiens on AnnotationHub
query(ah, "EnsDb.Hsapiens.v98")
## AnnotationHub with 1 record
## # snapshotDate(): 2025-10-29
## # names(): AH75011
## # $dataprovider: Ensembl
## # $species: Homo sapiens
## # $rdataclass: EnsDb
## # $rdatadateadded: 2019-05-02
## # $title: Ensembl 98 EnsDb for Homo sapiens
## # $description: Gene and protein annotations for Homo sapiens based on Ensembl version 98.
## # $taxonomyid: 9606
## # $genome: GRCh38
## # $sourcetype: ensembl
## # $sourceurl: http://www.ensembl.org
## # $sourcesize: NA
## # $tags: c("98", "AHEnsDbs", "Annotation", "EnsDb", "Ensembl", "Gene", "Protein", "Transcript")
## # retrieve record with 'object[["AH75011"]]'
ensdb_v98 <- ah[["AH75011"]]
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = ensdb_v98)
# change to UCSC style since the data was mapped to hg38
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"
# add the gene information to the object
Annotation(pbmc) <- annotations
计算质控指标
现在我们可以为 scATAC-seq 实验计算一些质控指标。目前建议使用以下指标来评估数据质量。与 scRNA-seq 类似，这些参数的预期范围会因生物学系统、细胞活力和其他因素而异。
核小体条带模式：DNA 片段大小（由双端测序读段确定）的直方图应显示出强烈的核小体条带模式，该模式对应于包裹在单个核小体周围的 DNA 长度。我们为每个单细胞计算这一指标，并量化单核小体片段与无核小体片段的近似比值（存储为 nucleosome_signal）。
转录起始位点（TSS）富集得分：ENCODE 项目定义了一个基于 TSS 中心片段与 TSS 侧翼区域片段比值的 ATAC-seq 靶向评分（参见 https://www.encodeproject.org/data-standards/terms/）。质量较差的 ATAC-seq 实验通常 TSS 富集得分较低。我们可以使用 TSSEnrichment() 函数为每个细胞计算该指标，结果将存储在元数据中的 TSS.enrichment 列名下。
峰内片段总数：衡量细胞测序深度/复杂度的指标。由于测序深度过低，需要排除读段数非常少的细胞。而读段数极高的细胞可能代表双细胞、细胞核团块或其他假象。
峰内片段比例：表示落在 ATAC-seq 峰区域内的所有片段的比例。该值较低的细胞（例如 <15-20%）通常代表低质量细胞或应去除的技术假象。需要注意的是，该值可能对所使用的一组峰区域敏感。
基因组黑名单区域中的读段比例：ENCODE 项目提供了一份黑名单区域列表，这些区域中的读段通常与假象信号相关。映射到这些区域的读段比例（与映射到峰区域的读段相比）较高的细胞通常代表技术假象，应予以去除。FractionCountsInRegion() 函数可用于计算每个细胞在给定区域集合中所有计数的比例。我们可以使用该函数和黑名单区域列表来找出每个细胞的黑名单计数比例。
请注意，后三个指标可以从 CellRanger 的输出中获取（这些输出已存储在对象元数据中），但也可以使用 Signac 为非 10x 数据集计算（更多信息请参见本文档末尾）。
# compute nucleosome signal score per cell
pbmc <- NucleosomeSignal(object = pbmc)
# compute TSS enrichment score per cell
pbmc <- TSSEnrichment(object = pbmc)
# add fraction of reads in peaks
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
# add blacklist ratio
blacklist_regions <- ah[['AH107305']] # blacklist regions for hg38
pbmc$blacklist_ratio <- FractionCountsInRegion(
object = pbmc,
assay = 'peaks',
regions = blacklist_regions
)
对象元数据中存储的变量之间的关系可以使用 DensityScatter() 函数进行可视化。通过设置 quantiles=TRUE，该函数还可以用于快速查找不同质控指标的合适筛选阈值：
DensityScatter(pbmc, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
我们还可以观察所有细胞的片段长度周期性，并根据核小体信号强度的高低对细胞进行分组。可以看到，在（基于上图）单核小体/无核小体比率上表现异常的细胞，其核小体条带模式也有所不同。其余细胞则展现出典型的成功ATAC-seq实验所应有的模式。
pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')
我们可以使用小提琴图分别绘制每个QC指标的分布情况：
VlnPlot(
object = pbmc,
features = c('nCount_peaks', 'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal', 'pct_reads_in_peaks'),
pt.size = 0.1,
ncol = 5
)
下期继续~~~~
欢迎大家交流补充~有问题请评论区或者私信留言~~~
