10X Cell Ranger ATAC 算法概述

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

ATAC-seq 技术简介 sc-ATAC-seq细胞类型注释策略

Barcode Processing

执行此步骤是为了修复条形码(barcode,细胞的标识)中偶尔出现的测序错误，从而使片段与原始条形码相关联，从而提高数据质量。16bp条形码序列是从“I2”索引读取得到的。每个条形码序列都根据正确的条形码序列的“白名单”进行检查，并计算每个白名单条形码的频率。我们试图纠正不在白名单上的条形码，方法是找出所有白名单上的条形码，它们与观察到的序列之间的2个差异(汉明距离（Hamming distance）<= 2)，并根据reads数据中条形码的丰度和不正确碱基的质量值对它们进行评分。如果在此模型中，未出现在白名单中的观察到的条形码有90%的概率是真实的条形码，则将其更正为白名单条形码。

Alignment

Cell Ranger ATAC执行基于参考（reference-based）的分析，并要求adapter 和引物寡核苷酸序列（primer oligo sequence）在确定映射之前进行修剪。在目前的策略中，如果读长度大于基因组片段的长度，读序列的3'端(读序列的末端)可能包含引物序列的反补序列。我们使用cutadapt工具在每次读取结束时识别引物序列的反向补码，并在比对之前从读取序列中对其进行修剪（trimmed ）。

然后，使用带有默认参数的BWA-MEM将修剪后的读对（read-pairs）对齐到指定的引用（reference ）。请注意BWA-MEM不会将读码小于25bp对齐（align ）。这些未对齐的读包含在BAM输出中，并标记为未映射（unmapped）。

Duplicate Marking

由于PCR扩增，一个条形码片段（fragment ）可能被测序多次。我们标记副本是为了识别构成库的原始片段（fragment ）并增加其复杂性。我们通过识别所有条形码上的一组读码对来发现重复的读码，其中R1和R2的5'端在参考上具有相同的映射位置，可以进行软裁剪校正。这些读对来自于同一个原始分子。在这些读取对中，最常见的条形码序列得到了识别。带有条形码序列的一个读对被标记为“原始的”，组中的其他读对被标记为BAM文件中该片段的副本。如果它通过了下一段描述的过滤器，这是片段文件中作为片段报告的唯一读对，并且它将被标记为最常见的条形码序列。

在处理如上所述的一组相同排列的读码对时，一旦标记了原始片段，我们将确定该片段是否在两次读码时都使用MAPQ > 30进行了映射，它不是线粒体，也不是嵌合映射。如果片段通过这些过滤器,我们创建一个条目在fragments.tsv.gz文件的开始和结束标记片段调整后5 '末端的read-pair占换位,在转座酶DNA占据了一个地区的9碱基对长(见图)。在这个条目中，我们关联了为读对组观察到的最常见的条形码，以及这个片段在库中观察到的次数(组的大小)。注意，作为这种方法的结果，基因组上的每个唯一间隔只能与一个条形码相关联。每个条目是用选项卡分隔的，文件是位置排序的，然后使用默认参数运行SAMtools tabix命令。

Peak Calling

由于每个片段的末端表示开放染色质区域，我们分析来自这些片段的组合信号，以确定为开放染色质而富集的基因组区域，从而推定具有调控和功能意义。使用位置排序的片段文件中片段末端所确定的位置，我们计算了基因组中每个碱基对上的转位事件的数量。

我们在每个基对周围用一个401bp的移动窗口和生成这些事件的平滑轮廓，并拟合一个津巴式的混合模型（ZINBA），该模型由几何分布对零充气计数、负二项分布对噪声和另一个几何分布对信号进行建模组成。拟合的方法是保证混合分量均值的排序是固定的:首先是负二项噪声均值，然后是几何零膨胀均值，最后是几何信号分布均值。我们根据1/5的赔率（odds-ratio）设置一个信号阈值，该阈值决定了在碱基对分辨率下，一个区域是峰值信号(为开放染色质而富集)还是噪声。因此，并不是所有的切割点都在一个峰值区域内。将彼此之间500bp以内的峰值合并到一起，生成一个位置排序的峰值BED文件。

窗口大小,odds-ratio和距离选择操作者曲线下的面积最大化高峰时相比,一个高质量的DNase hypersensitive sitesGM12878从ENCODE,并产生令人满意的聚类指标以及细胞类型识别一组外周血单核细胞(PBMC)库。这种识别峰的方法独立于条形码和它们的细胞(或非细胞)身份，这使我们能够包含所有由映射确定的真实基因组片段的信号。

