这是目前为止我看到过的关于ATAC-seq的最新综述，感兴趣的话值得一读。
2020.4.22更新：我看到一个公众号也翻译了这篇文章，而且正好跟我这篇翻译是在同一天发表https://mp.weixin.qq.com/s/7JAEPDuEEsmRxXI3UZDZHQ。他的文字比我流畅很多，排版也比较舒服，推荐给大家~

基本信息

From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis

• 2020年2月3日发表在Genome Biology，一作为Feng Yan，通讯作者是Nicholas C. Wong。

• 原文链接：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-1929-3#Abs1

image

摘要

转座酶可及染色质测序法（ATAC-seq）已广泛用于研究染色质生物学，但分析工具的全面性综述尚未完成。在这里，我们讨论ATAC-seq数据分析的主要步骤，包括预分析（质量检查和比对），核心分析（peak calling）和高级分析（peak差异分析和注释，motif富集，footprint分析，以及核小体定位分析）。我们还回顾了利用多组学数据重建转录调控网络的过程，并重点指出了每个步骤当前面临的挑战。最后，我们描述了单细胞ATAC-seq的潜力，并强调了开发ATAC-seq特定分析工具以获得有生物学意义的深入理解的必要性。

介绍

哺乳动物的DNA通过三个主要的层次尺度进行高度浓缩：第一层次是核小体，然后包装到染色质，再通向第三层次——染色体[ 1，2，3，4，5，6 ]。染色质可以在转录活跃的常染色质和不活跃的异染色质之间进行动态切换[7，8 ]。DNA压缩的三个尺度及其相互作用共同造就了基因的表达调控。

最近的基因调控研究集中在表观遗传学上。高通量测序技术的进步给我们提供了各种破译表观遗传学图谱的方法。其中包括测定染色质可及性的转座可接近的染色质测序（ATAC-seq）[9，10 ]，DNA酶I高敏位点测序（DNase-seq）[ 11，12，13 ]和甲醛辅助隔离调控元件测序（FAIRE-seq） [ 14 ]；其测量转录因子（TF）结合[ 15，16，17]和组蛋白修饰[ 18，19 ]的染色质免疫沉淀测序（ChIP-seq）; 检测核小体定位和占位[ 20，21 ]的微球菌核酸酶测序（MNase-seq）。这些测定的详细步骤不在本综述的范围之内，在其他文章中[ 22 ]进行了详细讨论。

自2013年发明以来，ATAC-seq在各种染色质可及性的检测方法中特别受欢迎。经过整理的ATAC-seq数据集和出版物呈指数增长，表明其在广泛的生物学问题中的价值（图1a），例如如描绘哺乳动物健康组织和细胞类型中的增强子图谱[ 23，24，25 ]，研究正常造血和白血病之间的可及性变化[26，27]，以及精神分裂症患者和癌症基因组图谱（TCGA）泛癌队列中的染色质状态[ 28，29]。图3a展示了这项尖端技术在基础和转化研究中的示意图 。简而言之，ATAC-seq整合了基因工程修饰的高活性Tn5转座酶，可以同时切割开放的染色质、留下9 bp的交错缺口，并将高通量测序接头连接到这些区域。在此过程中，切口被修复，留下了9-bp的重复序列[ 30，31 ]。然后进行双端测序以使这些开放区域有更高的非重复比对率[ 32 ]。

图1 ATAC-seq数据集增长，以及预分析和高级分析的样本数总览。a* 从2013年1月1日至2019年10月1日，PubMed中ATAC-seq数据集、ATAC-seq出版物、DNase-seq数据集、FAIRE-seq的数据集、MNase-seq的数据集在的数量 b典型片段大小分布曲线显示100bp和200 bp附近的富集，表明无核小体结合和单核小体结合的片段。c典型的TSS富集图显示，不含核小体的片段在TSS富集，而单核小体的片段在TSS处耗尽，但在侧翼区域富集。d典型的峰注释饼图显示，超过一半的峰落入增强子区域（远端基因间区和内含子区域），只有约25％的峰在启动子区域。TSS：转录起始位点*

Tn5转座酶的高活性使ATAC-seq protocol成为一种简单、省时的方法，需要500-50,000个细胞[ 9 ]。灵敏度和特异性与DNase-seq相当，但优于FAIRE-seq，这两种方法都需要数百万个细胞作为输入材料[ 9 ]。由于ATAC-seq在文库制备过程中不涉及严格的大小选择，因此它也可以使用代表核小体单体和多聚体的片段来鉴定核小体位置[ 9 ]。最近，单细胞ATAC-seq（scATAC-seq）已被报道，依赖的方法有流式细胞分选（FACS）、微流体、基于纳米孔等不同类型[ 33，34，35]。scATAC-seq可以在多种情况下（包括临床标本和发育生物学等）被应用于单细胞分辨率水平研究异质性的细胞群[23，29 ]。

尽管ATAC-seq简单且鲁棒，但它存在一个主要的障碍——专门为ATAC-seq数据开发的生物信息学分析工具很少[32，36 ]。ChIP-seq和DNase-seq中使用的分析工具已应用于ATAC-seq [ 37 ]，基于它们数据特征相似的假设。但是，此假设尚未得到系统地评估。

这篇综述的主要重点是讨论ATAC-seq分析的现有资源。我们旨在为ATAC-seq数据分析提供带注释的指南，而不是详尽的工具集。此前关于ATAC-seq数据分析的综述都集中在peak calling和调控网络建模[ 37，38 ]，但现在我们迫切需要一篇涵盖ATAC-seq数据分析各个主要部分的系统性综述。这篇综述将涵盖流程图（图 2）中列出的四个最重要的步骤。其中包括（1）预分析（质量控制（QC）和比对），（2）核心分析（peak calling），（3）peak，motif，核小体和TF footprint水平的高级分析，以及（4） 与多组学数据整合以重建调控网络。这些步骤将使研究人员能够对ATAC-seq数据进行鲁棒的分析，并产生更具生物学意义的结果。最后，我们将介绍ATAC-seq分析和scATAC-seq的挑战和机遇。

