单细胞笔记9-scATAC-seq测序方法和数据分析

作者: 江湾青年 | 来源:发表于2021-09-16 20:43 被阅读0次

引言

转坐酶（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。

1. sci-ATAC-Seq

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）进行定制和微调。

2. 基于微流体方法的scATAC-seq

我们以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）。

为了节省成本，我们往往将多个样本（Sample）混合同时测序，并用Index Adaptor Sequence标记出样本。因此，Demultiplexing要做的就是把样本的ID和细胞的ID（Barcode）拷贝到这个细胞下的每个Read中。我们通常用Illumina的bcl2fastq进行Demultiplexing。
Adapter Trimming的目的就是把测序中添加的Adaptor Sequence和Primer序列切除掉，工具是AdapterRemoval和Trimmomatic。
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注释包括以下7种：

TF模体（Motif）：也就是转录因子在Peak上的结合位点。斯坦福大学的William Greenleaf开发的chromVAR工具就做了这种注释。
k-mer：即4k种ACGT的字符串。chromVAR中也采用了这种注释。
TSS：基因的转录起始位点承载了较多的信息，因此有些研究顺式（cis-）转录调控的工具（例如Trapnell的Cicero）会采用这种注释。
Bin/Window：把基因组切割成不重叠的固定长度的区间（通常是5 kb），在下游分析中做降维或者聚类。UCSD的任兵开发的snapATAC和Greenleaf开发的ArchR都采用了这种注释，只是Bin的长度不同，ArchR用的长度是500 bp。这种小的Bin能够更精确地检测到TF的结合位点的位置范围（300-500 bp）。
Peak：这是最主流、最复杂的注释方式，需要Peak Calling工具（Peak Caller）完成。多数软件（cisTopic、snapATAC、SCALE、scATAC-pro、Signac和ArchR）都采用了这种注释。
基因：与TSS类似，区别就在于TSS只考虑聚集在TSS附近的Fragment，而基因则考虑整个基因的区域（尤其是编码区）检测到的Fragment。
Topic：类似于双聚类（Bicluster），与双聚类的区别就在于它描述的是与Peak或者细胞的概率上的关系，而不是硬性的0-1关系。目前只有cisTopic采用这种注释。
1-4）和6-7）等方法比较简单。但是，由于单细胞中Read数量太少，ChIP-Seq用的Peak Caller无法检测到Peak。因此Peak Calling都是在Bulk规模上进行，也就是大量细胞的混合物。

目前存在五种Peak Calling方法：

直接使用数据库中的DNase-Seq和ChIP-Seq的Peak，例如：cisTopic；
使用整套数据集上的所有Read进行Peak Calling，例如：CellRanger和Cicero；
使用一个细胞系（Cell Line）上所有的Read进行Peak Calling，例如：chromVAR；
使用一个细胞类型（Cell Type）下的所有Read，例如：snapATAC；
两阶段方法，即先用其他注释（例如：Bin）得到Feature-Cell矩阵，并作聚类，然后把每一个类中所有的Read汇总起来做Peak Calling，例如：scATAC-pro和ArchR。

需要强调的是，在目前的单细胞多组学（RNA+ATAC）数据中，也存在三种Peak Calling方法：

用所有细胞的Read；
用一个细胞系下的Read；
先用整合算法对RNA+ATAC数据进行聚类，或者只对scRNA-Seq进行聚类，然后在每个类中分别作Peak Calling。

五、矩阵构建

针对不同的注释，矩阵元素有不同的计算方法。

对TF Motif或者k-mer，我们首先要得到Peak，然后在Peak上搜索TF的Motif和k-mer。计算方法是先对TF Motif/k-mer所处的Peak的开放性数值相加，然后计算这个值在所有细胞中的z值。
对基因/TSS，我们需要计算基因活性得分（Gene Activity Score）。Cicero采用的是线性模型，把与这个基因具有顺式调控关系的Peak的开放性数值叠加。然而，scATAC-Seq矩阵具有很高的稀疏性，平均每个细胞只能覆盖1～10%的Peak。考虑到scATAC-Seq是二进制，因此矩阵的信息量很低。为了提高信息量，我们需要把scATAC-Seq矩阵的0-1值根据Peak的唯一性转化为浮点数，使得越罕见的Peak，分数越高。
目前存在三种方法：
- 文本挖掘算法Term-frequency inverse-document-frequency（TF-IDF）：这是最主流的算法，它会给罕见的Peak更高的分值，构造出新矩阵。转化后的矩阵有利于识别不同细胞类型之间高度可变的Peak。
- Jaccard指数：它通过计算两个细胞之间共有的Peak来发现某一个细胞所特有的Peak。
- 测序深度：另外一类算法不用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就随机选取一个细胞，我们需要1002 = 10000个细胞。以此类推，维度越高，样本就越难以覆盖全集。因此，全集中某些区域可能一个细胞都没被抽取。所以，降维的性能严重影响下游的分析。