Cell Calling

此步骤将库中观察到的条形码子集与从样本文库的cell相关联。这些细胞条形码的识别允许人们在单细胞分辨率下分析数据的变化。对于每个条形码，我们有通过所有（the fragments.tsv file)的映射高质量片段的记录。在此之前，我们已经确定了峰值，我们使用重叠于任何峰值区域的片段（fragments ）的数量，对于每个条形码，来将信号从噪声中分离出来。与使用每个条形码的片段数量相比，这在实践中效果更好。cell Calling 分两个步骤完成。首先，我们识别出有片段重叠部分的条形码，这些重叠部分称为峰值，低于基因组的峰值部分(仅为计算的目的，为了说明片段长度，峰的两边都填充了2000 bp)。我们发现，这些条形码的切割位点通常随机分布在基因组中，不以功能区域附近的富集为目标，也不表现出预期的ATAC-seq“峰值”信号。因此，我们屏蔽了这些“低目标”条形码，使其不出现在Peak Calling之前库中观察到的所有条形码。

10x Chromium 系统的凝胶珠多联率较低(主要是双联)，其中一个细胞共享一个以上的条形码凝胶珠。然后，这些cell在数据集中显示为相同单元类型的多个条形码。这些额外的条形码的存在并不影响二次分析，如聚类分析或差异分析，尽管它可能会增加对非常罕见细胞类型的丰度测量。我们通过观察这对条形码是否彼此共享更多基因组上相连的“连接”片段(共享一个移位事件的片段)(B1-B2)，而不是它们自己(B1-B1或B2-B2)，从而识别出推定凝胶珠双重态的一个小主条形码对(B1, B2)。次要条形码被标识为片段较少的条形码，并从cell calling中使用的总条形码集中丢弃。单细胞ATAC数据还有另一个来源，可以产生类似类型的额外细胞。这种现象被称为条形码多联，当一个细胞相关的凝胶珠不是单克隆的，并且存在一个以上的条形码时，就会发生这种现象。与这种多联体相关的条形码被识别为共享大量相互连接的片段以及具有公共后缀或前缀核苷酸序列的条形码。同样，我们掩盖了参与这些多联的“次要”条形码，同时保留了作为相关cell的唯一代表的主条形码。

然后，我们对剩余的条形码执行cell calling。我们从所有的条形码计数中减去与深度相关的固定计数，从而对白名单污染进行建模。这个固定的计数是来自不同GEM的每个条形码fragments 的估计数量，假设污染率为0.02。然后，我们拟合了两个负二项分布的混合模型来捕获信号和噪声。将比值比设置为100000(这似乎在内部测试中效果最好)，我们将与实际细胞相对应的条形码与非细胞条形码分开。

在参考中，cell calling 被限制在每个物种产生< 20k细胞，因为目前的实验设计支持500-10k细胞。如果 --force-cells=N作为Cell Ranger ATAC的参数提供，我们将确保按照上面的条形码等级图，将每个物种的前N个条形码报告为Cell。n> 20k将不被管道接受。如果没有提供--force-cells，在混合物种样本的情况下，我们进行第二次迭代，在这里我们掩盖了非细胞条形码，并将相同的混合模型适用于细胞条形码和ref中出现的两种物种分布.一般来说，要使用的--force-cells值应该与上面的条形码等级图中看到的“膝盖” (knee)下面的一点相对应。

Peak-Barcode Matrix

与我们对单细胞基因表达解决方案和单细胞免疫分析解决方案的分析管道类似，我们生成一个计数矩阵，其中包含每个条形码每个峰区域内的片段末端(或切割位点)计数。这是原始的峰条形码矩阵，它捕获每个条形码的开放染色质的富集。然后将矩阵过滤为只包含单元条形码，然后将其用于随后的分析，如降维、聚类和可视化。

Dimensionality Reduction, Clustering and t-SNE Projection

生物发现常常借助于可视化工具，这些工具允许一个人将一个细胞群与另一个细胞群进行分组和比较。为了实现发现，Cell Ranger ATAC执行集群和t-SNE投影。由于数据在单细胞分辨率下是稀疏的，我们首先进行降维，将其投射到更低的维度空间，这也具有去噪的优点。Cell Ranger ATAC通过主成分分析(PCA)、潜在语义分析(LSA)或概率潜在语义分析(PLSA)支持降维。所采用的默认方法是LSA，但是用户可以通过向Cell Ranger ATAC提供降维参数(—dim-reduce=< method >)来指定使用哪个方法。这些方法中的每一种都作用于经过过滤的峰条形码矩阵，该矩阵由称为峰的cell条形码的切割位点计数组成。每个方法都有一个在降维之前使用的相关数据归一化技术和一组接受降维后数据的聚类方法。我们还提供了Barnes Hut TSNE算法的优化实现(该算法与我们针对单细胞基因表达解决方案的分析管道中的算法相同)。维数固定为15，因为在外周血单核细胞(PBMCs)上进行测试时，发现它能够以视觉上和生物学上有意义的方式充分分离聚类。