图2 经典ATAC-seq分析的路线图。列出了四个主要步骤，包括预分析、核心分析、高级分析以及与多组学数据的集成。预分析包括比对前质量控制、比对和比对后处理以及质量控制。核心分析包括peak calling。高级分析包括peak，motif，footprint和核小体分析。多组学数据集成包括与ChIP-seq和RNA-seq数据整合以及调控网络重建。每个框中的文本强调每个分析步骤中的重要注意事项。我们建议研究人员用FastQC，trimmomatic和BWA-MEM进行预分析，用MACS2进行peak calling，用csaw进行peak差异分析，用ChIPseeker进行注释和可视化，用MEME系列进行motif检测和富集，以HMMRATAC进行核小体检测，HINT-ATAC用于footprint分析，用PCEA整合RNA-seq进行调控网络重建。QC：质量检查；TSS：转录起始位点；TF：转录因子；DEG：差异表达基因

预分析：质量控制和比对

ATAC-seq分析的第一步包括比对前QC，read比对到参考基因组，和比对后QC和处理（图 2 A）[ 32 ]。

比对前质量控制

比对前质量控制和read比对步骤是大多数高通量测序技术的标准配置。例如，FastQC [ 39 ]可用于在测序数据中可视化碱基质量得分、GC含量、序列长度分布、序列重复水平、k-mer过高以及引物和衔接子的污染。总体高的碱基质量评分下，read 3'端评分略有下降是可以接受的。与预期的GC含量和read序列长度之间不应该有明显的偏差。此外，在同一实验批次和测序操作的所有样品中，指标应均一。

当前，由于ATAC-seq普遍使用Illumina的Nextera文库，经常会观察到Nextera测序接头比例过高，应将其删除以进行准确的read比对。大多数去除接头的工具采用不同的动态编程，例如 cutadapt [ 40 ]，AdapterRemoval v2 [ 41 ]，Skewer [ 42 ]和trimmomatic [ 43 ] 都需要输入已知的接头序列。例如，对Nextera和Truseq文库使用trimmomatic和内置接头序列是一种直接简单的办法。使用这些工具也可以去除低质量的碱基。根据我们的经验，各种read过滤工具在有效去除低质量和污染接头序列的性能方面通常表现差不多。

比对

过滤后，可以再次执行FastQC，以检查接头和低质量碱基是否已成功移除。然后将过滤的read比对到参考基因组。BWA-MEM [ 44 ]和Bowtie2 [ 45 ] 对于短的双端read存储效率高且快速。两个比对工具的软限位策略允许在read的两端有突出碱基，这可以进一步提高unique mapping rate[ 46 ]。我们建议，unique mapping rate达到80％以上时认为ATAC-seq实验成功。对于哺乳动物物种，基于经验和计算估计，建议染色质开放区域检测和差异分析至少需要5000万mapped read，TF footprinting至少需要2亿[ 10，12，47，48，49 ]。

比对后处理和质量控制

序列比对后，就像大多数DNA测序数据一样，可以使用Picard [ 50 ]和SAMtools [ 51 ] 收集比对BAM文件的基本指标，例如unique mapping reads/rate，duplicated read的百分比以及片段大小分布。此外，如果read比对不正确或mapping质量不佳，则应将其删除。线粒体基因组（由于缺乏染色质包装而更可及 [ 52 ] ）和ENCODE列入黑名单的区域[ 53，54 ]通常具有非常高的read覆盖度，应该去除 [33]。重复的read（很可能已作为PCR产物出现）也应去除，以显着提高生物学的可重复性[ 48 ]。这些步骤将共同提高开放染色质检测的能力，并减少假阳性。

还有其他需要评估的ATAC-seq特定质量指标。通常，成功的ATAC-seq实验应生成片段大小分布图，其具有递减的和周期性的峰，对应于无核小体区域（NFR）（<100 bp）和单核、双核和三核小体（〜200， 400，600碱基对）（图 1 b）[ 9，55 ]。NFR的片段应该在基因的转录起始位点（TSS）周围富集，而核小体结合区域的片段应该在TSS处被形成低谷，TSS周围的侧翼区域会稍微富集（图 1 c）[ 55 ]。可以使用工具ATACseqQC [ 55 ]进行评估。最后，对于正链和负链，read应分别偏移 +4 bp和 -5 bp，以便实现TF足迹和基序的碱基对解析相关分析[ 9，33，56 ]，因为Tn5转座酶对缺口进行DNA修复产生了9 bp重复。大多数上述质量控制和分析报告可以使用MultiQC [ 57 ] 汇总以进行集成的、用户友好的交互式的呈现。

合适工具的选择主要是考虑计算出结果所需的时间。read过滤和比对可能很耗时，并且在速度和准确性之间始终要取舍。根据我们的经验，以下管道的性能相当好：FastQC➔trimmomatic➔BWA-MEM➔ATACseqQC，我们建议这是处理ATAC-seq数据的良好起点。

核心分析：peak calling

ATAC-seq数据分析的第二个主要步骤是识别可及区域（也称为peak），并且是进行高级分析的基础。对于ChIP-seq[ 58，59 ]和DNase-seq [ 60 ]，类似的过程已被全面综述。当前，MACS2是ENCODE ATAC-seq管道的默认call peaks工具。据我们所知，专门针对ATAC-seq开发的call peak工具只有一个[ 61 ]。其他都是从ChIP-seq和DNase-seq中借用过来的，这基于一个假设——ATAC-seq peak模式与ChIP-seq/DNase-seq有相同的性质。因此，我们将集中于当前用于ATAC-seq的工具，并提供潜在的替代品（图 4 a）。

