转坐酶（Transpose）检测染色质开放性（Accessibility）结合高通量测序（Assay for Transposase-Accessible Chromatin with high-throughput sequencing，ATAC-Seq）技术能够一次性检测样本中染色质的开放性，无需像ChIP-Seq技术一样需要选择某一个TF进行试验。另外，由于测序成本远远低于它的同行——DNase-Seq、FAIRE-Seq、DNase-Seq和MNase-Seq，目前ATAC-Seq几乎终结了染色质开放性领域的竞争。ATAC-Seq的技术几乎完全移植到了单细胞领域，实验包括三个过程：细胞裂解、转座和扩增。目前scATAC-Seq主要存在两类方法——Split-and-Pool方法和微流体方法。前者的代表是美国华盛顿大学的Cole Trapnell实验室的sci-ATAC-Seq技术；后者的代表是Integrated Fluidics Circuit（IFC）和已经商业化的10X Genomics Chromium Single Cell ATAC。

我们先分别简要介绍下两类方法的原理。sci-ATAC-Seq是先后用多个具有编码功能的Plate来标记核质所属的细胞ID，过程如下：首先，把裂解后的细胞核质均匀分配到第一个Plate（包含96个Well）中，每个Well拥有一个唯一的Barcode，用来标记这个Well下的转座酶。也就是说，一个Well下所有的转座酶具有同样的Barcode。然而，96个Barcode显然不能标记上千个细胞，因此我们在后面还需要第二个Barcode。然后，这96个Well中的核质充分混合后，我们用Fluorescence-Activated Cell Sorting（FCS）技术对细胞进行排序，均匀分配到第二个Plate中。第二个Plate也是由96个Well组成。随后，我们用聚合酶链式反应（Polymerase Chain Reaction，PCR）在第二个Plate中对核质进行扩增。扩增的引物（Primer）上含有第二个Barcode，与第二个Plate中的Well一一对应。两个Barcode能形成96 * 96 = 9216个组合，用来作为上千个细胞的ID。根据经验，sci-ATAC-Seq能够检测1500左右的细胞，其中每个细胞中所包含的Read含量中位数约为2500，ID中有11%的冲突率，也就是两个细胞共用了同一个ID。如果想蹭更大规模的细胞样本，我们只需要增加Plate的数量。sci-ATAC-Seq技术的缺点有二：1）随着Plate数量的增加，测序成本大大增加；2）sci-ATAC-Seq没有实现商业化，因此在实际操作中缺乏客观有效的操作流程，往往需要针对不同的批次（Batch）进行定制和微调。我们以10X Genomics Chromium Single Cell ATAC为例，介绍基于微流体方法的技术。首先细胞核质和标有Barcode的凝胶珠（Barcoded Gel Beads）分别保存在一个容器中。然后，两者分别从两条管道以垂直角度往关岛的十字交叉口运动。这样，有的凝胶珠会粘住一个核质，然后它们会被一个油滴包裹住，凝胶珠上的Barcode会与核质结合并作为其ID。然后，已经拥有Barcode的核质通过PCR进行扩增。最后，清洗油滴并且测序。详细的实验流程可以参考此文献。

以上是scATAC-Seq测序技术的介绍。与scRNA-Seq一样，scATAC-Seq的测序数据需要大量的预处理才能进行下游分析。关于scATAC-Seq的预处理和下游分析，这篇文献介绍得最详尽。scATAC-Seq的预处理包括六步——读段（Read）预处理、质量控制（Quality Control，QC）、峰（Peak）注释、矩阵构建、矫正批次效应（Batch Effect Removal）以及降维（Dimensionality Reduction）/可视化/聚类。

