引言
在本教学指南中，我们将探讨由10x Genomics公司提供的人类外周血单核细胞（PBMCs）的单细胞ATAC-seq数据集。本节中所涉及的所有文件均可在10x Genomics的官方网站上找到：
- 原始数据文件
- 相关的元数据文件
- 包含数据片段的文件
- 用于索引的片段文件
要下载所有必需的文件，您可以在 shell 中运行以下行：
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi
加载包
首先加载 Signac、Seurat 和我们将用于分析人类数据的其他一些包。
if (!requireNamespace("EnsDb.Hsapiens.v75", quietly = TRUE))
BiocManager::install("EnsDb.Hsapiens.v75")
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)
library(ggplot2)
library(patchwork)
计算 QC 指标
现在，我们能够对单细胞ATAC-seq实验进行一系列质量控制指标的计算。我们推荐以下几种指标来评估您的数据是否达到标准。与单细胞RNA测序技术类似，这些指标的理想数值范围会根据您的生物样本、细胞的生存状态以及其他一些因素而有所不同。
核小体条带模式：通过配对末端测序得到的DNA片段大小分布图应该清晰显示出核小体条带模式，这反映了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项目已经列出了一些黑名单区域，这些区域的读数常常与技术伪影相关。如果细胞有较高比例的读数定位在这些区域（与定位在峰区的读数相比较），那么这些细胞很可能代表技术伪影，应该被排除。人类（hg19和GRCh38）、小鼠（mm10）、果蝇（dm3）以及秀丽隐杆线虫（ce10）的ENCODE黑名单区域已包含在Signac软件包中。
最后请注意，上述提到的最后三个指标不仅可以从CellRanger的输出结果中获得（这些结果会存储在对象的元数据中），而且也可以使用Signac软件包来计算非10x科技公司的数据集（有关更多信息，请参见本文档的末尾部分）。
# compute nucleosome signal score per cell
pbmc <- NucleosomeSignal(object = pbmc)
# compute TSS enrichment score per cell
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)
# add blacklist ratio and fraction of reads in peaks
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
我们可以使用DensityScatter()函数来直观展示对象元数据中存储的变量之间的相互关系。此外，通过将quantiles参数设置为TRUE，这个函数还能够帮助我们迅速确定不同质量控制指标的适宜阈值。
DensityScatter(pbmc, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
我们可以通过将细胞按照TSS富集分数进行分类，并在所有转录起始位点（TSS）上绘制其可访问性信号来审视这些分数。在调用TSSEnrichment()函数时，如果设置fast=TRUE选项，将只计算TSS富集分数，而不会存储每个细胞的Tn5插入频率的完整细胞位置矩阵，这样做可以节省内存空间。但是，请注意，如果设置了fast=TRUE，我们将无法使用TSSPlot()函数对不同组的细胞进行TSS富集信号的后续绘图分析，这一点在下面的示例中有所展示。
pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 3, 'High', 'Low')
TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()
我们还可以通过观察所有细胞的DNA片段长度的周期性变化，并根据细胞的核小体信号强度高低进行分组，来进一步分析。可以观察到，那些在单核小体与无核小体比例上表现异常的细胞（如前述图表所示）呈现出不同的核小体条带模式。而其他细胞则展现出一种符合成功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
)
最后，我们删除这些 QC 指标的异常值细胞。
pbmc <- subset(
x = pbmc,
subset = nCount_peaks > 3000 &
nCount_peaks < 30000 &
pct_reads_in_peaks > 15 &
blacklist_ratio < 0.05 &
nucleosome_signal < 4 &
TSS.enrichment > 3
)
pbmc
## An object of class Seurat
## 87561 features across 7307 samples within 1 assay
## Active assay: peaks (87561 features, 0 variable features)
## 2 layers present: counts, data