与ChIP-seq不同，ATAC-seq通常没有input control（Tn5转座酶随机切割没有蛋白结合的DNA），因为获得与样本相同覆盖率的input control测序成本较高。因此，需要input control的call peak工具对于ATAC-seq是不切实际的。此外，来自ATAC-seq的双端片段的直接堆积代表了无核小体区（NFR）和核小体结合区（图 3a）。开放染色质区可以通过NFR的短片段堆叠来检测，或使用一个移位延伸的方法——尝试对通过延伸尺寸来平滑化的切割事件进行计数（图 3 B，右框）[ 61，62]。此方法更为通用，因为它可以应用于几乎所有的ChIP-seq call peak工具，并且不受数据片段大小的影响。

图3 核心和高级分析的示意图和实际ATAC-seq数据。a在ATAC-seq实验中，Tn5结合并切割开放的染色质，同时加上接头。对片段进行测序以鉴定开放的染色质区域（黑色）和footprint（蓝色）。NFR片段代表开放的染色质，而核小体结合的片段则反映了核小体的位置（灰色阴影轨迹）。b实际的ATAC序列数据。从BAM文件（原始）生成信号轨，并通过HINT-ATAC校正偏差（校正了偏差）。峰集由三种类型的call peak工具生成：基于计数的（红色），基于形状的（蓝色）和基于HMM（黑色）。对于MACS2，使用两种策略（双端和移位扩展）。对于HMMRATAC，两侧的扩展范围表示核小体。HINT-ATAC轨是由HINT-ATAC检测到的footprint，而RUNX1 motif轨是与JASPAR数据库中的RUNX1 motif匹配的footprint。K562 ChIP-seq轨是ENCODE的RUNX1 ChIP-seq，表明足迹检测可以重现真实的TF结合。右边的方框说明了移位扩展方法。首先，它将两端s-bp移至外部，然后将2s-bp移至内部。C通过ATAC-seq数据进行网络重建的图示。TF的存在可以通过上述方法检测到的motif或footprint来表示。NFR：无核小体区域；TF：转录因子；HMM：隐马尔可夫模型

针对ATAC-seq的热门call peak工具可以分为两大类：基于计数的或基于形状的。基于计数的call peak工具采用不同的统计方法将候选区域中的read分布形状与随机背景进行比较。MACS2 [ 63 ]，HOMER [ 64 ]，和SICER / epic2 [ 65，66，67 ]假定泊松分布，而ZINBA [ 68 ]假设零膨胀的负二项式分布。F-seq [ 69 ]和PeakDEck [ 70 ]使用核密度估计来分析片段分布。SPP [ 71]没有片段分布的假设，但是使用滑窗，根据来自上下游侧翼窗口的片段counts来计算分数。请记住，某些工具（例如F-seq和ZINBA）并未得到积极维护，因此应谨慎使用。当将混合模型聚类应用于生物重复时，JAMM可以更准确地确定峰宽和边界[ 72 ]。通常，基于计数的方法更易于解释和广泛使用。

基于形状的call peak工具当前没有在ATAC-seq中使用，但它们直接或间接利用读取的密度分布信息，并被认为可以改善ChIP-seq中的peak calling [ 73 ]。PICS [ 74 ] 对除counts以外的片段位置进行建模，并计算每个候选区域的富集得分。PolyaPeak [ 75 ] 使用描述峰形的统计数据对峰进行排名。CLC [ 76 ]从正峰和负峰中学习了一种用于峰形的高斯滤波器。

当前，HMMRATAC是ATAC-seq [ 61 ] 独有的call peak工具。它采用三态半监督隐马尔可夫模型（HMM），分别将基因组分割为信号强度高的开放染色质区域，信号强度适中的核小体区域和信号强度低的背景区域。尽管HMMRATAC在计算上更加密集，但其性能优于MACS2和F-seq，并同时提供核小体位置信息。

其他考虑因素应包括call peak工具是否解释了Tn5切割的偏好性以及如何处理生物学重复。类似于DNase-seq，Tn5酶促切割将引入偏误，因为它有结合偏好性[ 30，31，77 ]。这与GC含量相关，并且在call peak时应调整[ 22，56 ]。生物学重复可以提高可重复性并减少假阳性峰。通过合并原始read或合并各个样品的峰集，可以将大多数工具应用于重复。重复也可以使用混合模型来集成[ 72 ]。

这些工具所产生的peak tracks可以在图 3b 中可视化。基于计数的工具表现类似，但是与基于形状的工具相当不同。此外，使用神经网络提取了这些peak的潜在序列特征，并展示出来以概括已知的TF motif。这证实了转录因子通过开放的、可接近的染色质在基因调控中起重要作用 [78，79 ]。参数微调是所有上述工具必不可少的[ 9，33 ]，因为开放染色质的宽度会发生变化[ 32]。将邻近的窄峰接合以产生宽峰的工具（例如MACS2，HOMER和SICER / epic2）也被认为可提供更有意义的结果。但是，迄今为止，尚无关于ATAC-seq call peak工具的全面基准性研究。我们建议使用积极维护的工具（例如MACS2和HOMER）进行call peak。并且如果计算资源足够，HMMRATAC可用于ATAC- seq call peak。

进阶分析

Peak

由于其性质，ATAC-seq揭示了转录调控的多个方面，因此第三个主要步骤涉及四个不同水平的解析：峰，基序，核小体和TF足迹。但是，只有少数工具是专门为ATAC-seq设计的。

Peak差异分析