PCA

对于PCA，我们首先将数据归一化为每个条形码的中间切割点计数，并对其进行log转换。我们使用了一种快速、可伸缩和内存有效的IRLBA实现(增强的、隐式重新启动的Lanczos双对角化算法)，它允许原地定心和特征缩放，并生成转换后的矩阵以及主成分(PC)和奇异值，这些奇异值对每台PC解释的方差进行编码。针对PCA，我们提供了k-means聚类，可以生成2到10个用于可视化和分析的聚类。我们还提出了一种基于社区检测的k近邻图聚类方法，该方法采用louvain模块化优化算法。转换后的矩阵由默认参数的t-SNE算法操作，并为每个条形码提供二维坐标以进行可视化。使用我们的单细胞基因表达解决方案的用户可能会认识到，使用PCA进行的分析类似于运行Cell Ranger (cellranger count)。

LSA

灵感来自于大量的信息检索领域的工作,我们通过文件规范化数据频率(idf)变换,每个峰值计数是缩放的日志的数量的比率矩阵条形码,条形码的数量的峰值有非零的数。这为出现在更少条形码中的峰值计数提供了更大的权重。利用不定标、不定心的IRLBA对该归一化矩阵进行奇异值分解(SVD)，生成低维空间的变换矩阵，以及表示各分量重要性的分量和奇异值。在聚类之前，我们通过在低维空间中将每个条形码数据点缩放到单位L2-norm来对深度进行归一化。我们发现这些标准化技术的组合避免了删除第一个PC的需要。针对LSA，我们提供了球形k-means聚类，可以产生2到10个用于下游分析的聚类。通过在l2归一化的球形流形数据上使用k-means识别簇，球形k-means的性能优于普通k-means。与PCA类似，我们还通过t-SNE提供了一个基于图的集群和可视化。但是，与球形k-means聚类相似，在进行基于图的聚类和t-SNE投影之前，我们将数据归一化为单位范数。

PLSA

PLSA是一种特殊类型的非负矩阵分解，起源于自然语言处理。在PLSA中，通过期望最大化算法，我们最小化了经验确定的条形码中观察到峰值的概率与该峰值的低秩近似之间的kl -散度。在通过PLSA降维之前，我们不会对数据进行归一化处理。与LSA和PCA类似，我们生成一个转换矩阵、组件向量和一组值来解释每个组件的重要性。PLSA提供了组件和转换矩阵的自然解释。每个组件都可以被解释为一个隐藏的主题，转换后的矩阵就是从给定主题观察到条形码的概率，即Prob(条码|主题)。分量向量是一个给定主题(Prob(peak|topic))观察到一个峰值的概率，LSA/PCA的奇异值对应的是数据中观察到的每个主题(Prob(topic))的概率。与LSA类似，我们将变换后的矩阵归一化为单位l2范数，并进行球形k-means聚类，生成2到10个聚类，并通过t-S实现基于图的聚类和可视化.虽然PLSA在低维空间的可解释性方面提供了巨大的优势，但它比PCA和LSA都要慢得多，而且在大型数据集上也不能扩展超过20个组件。为了在一定程度上改进这一点，PLSA的内部实现是多线程的(计算集群上的4个线程)，用c++编写和编译。为了确保一个合理的运行时间，如果不首先收敛，则算法的迭代次数上限为3000次。

下面总结了降维技术以及相关的聚类和可视化方法。

Peak Annotation