一、Read预处理。这一步包括三小步——Demultiplexing、Adaptor Trimming和比对（Alignment）。1）为了节省成本，我们往往将多个样本（Sample）混合同时测序，并用Index Adaptor Sequence标记出样本。因此，Demultiplexing要做的就是把样本的ID和细胞的ID（Barcode）拷贝到这个细胞下的每个Read中。我们通常用Illumina的bcl2fastq进行Demultiplexing。2）Adapter Trimming的目的就是把测序中添加的Adaptor Sequence和Primer序列切除掉，工具是AdapterRemoval和Trimmomatic。3）Alignment的作用是把测序Read比对到参考基因组上，常用工具是Bowtie2、BWA和STAR。并不是所有Read都能比对到基因组上，比对成功的一对Read称为片段（Fragment）。sci-ATAC-Seq数据预处理后得到的是二进制BAM文件（参考文献），10X Genomics Chromium Single Cell ATAC数据预处理后的得到的是Fragments文本文档。10X Genomics Chromium Single Cell ATAC配套有“御用”预处理工具箱CellRanger，可以完成所有Read的预处理步骤。

二、QC。这一步最主流的工具是Seurat工具箱下的Signac软件，目标是通过五个标准删减质量低的细胞。1）Nucleosome Banding Pattern：染色质是缠绕在组蛋白上形成核小体（Nucleosome）的。核小体的体积是固定的，因此转座酶切割染色质的时候是避开核小体区域的，所以核小体附近切割下来的Fragment长度（Insert Size）是200 bp的倍数，也就是200，400和600 bp。Signac能够根据这些模式识别核小体附近的Fragment和染色质开放区域的Fragment。对每一个细胞，Signac计算核小体附近Fragment的数量与开放区域Fragment的数量的比值。这个比值越小越好。2）转录起始位点丰度得分（Transcription Start Site Enrichment Score，TSS Enrichment Score）：ENCODE数据库提供了一个得分，计算TSS附近区域Fragment的数量和TSS侧翼区域Fragment的数量的比值，作为TSS丰度得分。通常，测序质量低的ATAC Fragment的TSS丰度得分会很低。因此，TSS丰度得分越高越好。3/4）Peak中的Fragment数量/比例：这一步在完成Peak注释后才能做。类似地，我们也可以计算Peak中的Read数量/比例。比例越大越好，数量适中即可。如果Fragment或Read的数量太低，表示这个细胞中的Read没有被充分测得，是低质量细胞；如果Fragment或Read的数量太高，表示我们可能把多个细胞核质错当作一个核质对待了，这种情况也要去除。5）Blacklist区域：ENCODE数据库提供了Blacklist区域列表（人类、小鼠、果蝇和线虫）。这些区域倾向于被大量的Read覆盖，是需要被去除的技术假阳性。这些区域也被移植到了Signac工具箱中。