当前，尚未专门开发用于ATAC-seq数据分析的差分峰分析工具。一种简单的方法是找到候选区域（共有峰或分箱的基因组），进行normalize，对这些区域中的片段进行计数，并与其他条件进行统计学比较[ 80 ]。这可以手动地实现，或使用自动化工具，如基于共有峰或基于滑窗的工具（图 4 b）。

图4 call peak和peak差异分析工具总结。acall peak工具可分为基于计数的，基于形状的，和马尔可夫模型方法。根据使用的统计学方法或模型不同，它们可以进一步划分。b峰差分分析工具可以分为基于峰集的方法和滑动窗口方法。基于峰集的方法根据外部call peak工具和RNA-seq DE包的使用情况进行划分。滑动窗口方法根据使用的统计方法或模型进行划分。ZINB：零膨胀负二项式；HMM：隐马尔可夫模型；DE：差异表达；NB：负二项式

在基于共同峰的工具中，HOMER，DBChIP [ 81 ]和DiffBind [ 82 ]依赖于RNA-seq差异表达（DE）分析包，例如edgeR [ 83 ]，DEseq [ 84 ]或DEseq2 [ 85 ]。 因此，它们都假设负二项式（NB）分布，并且需要生物学复制以估计分散度。有人建议通过合并所有样本来call consensus peak，以减少假阳性差异peak，这是HOMER的默认行为[ 86 ]。但是，DBChIP和DiffBind通过交集或并集操作生成共识峰。但是，交集操作会忽略样本或条件特定性的峰，而并集操作通常导致较低的p值和更多的假阳性[ 86 ]。

滑动窗口方法不需要预先生成的峰集。相反，他们对分箱基因组的所有窗口进行评估，并倾向于产生更多的假阳性，并需要严格的过滤和错误发现率（FDR）控制。PePr [ 87 ]和DiffReps [ 88 ]根据重复的情况决定使用NB测试、G检验或卡方检验。对于更宽的峰，ChIPDiff [ 89，90 ]采用HMM来解释相邻窗口之间的相关性。这三个工具独立于RNA-seq DE分析软件包。相反，csaw是将edgeR框架扩展到分箱基因组[ 91]。滑动窗口方法被认为可以更无偏地估计全基因组的read count，但是需要严格的FDR控制才能正确合并相邻窗口。

当前，大多数研究假设峰区域中的ATAC-seq读数遵循NB分布，RNA-seq数据也是如此。但是，基于形状的ATAC-seq数据差异分析工具并不存在。峰不仅包含read count信息，还包含分布形状轮廓。这对于宽峰尤为重要，因为宽峰可能包含多个局部最大值，并且这些偏移可以指示生物学上相关的扰动，这可以在滑动窗口或基于形状的方法中检测到。尽管尚未进行系统性研究，但我们认为合并形状信息将改善差分峰分析。但是，考虑到对于重复的处理、外部call peak工具依赖性以及后端统计方法，由于csaw的edgeR框架易于解释，因此值得一试。

Peak注释

获得峰集后，峰的注释可将染色质的可及性与基因调控联系起来。通常，峰会被注释到最接近的基因或调控元件。HOMER，ChIPseeker [ 92 ]和ChIPpeakAnno [ 93 ]被广泛用于将peak注释到最接近或重叠的基因、外显子、内含子、启动子、5'非翻译区（UTR），3'UTR和其他基因组特征。ChIPseeker和ChIPpeakAnno还具有丰富的可视化功能来解释注释结果，例如带注释的基因组特征的饼图（图1 d中的示例 ）。通常，来自ATAC-seq的峰将代表包括增强子和启动子在内的不同顺式调控元件的混合[ 12]。获得诸如最接近的基因之类的基因组特征列表后，还可以使用GO[ 94 ]，KEGG [ 95 ]和Reactome [ 96 ] 等数据库进行功能富集分析。通常，峰注释会产生生物学和功能上有意义的结果，以供进一步研究。

Motif

尽管峰注释提供了功能解释，但它不能直接解释潜在的机制。开放染色质可以通过TF的影响转录，即TF通过识别并结合到DNA上的特定序列来促进转录。该序列称为基序（motif），结合位置称为TF结合位点（TFBS）。人类中大约有1600个TF，一半以上具有通过实验或计算得到的基序[ 97 ]。大多数转录因子需要染色质开放才能结合，但某些先驱转录因子可以结合不开放的核小体DNA [ 98，99 ]。转录因子通过与组蛋白或非组蛋白竞争[100，101]以及与辅助因子配合[ 102 ]来调节转录。这些染色质可访问性重塑过程已由Klemm，Shipony和Greenleaf等人在最近的出版物[ 103 ]中详细综述。因此，了解motif使用或活性变化可能有助于破译潜在的调控网络，并确定关键调控元件[ 104 ]。有两种类型的基于motif或基于TF的分析方法：基于序列的motif频率或活性预测，以及针对TF占位的footprint分析（在下一节中讨论）。

Motif数据库和扫描

为了利用基序信息，已经作出很大的努力来从实验方法或计算预测中编译基序序列的数据库。流行的数据库如JASPAR [ 105 ]包含多个种类，并且可以使用应用程序编程接口（API）或Bioconductor的封装[容易地检索到106，107 ]。仅举几个数据库为例，CIS-BP [ 108 ]和TRANSFAC [ 109 ]包含真核TF基序，HOCOMOCO [ 110 ]专注于人和小鼠，RegulonDB [ 111 ]专用于大肠杆菌。。但是，没有中央数据库，该数据库包含全面且一致的基序信息，并且差异可能源自原始ChIP-seq实验和用于进行从头基序发现的软件的差异。