由于峰是富含开放染色质的区域，因此具有潜在的调节功能，因此观察峰相对于基因的位置是很有意义的。我们使用最接近-D=b的工具将每个峰与基于最接近转录起始位点的基因联系起来(包装在参考文献中)，这样，峰在TSS上游的1000个碱基或下游的100个碱基内。此外，我们还将基因与假定的远端峰值相关联，这些远端峰值距离TSS远得多，且位于转录本末端上游或下游的距离小于100kb。这种关联被我们的可视化软件(Loupe Cell Browser)采用，用于构建和可视化派生的特性，比如启动子和，它们从与基因相关的峰值累积计数。转录因子(TF)结合基序的峰被富集，某些基序的存在表明了转录因子的活性。为了识别这些基序，我们首先计算峰值的GC%分布，然后将这些峰值分配到GC内容分布中的相等分位数范围。我们使用包装在Cell Ranger ATAC中的Python库来扫描每个峰，寻找与motif位置权重矩阵(motif position-weight-matrices, PWMs)匹配的转录因子，这些转录因子来自直接构建在参考包中的JASPAR数据库。我们将p值阈值设置为1E-7，背景核苷酸频率设置为每个GC桶中峰值区域内观察到的核苷酸频率。在这些bucket上统一了motif-peak匹配列表，从而避免了扫描过程中的GC偏差。

TF Motif Enrichment Analysis

由于转录因子(TF)倾向于在包含其同源基序的位点上结合，因此将可达性测量值与常见基序进行分组，可以在单个细胞间对TFs进行有益的富集分析。我们按照以下方式为每个细胞条码的每个TF构建一个整数计数:我们考虑所有与给定TF匹配的峰值，正如在TF motif检测中发现的那样。然后，对于每个条形码，我们将经过筛选的峰值条形码矩阵中这些匹配峰的切割点计数汇集在一起。它计算细胞条形码中共享TF基序的峰值的总切割点。我们计算了一个条形码中TF的切割位点占该条形码所有切割位点的比例，从而使其归一化到深度。我们通过对给定TF的这些比例值在条形码上的分布进行z分数来检测TF的富集。为了使它对异常值具有鲁棒性，我们使用修改后的z分数，该分数使用中位数和中位数与中位数的比例绝对偏差(MAD)计算，而不是使用平均值和标准偏差。当您在Loupe中加载数据集并将差异分析内置于Loupe中时，这些z分数值是可见的。

Differential Accessibility Analysis

为了识别每个簇可达性不同的转录因子motifs, Cell Ranger ATAC对每个motif和每个cluster进行测试，看簇内均值和簇外均值是否存在差异。熟悉我们的单细胞基因表达解决方案的用户可能会认识到，这与Cell Ranger用于识别差异基因表达的工作是相同的。为了发现细胞群之间的差异可达基序，Cell Ranger ATAC使用负二项(NB2)广义线性模型来发现集群的特定均值及其标准差，然后使用Wald检验进行推理。对于每个集群，相对于所有其他cell，该算法在该集群上运行，生成一个TF基序列表，这些TF基序在该集群中相对于样本的其余部分有差异的表达。使用GLM框架，我们可以将每个细胞的测序深度和每个细胞的GC含量峰值直接作为协变量进行建模。这使我们能够将它们自然地归一化，作为模型估计和推断过程的一部分。我们还使用泊松广义线性模型对峰中的可达性进行微分富集分析，这与我们对TF基序的分析方法非常相似。在这种情况下，我们只将每个单元深度作为协变量进行建模。

Aggregation

在aggr管道中，用户可以提供要聚合的库列表。聚合管道根据运行时指定的规范化模式，执行将每个列出的库中的片段合并到一个聚合文件中的任务。合并是通过向下采样每个库来执行的，速率由标准化模式决定。如果规范化模式为“None”，则保留所有片段并合并在一起。如果归一化模式是“深度”，则每个库都向下采样以具有相同的灵敏度(定义为每个单元格片段的中位数)。如果归一化模式是“信号”，则下采样率是利用每个文库中沿基因组分布的剪切位点的信息来确定的。具体来说，对于每个库，我们构建了一个窗口分割站点计数的分布，并拟合了3个组件的混合模型，这与我们在峰值调用中所做的工作是相同的。下采样率是通过匹配每个库的信号分量的平均值来设置的。一旦这些fragments 合并在一起，它们就按照位置进行排序，并被制成表格以供后续使用，如降维、聚类、可视化和差异分析。

References

[1] ATAC-seq 技术简介: https://www.jianshu.com/p/45be11b879e3 [2] ZINBA: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3218829/ [3] DNase hypersensitive sites: https://links.jianshu.com/go?to=http%3A%2F%2Fhgdownload.cse.ucsc.edu%2Fgoldenpath%2Fhg19%2FencodeDCC%2FwgEncodeAwgDnaseUniform%2FwgEncodeAwgDnaseUwdukeGm12878UniPk.narrowPeak.gz [4] single-cell-atac-algorithms-overview: https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/algorithms/overview