三、Peak注释。Peak注释不等于Peak Calling，而是包含Peak Calling。这一步是为第四步的矩阵构建做准备。Peak注释包括以下六种：1）TF模体（Motif）：也就是转录因子在Peak上的结合位点。斯坦福大学的William Greenleaf开发的chromVAR工具就做了这种注释。2）k-mer：即4^k种ACGT的字符串。chromVAR中也采用了这种注释。3）TSS：基因的转录起始位点承载了较多的信息，因此有些研究顺式（cis-）转录调控的工具（例如Trapnell的Cicero）会采用这种注释。4）Bin/Window：把基因组切割成不重叠的固定长度的区间（通常是5 kb），在下游分析中做降维或者聚类。UCSD的任兵开发的snapATAC和Greenleaf开发的ArchR都采用了这种注释，只是Bin的长度不同，ArchR用的长度是500 bp。这种小的Bin能够更精确地检测到TF的结合位点的位置范围（300-500 bp）。5）Peak：这是最主流、最复杂的注释方式，需要Peak Calling工具（Peak Caller）完成。多数软件（cisTopic、snapATAC、SCALE、scATAC-pro、Signac和ArchR）都采用了这种注释。6）基因：与TSS类似，区别就在于TSS只考虑聚集在TSS附近的Fragment，而基因则考虑整个基因的区域（尤其是编码区）检测到的Fragment。7）Topic：类似于双聚类（Bicluster），与双聚类的区别就在于它描述的是与Peak或者细胞的概率上的关系，而不是硬性的0-1关系。目前只有cisTopic采用这种注释。1-4）和6-7）等方法比较简单。但是，由于单细胞中Read数量太少，ChIP-Seq用的Peak Caller无法检测到Peak。因此Peak Calling都是在Bulk规模上进行，也就是大量细胞的混合物。目前存在五种Peak Calling方法：1）直接使用数据库中的DNase-Seq和ChIP-Seq的Peak，例如：cisTopic；2）使用整套数据集上的所有Read进行Peak Calling，例如：CellRanger和Cicero；3）使用一个细胞系（Cell Line）上所有的Read进行Peak Calling，例如：chromVAR；4）使用一个细胞类型（Cell Type）下的所有Read，例如：snapATAC；5）两阶段方法，即先用其他注释（例如：Bin）得到Feature-Cell矩阵，并作聚类，然后把每一个类中所有的Read汇总起来做Peak Calling，例如：scATAC-pro和ArchR。需要强调的是，在目前的单细胞多组学（RNA+ATAC）数据中，也存在三种Peak Calling方法：1）用所有细胞的Read；2）用一个细胞系下的Read；3）先用整合算法对RNA+ATAC数据进行聚类，或者只对scRNA-Seq进行聚类，然后在每个类中分别作Peak Calling。

四、矩阵构建。针对不同的注释，矩阵元素有不同的计算方法。例如：对TF Motif或者k-mer，我们首先要得到Peak，然后在Peak上搜索TF的Motif和k-mer。计算方法是先对TF Motif/k-mer所处的Peak的开放性数值相加，然后计算这个值在所有细胞中的z值。对基因/TSS，我们需要计算基因活性得分。Cicero采用的是线性模型，把与这个基因具有顺式调控关系的Peak的开放性数值叠加。然而，scATAC-Seq矩阵具有很高的稀疏性，平均每个细胞只能覆盖1～10%的Peak。考虑到scATAC-Seq是二进制，因此矩阵的信息量很低。为了提高信息量，我们需要把scATAC-Seq矩阵的0-1值根据Peak的唯一性转化为浮点数，使得越罕见的Peak，分数越高。目前存在三种方法：1）文本挖掘算法Term-frequency inverse-document-frequency（TF-IDF）：这是最主流的算法，它会给罕见的Peak更高的分值，构造出新矩阵。转化后的矩阵有利于识别不同细胞类型之间高度可变的Peak。2）Jaccard指数：它通过计算两个细胞之间共有的Peak来发现某一个细胞所特有的Peak。3）测序深度：另外一类算法不用0-1值构造矩阵，而是用一个细胞内的测序深度作为矩阵的值。

五、矫正批次效应。矫正批次效应的问题等价于Intra-Modality的多数据整合。目前的算法大多基于这样一个假设：不同的Batch之间至少共享一个细胞类型；另外，Batch之间的差异要小于细胞类型之间的差异。这些算法在矫正批次效应的过程中，有时候会把一些生物本身的特征消除，从而导致过度矫正。与scRNA-Seq不同，目前还不存在整合scATAC-Seq数据的算法。一项Benchmark的研究比较了不同的scRNA-Seq整合工具在scATAC-Seq数据上的性能。结果显示，大多数工具的表现都不尽如人意。只有Harmony、Seurat v3和scVI相对来说表现较为突出。