图案信息主要以文本格式存储，例如，作为位置权重矩阵（PWM）。HOMER和Bioconductor软件包TFBSTools [ 112 ]和基序匹配器[ 113 ]能够使用PWM在给定的核苷酸序列中寻找推定的TFBS。PWMScan [ 114 ]提供了一个Web服务器，用于使用Bowtie索引基因组进行快速基序扫描。另一广泛使用的工具是MEME套件[ 115，116 ]，其包括FIMO [ 117 ]来搜索个体基序，MAST [ 118 ]用于从多个基序聚集搜索结果，并且MCAST [ 119]来推断由多个基序形成的调控模块。这些工具基于统计匹配生成推定的TFBS列表。其中，由于MEME套件和PWMScan具有Web应用程序界面，因此更易于访问。

Motif富集和活性分析

基于上述主题搜索工具，可以获得每个峰值区域中主题的位置和频率，并将其与随机背景或其他条件进行比较。HOMER使用超几何检验，而MEME-AME [ 120 ]使用秩和检验来比较峰内的基元频率。MEME-CentriMo [ 121 ]进一步确定了峰中心附近富集的基序。DAStk [ 62 ]生成MD分数（基序位移分数）[ 122 ]。这是通过计算从每个峰中心开始的小窗口（150 bp）与大半径（1500 bp）内基序出现的比率来实现的。MD得分也可以在不同条件下用Z进行比较-测试。这些方法采用不同的统计测试来比较峰和背景区域中的图案频率。

除了超额测试外，每个假定的TFBS的可访问性都被认为与TF活动相关联，并且可以通过片段计数来衡量。ChromVAR [ 56 ]使用Z得分计算每个条件在多个条件下的可及性偏差，并针对已知技术偏差（GC偏差，平均可及性和峰读取分数）进行调整。它是专门针对具有大量可被视为重复单元的scATAC-seq数据而设计的。但是，尚未评估其在批量ATAC-seq中的性能。DiffTF会针对所有TFBS生成可访问性倍数变化的分布，并针对每个基序的GC含量进行调整，然后将其与重排的零背景进行比较以评估显著性[ 123，124 ]。总之，MEME-CentriMo是一种广泛使用的Web应用程序，可以生成可视报告，而在scATAC-seq中chromVAR可作为其替代。

到目前为止提到的所有工具都根据峰区内发现的序列间接预测推定TFBS。此类TFBS可能包含很大比例的假阳性，并且可能不完整且令人困惑。这是因为并非所有TF都具有确定的基序，并且同一家族的TF可以共享非常相似的基序[ 125 ]。而且，预测的富集或活性变化可能没有显著的生物学意义，这妨碍了对基于序列的motif分析结果的解释。

Footprint

解密TF调控的另一种方法是使用footprint。ATAC-seq中的足迹是指活性TF与DNA结合并阻止Tn5在结合位点切割的图式。这使得开放的染色质区域内出现一个信号低谷（图 3 A）[ 47，126，127 ]。因此，活跃结合的TF的足迹可以用于重建专门针对某些样本的调控网络。

但是，ATAC-seq足迹分析存在一些障碍。首先，在预处理步骤中，由于9-bp的重复序列的存在，原始read的shift对于足迹的精准检测很重要9，33 ]。其次，由于Tn5具有结合偏好性[ 32，128 ]，且瞬时TF结合的信号弱[ 129 ]，足迹检测在实验上和计算上都很困难[ 130 ]。在DNase-seq中，人们在footprinting方面已经做出了巨大的努力。除了酶促偏好性不同外，其他方面两者面临着类似的挑战。然而，只有少数足迹工具已经对ATAC-seq进行过测试，尚未有系统性的评估完成[ 48，131，132 ]。

足迹分析工具主要分为两类：从头和以motif为中心的方法。从头计算方法根据典型足迹模式（峰-谷-峰）的特征，预测峰上的所有足迹位置。然后，将这些假定的足迹位置用于匹配已知motif或识别新D的motif。而以motif为中心的方法则需要先验TFBS的输入，并使用监督或无监督的方法将这些位点区分为结合或未结合（表 1）。

表1 footprinting工具的总结，包括软件类别，编程语言，算法或统计方法，DNase-seq或ATAC-seq的偏差校正以及输出的统计量。此外，倒数第二栏举例说明了工具在ATAC-seq数据中的应用

从头(de novo)工具

对于从头方法，在数学上定义什么是足迹并对Tn5的切割偏好性进行去噪非常重要 [128，134 ]。Boyle等人[ 135 ]提出了一种HMM，在每个碱基处使用标准化和平滑化的fragment counts来检测不同的状态，如足迹、侧翼和背景。HINT，HINT-BC（用于DNase-seq偏误校正），和最近HINT-ATAC也采用HMM，但只有HINT-ATAC对链特异性的Tn5切割偏好进行了校正（图 3 b）[ 130，133，134 ]。一个例子如图3b所示：通过HINT-ATAC在白血病样本中检测到的足迹在K562细胞系中通过RUNX1 ChIP-seq得到了验证。由于这些基于HMM的方法需要使用人工注释的基因组区域进行监督式训练，因此需要进一步评估它们在较大数据集中的通用性。Wellington和Wellington-bootstrap[ 136，137 ]比较侧翼和候选足迹区域的Tn5切割数以找到局部极小值。Neph’s method, Boyle’s method, HINT, and Wellington都没有考虑偏好校正，而DNase2TF和HINT-BC按照DNase-seq的模式进行偏好矫正[ 47，129]。参数调整是一个关键的考虑因素，它将影响最终的调用。一种使用HINT和Wellington的优化管线已经被报道，该管线将ChIP-seq结合位点视为真阳性，使用曲线下面积（AUC）分析评估结果[ 48 ]。总之，目前只有HINT-ATAC能够处理ATAC-seq特易的偏好性。