目前scATAC_Seq数据的降维算法有五种：

主成分分析（Principal Component Analysis，PCA）：是一种线性降维算法，计算速度快，但是难以反映数据内部的非线性关系。BROCKMAN、SnapATAC和Cusanovich2018都采用了PCA。但是只用PCA会造成大量细胞之间具有很高的相似性，因为每个细胞的剖面（Profile）中都含有大量的零值。因此，PCA通常与非线性降维算法一起使用。
Topic：它是基于Latent Dirichlet Allocation（LDA）的降维算法。它运算速度很慢，但是能够识别出具有细胞类型特异性的特征，能够显著提高其他下游分析（例如：聚类）的准确性。cisTopic采用了该算法。
Latent Semantic Indexing（LSI）：它结合了TF-IDF和奇异值分解（Singular Value Decomposition，SVD）。Seurat v4、Signac、ArchR、BROCKMAN、Cusanovich2018和Signac都采用了这种方法。需要注意的是，LSI的第一个维度会高度地与测序深度相关，因此在下游分析中我们会丢弃第一个维度。
Multimensional Scaling（MDS）：原理很简单，它通过描述细胞之间的剖面的相似性完成降维。Scasat采用了MDS算法。
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算法具有最高的鲁棒性。

八、细胞类型注释

尽管有许多工具可以对scRNA-seq数据自动进行细胞类型注释，还可以从各种数据库中获得细胞marker基因列表，但对于 scATAC-seq 数据，仅有有限的工具和特定细胞类型染色质可及性的参考数据集。因此，对于 scATAC-seq 数据，必须结合使用补充方法进行细胞群注释。

目前有两种方法进行细胞类型注释：

1. 基于ATAC peak的特征注释

细胞聚类后，每个细胞群的差异可及性区域可能包含不同的调控元件。细胞身份注释的第一种方法使用细胞群特异性的 peak 进行注释，监督或手动注释细胞群身份需要参考数据库或有关细胞类型特定基因组特征（例如TF motif，增强子，启动子和TSS）的文献。基于细胞类型特异的基因列表，启动子和 TSS 被最广泛地用于细胞群注释。一些简易的方法通过启动子或 TSS 上游一定距离内 peak 的存在来定义细胞类型特异性基因的可及性，而高级的分析则考虑了远端和近端调控因子的影响。“基因活性分数”对与基因启动子区共开放元件给予不同权重，从而可以更准确地利用染色质可及性推断基因表达水平。与简单的使用启动子区可及性相比，基因活性分数能更好的表征基因表达。Garnett 软件利用基因活性分数和已知细胞类型的先验特征及标记基因对细胞类型进行监督分类。

2. 利用与参考 scRNA-seq 数据的整合进行注释

这种方法将来自 scRNA-seq 数据的基因表达矩阵与来自相同细胞类型的 scATAC-seq 数据的基因活性矩阵整合在一起。将它们投影到最大相关维度后，使用 MNN 算法将细胞标记从 scRNA-seq 数据转移到 scATAC-seq 数据。尽管具有高度主导的细胞类型或与其他组学数据不匹配的细胞类型的样本显示出准确性方面的局限性，但细胞身份注释的总体结果与匹配的数据集一致。通过对 scATAC-seq 数据中的细胞群体进行半监督识别，现有的参考 scRNA-seq 和 bulk ATAC-seq 数据可用于生成 scATAC-seq 样本的网络，进而将标签进行转移。

以上便是关于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覆盖情况。

参考

http://blog.sciencenet.cn/home.php?mod=space&uid=3447504&do=blog&id=1301515