六、降维/可视化/聚类。scATAC-Seq具有比scRNA-Seq更高的稀疏性和维度。维度越高，我们越难以对细胞进行分类。举个例子，我们可以把测序数据理解为从全体细胞中随机均匀抽取的一个样本。如果只有一个Peak，也就是说数据是一维的（一个0-1之间的浮点数就是特征的所有信息），如果要保证每隔一段距离（0.01）就选一个细胞，我们只需要选取100个细胞。如果有两个Peak（一对0-1之间的浮点数就是特征的全部信息），全集就是一个边长为1的正方形。如果保证每隔0.01就随机选取一个细胞，我们需要100² = 10000个细胞。以此类推，维度越高，样本就越难以覆盖全集。因此，全集中某些区域可能一个细胞都没被抽取。所以，降维的性能严重影响下游的分析。目前scATAC_Seq数据的降维算法有五种：1）主成分分析（Principal Component Analysis，PCA）：是一种线性降维算法，计算速度快，但是难以反映数据内部的非线性关系。BROCKMAN、SnapATAC和Cusanovich2018都采用了PCA。但是只用PCA会造成大量细胞之间具有很高的相似性，因为每个细胞的剖面（Profile）中都含有大量的零值。因此，PCA通常与非线性降维算法一起使用。2）Topic：它是基于Latent Dirichlet Allocation（LDA）的降维算法。它运算速度很慢，但是能够识别出具有细胞类型特异性的特征，能够显著提高其他下游分析（例如：聚类）的准确性。cisTopic采用了该算法。3）Latent Semantic Indexing（LSI）：它结合了TF-IDF和奇异值分解（Singular Value Decomposition，SVD）。Seurat v4、Signac、ArchR、BROCKMAN、Cusanovich2018和Signac都采用了这种方法。需要注意的是，LSI的第一个维度会高度地与测序深度相关，因此在下游分析中我们会丢弃第一个维度。4）Multimensional Scaling（MDS）：原理很简单，它通过描述细胞之间的剖面的相似性完成降维。Scasat采用了MDS算法。5）Diffusion Map：是一种非线性降维算法。它对噪音具有较高的鲁棒性（Robustness）。snapATAC采用了该算法。通常，先用线性降维算法然后在其基础上使用非线性降维，会有更好的聚类性能。目前最流行的非线性降维算法是t-分布随机邻近嵌入（t-distributed stochastic neighbor embedding，t-SNE)和统一流形逼近与投影(Uniform Manifold Approximation and Projection，UMAP)。t-SNE适合发觉细胞之间的局部邻近信息，UMAP擅长发现全局邻近信息。但是t-SNE和UMAP只能用前两个维度做可视化，不适合把降维信息用于下游分析。Cicero使用t-SNE或UMAP的降维信息构造k最邻近图（k Nearest Neighbor Graph），实际上是不恰当的。关于聚类，这项Benchmark的研究的结果显示，Louvain算法具有最好的性能；k-medoids算法具有最高的鲁棒性。

以上便是关于scATAC-Seq的预处理分析和主要工具。实际上，我们一般选用一个集成化的工具箱完成所有分析，从而避免选择困难。目前功能最全的scATAC-Seq分析工具有两个：哈佛大学Shirley Liu团队的MAESTRO和Greenleaf的ArchR。尤其是ArchR，它能够完成：Read预处理、QC、Peak Calling、矩阵构建、降维、可视化、聚类、足迹（Footprinting）分析、共开放性（Co-accessibility）分析、轨迹（Trajectory）分析、整合分析和连接Peak和基因，等等。注意：ArchR除了采用MACS2对Fragments坐Peak Calling之外，它还具备一个基于Tile矩阵（即Bin = 500 bp的Bin-Cell矩阵）新的Peak Calling算法。两者性能类似，但是作者更推荐用MACS2。关于足迹分析，ArchR能构造出每个TF所结合的位点附近的Fragment覆盖情况。

上述就是scATAC-Seq的所有预处理和初步分析的介绍。其实在我用Seurat v4期间，它对scATAC-Seq数据的处理流程就发生过多次变化。可以预见，未来的预处理会包含更多的步骤、更全的角度以及更快的速度。