以motif为中心的工具

相比于从头方法，以motif为中心的方法着重于先验TFBS，考虑了TF特异性的足迹轮廓。它面临的挑战是在具有高质量motif的TF丰富时避免确定性偏倚。

无监督的motif中心方法根据基因组区域中提取的特征（例如到TSS距离、PWM匹配分数、序列保守得分[145，146 ]），以及从测序reads中提取的特征（例如putative TFBSs附近的read数量和形状分布[139，140，141，147 ]，将推定的TFBSs分为结合或不结合。其中，CENTIPEDE模型假定read符合多项式分布，其性能对于不同的TF和细胞类型具有特异性的参数敏感性 [133，139，143 ]，而msCentipede和Romulus考虑了这些具有异质性的足迹轮廓[ 140，141 ]。此外，msCentipde可以对Tn5偏好性进行建模，而Romulus可以改善低深度数据和低质量motif的性能。PIQ [ 147 ]使用高斯过程对read分布进行建模，并且在有重复时可以进一步提高鲁棒性。无监督工具的准确性在很大程度上取决于特征选择和构建，因此可以尝试进行特征工程和选择技术，例如单热编码、分箱和聚类，可以进一步提高性能。

相比之下，监督式的以motif为中心的工具需要高质量的ChIP-seq才能注释真实的TFBS，使其成为训练数据。MILLIPEDE和BinDNase都使用logistic回归[ 142，143 ]，而DeFCoM使用支持向量机（SVM），BPAC使用随机森林分类[ 131，144 ]。具体来说，BinDNase为每个TF分别训练一个模型，以解释TF特定的足迹模式。与逻辑回归相比，DeFCoM中使用的SVM方法对异常值的鲁棒性更高[ 131]。此外，在ATAC-seq数据上对DeFCoM进行了测试，与DNase-seq相比，DeFCoM的性能略有下降，读取次数是原来的两倍。对于所有监督式的工具，由于可变的足迹模式，跨TF /细胞类型验证的性能会降低[ 142 ]。这可能会阻碍其在稀有细胞群体或异质癌症样本中的应用。更大量、更多样化的训练数据合在一起可以改善足迹分析的效果[ 144]，并且我们猜测集成学习是有益的，因为培训了多个学习器以进行集体预测。此外，所有这些工具都是使用DNase-seq数据进行训练的，所以应使用ATAC-seq数据对其进行重新训练，以解决不同数据的内在偏差。通常，由于TF和细胞类型特异性的足迹模式具有很大的可变性，所以建模仍然很困难。

如果你对全局TF足迹模式在不同条件之间的变化感兴趣，则可以使用BaGFoot [ 132 ]。在序列深度归一化和偏差校正之后，它会计算所有TF的足迹深度和侧翼可及性。该方法对测序类型（DNase-seq或ATAC-seq）、call peak工具以及偏好校正方法[ 132 ]均很鲁棒。

关于足迹分析的评论

足迹分析有几个注意事项。首先，监督的motif为中心的足迹工具通常优于无监督的对应工具以及从头方法，但是其通用性就稍逊一筹 [130，131 ]。目前已经使用特定细胞类型中特定TF的ChIP-seq和DNase-seq数据对它们进行了训练。因此，它们的context可能无法推广应用于ATAC-seq。此外，感兴趣样品的训练数据不一定有，并且跨TF /细胞类型进行预测应谨慎131，144]。这些工具对ATAC-seq的通用性仍然需要广泛的评估。其次，偏差校正在DNase-seq和ATAC-seq足迹检测中都很重要。最近，Tn5偏好的基序已被鉴定出来，它显示与一些C2H2锌指TF容易混淆 [ 128 ]。第三，能够有效实现足迹分析的ATAC-seq最小测序深度是多少，尚未有通用的指南。虽然建议每个样品read数超过2亿，但有报道称DeFCoM对于更少的测序read数也能有相当的性能 [ 10，48，131 ]。footprinting的性能随着深度的增加而改善，但改善的程度在不同的TF和细胞类型之间不同，因为它们结合的亲和力和脱离力不同 [ 131]]。但是，我们需要进行饱和度分析才能为每种样品的测序深度提供合理的建议。第四，对于低质量motif和novel motif，从头方法仍然具有优势。尽管不同的研究中对足迹方法的评价不一致（因为分析工具的选择、参数设置和评价标准不同），我们认为HINT-ATAC可以是一个不错的选择，因为它具有ATAC-seq特异性的偏好校正[ 130，131]。此外，研究人员可以结合多种工具的结果来获得高度可靠的足迹。尽管如此，ATAC-seq中的足迹分析对于理解TF调控和进一步重建细胞特异性的调控网络很有用，因此需要针对特定情况下的软件比较和开发进行广泛的基准测试。

核小体定位

核小体由组蛋白八聚体和大约147bpDNA组成，通过改变染色质开放性来影响TF结合（图 3 a）2，103，148 ]。在标准ATAC-seq文库中，较长片段对应于核小体相关区域（图 3a）[ 9 ]。已有分析工具可以检测核小体片段富集的区域。但是Schep等人报道，ATAC-seq中的核小体检测要比MNase-seq中更困难，因为染色质开放区以外的read覆盖度更低[ 149 ]。

针对MNase-seq开发的软件如DANPOS2，PuFFIN，INPS，和NucTools，可以在ATAC-seq数据过滤得到核小体相关片段后使用 [149，150，151，152，153 ]，而NucleoATAC和HMMRATAC是专为ATAC-seq开发的。通过将位置信号与V-plots互相关联来计算每个碱基的信号得分，NucleoATAC优于DANPOS2。V-图是对片段大小和中点位置进行可视化的点图，和跨物种[是保守的149，154，155]。信号分数被归一化和平滑化，并且通过对数似然来找到局部最大值。HMMRATAC可以同时检测开放的染色质和核小体相关联的区域，如前文所述（图 3 b）中[ 61 ]。此外，DANPOS2和NucTools可以检测不同条件之间的核小体占位变化和位置偏移[150，151 ]。INPS采用小波去噪方法，而PuFFIN通过片段大小加权累加的核小体片段分布来识别核小体[152，153，156 ]。

但是，所有这些工具都具有典型ATAC-seq实验的相同潜在缺点，即染色质开放区之外的覆盖率较低。将来，将需要新的实验方案以及用于ATAC-seq的生物信息学方法，以更有效和精确地捕获核小体的占位。在这里，我们认为HMMRATAC和NucleoATAC是用于ATAC-seq核小体检测的两个有用且特异性的工具。

与多组学数据集成以重建调控网络

到目前为止，我们已经阐明了ATAC-seq数据分析的特定要求，将ATAC-seq与其他高通量测序技术（如RNA-seq和ChIP-seq）的集成越来越引起人们对基因调控的兴趣。

与ChIP-seq整合

由于开放染色质是大多数TF结合的先决条件，因此ATAC-seq峰通常与TF ChIP-seq峰重叠，但通常更宽。因此，TF ChIP-seq和ATAC-seq可以在同一实验系统中相互验证彼此的质量和可靠性[ 157 ]。TF ChIP-seq的unique peaks可能指示了先驱转录因子，它结合到封闭染色质，然后招募染色质重塑因子或其他转录因子并起始转录[ 98，103 ]。通过结合真实的TF ChIP-seq峰以减少假阳性，可以进一步改善基于推定TFBS的分析，例如motif富集和footprint检测[ 54]。ATAC-seq也可以与组蛋白标记ChIP-seq集成，并发现与活跃染色质标记（H3K4me3的，H3K4me1，H3K27ac等）正相关，与不活跃的染色质标记（的H3K27me3）负相关[9，157，158 ] 。总之，整合ChIP-seq和ATAC-seq有助于理解TF和组蛋白促进了染色质可及性的变化。由于protocol的简便性和样品需求较少，我们预见在特定的TF ChIP-seq实验之前，ATAC-seq可以成为一种预实验方法。

与RNA-seq整合

研究人员还对染色质可及性与RNA-seq基因表达变化的定性或定量关联感兴趣。直观地，研究人员可以发现DE基因在各自的TSS周围是否也具有明显的染色质可及性差异[ 159 ]。此外，可以推定DE基因受到开放染色质中特定基序或足迹相关的TF的调控。在单细胞水平，Litzenburger等试图结合scRNA-seq和scATAC-seq来鉴定TF靶基因（当GATA结合位点可及性改变时，其靶基因的表达量会随之改变）[ 160 ]。Cao等使用LASSO回归模型来识别导致目标基因表达变化的远端峰[ 161]。结合scATAC-seq和scRNA-seq的耦合聚类已被验证可提高亚群检测的准确性[ 162 ]。ATAC-seq与RNA-seq的整合有助于破译基因调控和细胞异质性。

重建调控网络

尽管ATAC-seq可以同时检测数百个TF基序的出现或足迹，但可以通过将足迹/基序与下游基因联系来重建细胞特异性的调控网络。类似的方法已在DNase-seq中被报道（图 3 C）[ 104，163 ]。但是，以前在DNase-seq中的尝试仅限于启动子区域，仅研究TF-TF调控[ 104 ]。启动子内的峰仅占所有ATAC-seq峰的一小部分，而大多数峰位于远端增强子中，从而降低了推断调控网络的能力[ 9]。增强子在线性基因组中可能非常遥远，但在空间上接近其目标基因（3D模式）。这导致增强子的直接靶基因难以预测。许多研究认为远端峰是增强子，并像启动子分析一样将远端峰联系到最近的基因[ 164，165，166 ]。对于scATAC-seq，Pliner等人推出了Cicero，准确地概括了可开放的峰，并将增强子和启动子联系到同一靶基因。该方法已通过正交方法验证[ 167]。尽管已证明Cicero可以用于scATAC-seq，但尚不清楚此方法是否适用于样本量小得多的常规ATAC-seq。尽管如此，Cicero是使用ATAC-seq将远端增强子与基因调控相联系的先驱。

尽管可以单独使用ATAC-seq重建无方向的TF基因调控网络，但当整合RNA-seq时，有向调控可以进一步推断为激活或抑制。Duren等提出了一个基于基因表达和染色质可及性配对（PECA）数据的模型，以预测目标基因的表达与TF表达、染色质重塑因子表达和染色质可及性的关系[ 168 ]。Miraldi等使用ATAC-seq衍生的二元TF-基因相互作用作为先验网络，以进一步完善从RNA-seq数据推断出的调控网络[ 166 ]。Berest等根据全基因组中TFBS的开放性和TF的表达量，将TF分类为激活因子和抑制因子[ 124]]与该可访问性。该分类依赖于一个假设——类似于组蛋白标记，染色质开放性与激活型转录因子的表达量呈正相关，与抑制型TF表达呈负相关[124，169 ]。此方法仅允许以全局方式进行分类。

为了进一步优化网络重建，可以集成可用的ChIP-seq公共数据集以提高footprinting的准确性。从染色质构象数据中整合已知的增强子-启动子相互作用也将有所帮助。随着深度学习的兴起，为了构建有效的算法来预测转录调控网络，需要在特征构建和选择方面进行更多工作。总而言之，将ATAC-seq与多组学数据整合在一起会产生生物学上有意义的结果，从而可以揭示基因调控的潜在机制。

ATAC-seq数据的管道

处理ATAC-seq数据的集成管道日益重要。目前已经开发了几种工具，通过将前面讨论的工具组合在一起，可以将重点放在下游分析的不同方面。

仅举几个例子，esATAC [ 170 ]和CIPHER [ 171 ]专注于peak注释，而图形用户界面（GUI）工具GUAVA [ 172 ]则专注于差异peak检测以及功能注释。ATAC2GRN [ 48 ]是另一种专门为footprinting而优化的管道。

这些管道将为研究人员提供便利而方便的入门，使研究人员可以用最少的编程技能来探索ATAC-seq数据。但是，这些管道的普遍缺乏参数调整的灵活性。大多数参数都是凭经验进行硬编码的，因为它们的组合会随着工具数量的增加而呈指数增长，因而难以在特定的条件下修改管道。总体而言，具有可视化和用户界面的管道将更适合非程序员探索数据。

单细胞ATAC-seq

通过微流体、纳米孔或组合索引技术，scATAC-seq现在能够以低成本和方便的protocol测量数千细胞的染色质可及性[33，34，35 ]。由于每个碱基的染色质可及性将是二元的，scATAC-seq数据将是稀疏的，因为在二倍体生物中只有两个DNA拷贝。这是分析scATAC-seq数据的挑战。尽管针对常规ATAC-seq的分析已列出，单细胞的还有一个重要分析是聚类。Chen等人最近的scATAC-seq聚类方法基准研究表明，SnapATAC，Cusanovich2018和cisTopic优于其他方法[ 23，173，174，175 ]。这三种方法的特点是结合了基于window的基因组分箱、可及性的二值化、覆盖率偏好校正和使用主成分分析降维的流程，专门用于处理稀疏的scATAC-seq数据[ 175 ]。这项研究为将来的scATAC-seq软件开发提供了有用的见解。

最近开发了同时测定同个细胞中的染色质可及性、转录组以及蛋白质组的新技术，如scNMT-seq，sci-CAR，和Pi-ATAC[161，176，177 ]。这些实验的数据可以帮助推断表观基因组、转录组和蛋白质组之间的复杂相互作用，并有助于我们理解为什么不同的细胞表现出不同的行为。尽管单细胞分析的优势很明显，但仍然存在挑战。单细胞技术的成本和时间效率提升以及生物信息学工具开发仍然是活跃的研究和开发领域。

未来展望和总结

ATAC-seq近年来发展迅速，已成为研究染色质可及性的一种可选方法。对单细胞、血液样品和冷冻组织，现在有优化的protocol[26，33，34，35，178 ]。尽管protocol取得了进展，但生物信息学分析工具的进展缓慢，没有定义全面的分析管道。这对今后ATAC-seq结果的解释造成了障碍。

在这篇综述中，我们系统地讨论了ATAC-seq分析流程中的所有主要步骤，从原始read读取开始，到最终生物学意义的解释，供读者参考。在这里，我们提供了可用工具的指南，并提出了可以考虑的分析步骤的建议，以促进对ATAC-seq数据进行合理的生物学解释。比对和QC步骤类似于RNA-seq和ChIP-seq。至于call peak，大多数ChIP-seq衍生工具都与ATAC-seq数据兼容，但是全面的基准测试将有助于选择合适的工具，并指导未来ATAC-seq特异性的call peak工具发展。越来越多的证据表明，当前工具改进或参数化之后可以适用于ATAC-seq数据。

对于下游分析，peak差异分析可以概述染色质可及性的变化。但是，这些变化可能来自read数和峰的形状，并且可以通过基于计数或滑动窗口的方法来检测。这两种方法在ATAC-seq中的性能仍需要进行进一步评估，并且可能对于特定环境具有特异性。为了推断生物学功能和相关的TF，peak注释和motif富集分析是初步了解的首选。

motif和footprint分别是调控事件的直接和间接指标。检测足迹的困难来自酶切偏倚和瞬时TF结合信号微弱。最近的出版物第一次尝试拥抱快速发展的监督式机器学习算法[131，144 ]，以取代用数学公式定义footprint。此外，由于ATAC-seq数据固有的弱点（峰以外的区域read覆盖率很低），核小体检测仍然很困难。NucleoATAC和HMMRATAC曾尝试这样做，但是该领域的检测方法仍存在巨大的空白。

单独用ATAC-seq数据或与多组学数据整合来重建是基因调控网络分析的另一个考虑因素。由于ATAC-seq所需细胞量可以低至500个，尤其是在发育生物学和临床样品中，可以研究明确定义的亚群，这特别诱人。scATAC-seq为研究异质细胞群体中的染色质生物学提供了另一种选择。

综上所述，作为一种信息丰富的研究方法，ATAC-seq十分需要特定的生物信息学分析工具，以便进一步分析染色质状态、TF足迹、核小体位置以及调节网络重建。对于初级分析，我们建议研究者不妨通过组合如下软件来搭建有效的分析流程：用FastQC，trimmomatic和BWA-MEM进行预分析，用MACS2进行peak calling。对于高级分析，我们建议使用csaw进行peak差异分析，使用MEME系列软件进行motif检测和富集，使用ChIPseeker进行注释和可视化，使用HMMRATAC进行核小体检测，使用HINT-ATAC进行足迹分析。如果可获得RNA-seq数据，则可以使用PECA方法重建调控网络。不过，研究者可以根据这篇综述对每个步骤所需的工具进行替换，我们推荐根据实验体系和数据的具体情况选择工具。

我们设想，这篇综述将鼓励研究人员认识到ATAC-seq数据分析的复杂性和当前的主要障碍。新的ATAC-seq特异性的工具和全面的基准研究将使ATAC-seq在不久的将来能够回答更多的生物学问题。

参考文献

见原文，此处略。