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

更好体验点击此：翻译

很长不看系列

作者建议研究人员可以建立一个有效的工作流程，结合 FastQC、trimmomatic 和 BWA-MEM 进行预处理，MACS2 进行 Peak calleing。对于高级分析，作者建议使用 csaw 进行差异 Peak 分析，使用 MEME suite 进行 motif 检测和富集，使用 ChIPseeker 进行注释和可视化，使用 HMMRATAC 进行核小体检测，使用 HINT-ATAC 进行足迹分析。如果 RNA-seq 数据可用，可以使用 PECA 方法重建调控网络。

摘要：

ATAC-seq（Assay of Transposase Accessible Chromatin sequencing）被广泛应用于染色质相关生物学的研究，但对其分析工具的全面综述却很少。本文讨论了关于 ATAC-seq 分析大的主要步骤：预处理（质量控制与比对）、核心分析（Peak calling）、高级分析（Peak 的差异分析以及注释、motif 富集分析、footprnting 足迹分析、核小体位置分析）。本文也回顾了用多组学数据重建转录调控网络的工作，并强调了目前每一步所面临的挑战。最后，本文描述了单细胞 ATAC-seq 的潜力，并强调了开发 ATAC-seq 特异性分析工具以获得生物学上有意义的见解的必要性。

引言：

哺乳动物的 DNA 通过三个主要的等级尺度高度浓缩：第一个是核小体，然后包裹成染色质，形成第三个层次，即染色体[1,2,3,4,5,6]。染色质可以在转录活性的常染色质和非活性的异染色质之间动态切换[7, 8]。三种尺度的 DNA 浓缩及其相互作用都有助于基因调控。

最近的基因调控研究集中在表观遗传学上，并且随着高通量测序技术的进步，已经开发了各种方法来破译表观遗传学图谱。其中包括 ATAC-seq [9, 10]、DNase-seq（DNase I hypersensitive sites sequencing）[11,12,13]、FAIRE-seq（Formaldehyde-Assisted Isolation of Regulatory Elements sequencing）[14] 等这些技术都是用于研究染色质可及性的。ChIP-seq (Chromatin Immuno-Precipitation sequencing)可用于转录因子（TF） [15,16,17] 的结合和组蛋白修饰 [18, 19] 的研究。MNase-seq（Micrococcal Nuclease sequencing）用于检测核小体的定位与占位（nucleosome positioning and occupancy）[20, 21]。这些实验的详细说明见 [22]。

以上这些用于检测染色质可及性的方法，自从 2013 以来，ATAC-seq 成为最普遍的方法。近年来 ATAC-seq 相关的文章与数据呈指数增长，说明了这个技术在研究生物学问题中发挥的重要性 (Fig. 1a)。例如绘制正常哺乳动物组织和细胞类型中的增强子（enhancer）图谱 [23,24,25]、研究正常造血（hematopoiesis）与白血病（leukemia）之间的可及性变化 [26, 27] 以及精神分裂症患者和癌症基因组图谱（TCGA）泛癌队列中的染色质状态 [28, 29]。这一前沿技术在基础和转化研究中的示意图如 Fig. 3a 所示。简而言之，ATAC-seq 整合了基因工程化的高活性 Tn5 转座酶，可同时切割开放的染色质，留下9bp 的交错缺口，并将高通量测序接头连接到这些区域 [ATAC-seq 原理30, 31]。然后进行配对末端测序以促进这些开放区域的更高的唯一比对率（uniq alignment rates） [32]。

Fig. 1 ATAC-seq 数据集增加的概述以及用于预处理和高级分析的样本输出结果。
• a：2013 至 2019 年在 PubMed 上 ATAC-seq 数据集、ATAC-seq 发表文章、DNase-seq 数据集、FAIRE-seq 数据集、MNase-seq 数据集的数目。
• b: 片段大小在 100bp 和 200bp 左右有明显的富集，表示没有核小体结合和单核小体结合的片段。
• c：TSS 富集可视化可以看出，没有核小体结合的片段在 TSS 处富集，而但核小体结合的片段在 TSS 上缺失，在 TSS 两侧富集。
• d：Peak 注释饼图显示，半数以上的 Peaks 位于增强子区域（远端的基因间区和内含子区），只有约 25% 的 Peaks 位于启动子区
• TSS： transcription start site

Tn5 转座酶的高活性使 ATAC-seq 流程成为一种简单、省时的方法，需要 500-50,000 个细胞 [9]。敏感性和特异性与 DNase-seq 相当，但优于 FAIRE-seq，这两种方法都需要数百万个细胞 [9]。由于 ATAC-seq 在制备过程中不涉及严格的大小选择，因此它也可以使用代表核小体单体和多聚体的片段来识别核小体的位置 [9]。近年来，利用FACS（fluorescence-activated cell sorting）、微流体（microfluidic）和 nano-well-based 等方法来开发单细胞 ATAC-seq (scATAC-seq) 技术 [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 callers 和建立调控网络模型 [37, 38]，但迫切需要对 ATAC-seq 数据分析的主要部分进行系统综述。这篇综述将涵盖流程中列出的四个最重要的步骤 (Fig. 2)。这些包括（<u>1）预处理（质量控制 QC 和比对）；（2）核心分息（Peak calling）；（3）对于 Peaks、motifs、nucleosomes、TF footprints 的高级分析；（4）整合多组学数据，重建调控网络。</u> 这些步骤将使研究人员能够对 ATAC-seq 数据进行精确分析，并产生更有生物学意义的结果。最后，我们将介绍 ATAC-seq 分析和scATAC-seq 的挑战和机遇。

Fig. 2 一个典型的ATAC-seq分析的路线图。列出了四个主要步骤，包括预处理、核心分析、高级分析和与多组数据的整合。
• 预处理（Pre-analysis）包括比对前的质量控制 QC（Pre-alignment QC）、比对（Alignment）、比对后处理（Post alignment processing）、QC。
• 核心分析（Core analysis）包括 Peak calling。
• 高级分析（Advance analysis）包括 Peak、motif、footprint、nucleosome 分析。
• 多组学整合包括与 ChIP-seq、RNA-seq 数据的整合以及调控网络的重建。
每个框中的文本强调了每个分析步骤中的重要注意事项。本文建议研究人员使用 FastQC、trimmomatic 以及 BWA-MEM 进行预处理、MACS2 进行 Peak calling、csaw 进行差异 Peak 分析、ChIPseeker 进行注释和可视化、MEME suite 进行 motif 查找和富集分析、HMMRATAC 进行核小体检测分析、HINT-ATAC 进行 footprint 分析、PCEA 进行与 RNA-seq 数据的网络重建分析。
QC: quality check; TSS: transcription start site; TF: transcription factor; DEG: differentially expressed gene

预处理：质量控制和比对

ATAC-seq 分析的第一步包括比对前的 QC、讲 reads 回比到参考基因组、比对后的 QC 和处理 (Fig. 2a) [32]。

比对前的质控控制

比对前 QC 和 reads 比对是大多数高通量测序技术分析的标准步骤。比如，FastQC [39] 可以用来可视化测序数据中碱基质量得分、GC 含量、序列长度分布、序列重复水平、k-mer 过表达以及引物和接头的污染情况。3 '端序列 reads 整体的高碱基质量分数稍微下降是可以接受的。与预期的 GC 含量和序列 reads 长度没有明显的偏差。此外，来自相同的实验批次和测序运行的所有样本之间的 metrics 应该是均匀的。

目前，由于 Illumina 的 Nextera 文库在 ATAC-seq 上的广泛使用，Nextera 测序接头的过表达经常被发现，为了获得准确的读校准，应该将其移除。大多数接头去除工具使用不同的动态编程变体，比如：cutadapt [40], AdapterRemoval v2 [41], Skewer [42] 和 trimmomatic [43] 都需要输入一直的接头序列。：例如，为 Nextera 和 TruSeq 文库使用带有内置接头序列的 trimmomatic 将是一个简单的步骤。使用这些工具也可以消除低质量的碱基。根据作者的经验，reads 修整工具在有效去除低质量和污染的接头序列方面的性能大致相当。

比对

reads 修剪后，可以再次执行 FastQC，检查接头和低质量的碱基是否被成功移除。然后将修剪后的 reads 序列回比到参考基因组。BWA-MEM [44] 和 Bowtie2 [45] 比对软件对于短的双端数据是更节省内存和更快的。来自两个比对软件的 soft-clip 策略允许在 reads 的两端都有碱基突出，这可以进一步提高惟一回比率[46]。作者认为，一个唯一回比率超过 80% 是一个成功的 ATAC-seq 实验的典型。对于哺乳动物物种，根据经验和计算估计，建议的最小回比上的 reads 量为 5,000万 用于开放染色质检测和差异分析，2亿 用于 TF 足迹分析 [10, 12, 47,48,49]。

比对后处理和质量控制

在序列回比之后，与大多数 DNA 测序数据一样，可以使用 Picard [50] 和 SAMtools [51] 收集回比上的 BAM 文件的基本度量，如唯一比对 reads 数目/率、重复 reads 百分比和片段大小分布。另外，如果 reads 不正确或比对质量低，应该删除 read。由于缺乏染色质包装 [52]，线粒体基因组更容易获得，而 ENCODE blacklisted 区域 [53, 54] 通常具有极高的 reads 覆盖，也应该丢弃 [33]。可能作为 PCR 产物出现的重复 reads 也应该被移除，以显著提高生物的可重复性 [48]。这些步骤将共同提高开放染色质检测的能力，并产生更少的误差。

还有其他需要评估的特定于 ATAC-seq 的质量度量。通常，一个成功的 ATAC-seq 实验应该生成一个片段大小分布图，其峰值与无核小体区域 (nucleosome-free regions: NFR) (< 100 bp)和单、二、三核小体(~ 200、400、600 bp) (Fig. 1b) 相对应，呈递减和周期性 [9, 55]。来自 NFR 的片段预计会在基因的转录起始位点 (transcription start site, TSS) 附近富集，而来自核小体结合区域的片段预计会在 TSS 附近被耗尽，在 TSS 附近的侧翼区域会有少量富集 (Fig. 1c) [55]。这些可以通过 ATACs****eqQC [55] 工具进行评估。最后，分别对正链和负链的 reads 进行+4bp和 -5bp的移位，以解释 Tn5 转座酶修复损伤 DNA 所产生的9bp 的重复，并实现TF footprint 和 motif 相关分析的碱基对分辨率 [9, 33, 56]。前面提到的大多数 QC 和分析报告都可以使用 MultiQC [57] 进行集成，以实现聚合的、用户友好的、交互式的表示。

选择合适的工具的主要考虑因素通常是时间。reads 修剪和比对可能会耗费时间，而且速度和准确性之间总是存在取舍。根据作者的经验，以下的性能相当好：FastQC➔ trimmomatic➔BWA-MEM➔ATACseqQC。 作者建议这作为一个处理 ATAC-seq 数据的良好起点。

核心分析：Peak calling

ATAC-seq 数据分析的第二个主要步骤是识别开放区域（也称为 Peak），后续高级分析以此为基础。类似的过程已经在 ChIP-seq [58, 59] 和 DNase-seq [60] 中得到了全面的综述。目前，MACS2 是 ENCODE ATAC-seq 流程的默认 Peak caller 程序。据作者所知，只有一个 Peak caller（HMMRATAC）是专门为 ATAC-seq 开发的 [61]。其他的都是从 ChIP-seq 和 DNase-seq 中采用的，并假设 ATAC-seq Peak 模式具有相同的属性。因此，本文将重点介绍ATAC-seq 中当前使用的工具，并提供潜在替代方案的概述 (Fig. 4a)。

与 ChIP-seq 不同的是，ATAC-seq 的 input controls（其中 Tn5 转座酶随机裂解无蛋白 DNA ）由于测序成本高，通常无法获得可比较的覆盖范围。因此，需要 input controls 的 Peak caller 对于 ATAC-seq 是不切实际的。此外，ATAC-seq 配对末端片段的直接堆积代表无核小体和核小体结合区域 (Fig. 3a)。开放染色质可以通过堆积来自NFRs 的短片段或使用 shift-extend 方法来检测，这种方法试图计算被扩展大小平滑的切割事件 (Fig. 3b, 红色) [61, 62]。这种方法更加通用，因为它可以应用于几乎所有的 ChIP-seq Peak caller 软件，并且不受数据片段大小的影响。

Fig. 3 原理图和来自实际的 ATAC-seq 数据的核心和高级分析。
• a：在 ATAC-seq 实验中，Tn5 结合并切开染色质，同时连上接头。这些片段被排序以识别开放染色质区域（黑色）和足迹 footprint（蓝色）。NFR 片段代表开放染色质，而核小体结合片段反映核小体的位置（灰色阴影轨迹）
• b：真实的 ATAC-seq 数据。信号轨迹由 BAM 文件（ Raw）和 HINT-ATAC（bias corrected）生成。Peak 集是由三种类型的 Peak caller 生成的，基于 count-base（红色）、基因形状（蓝色）以及基于 HMM（黑色）。
• 对于 MACS2，使用了两种策略（paired-end 和 shift-extend）
• 对于 HMMRATAC 两边拓展范围表示核小体
• HINT-ATAC 轨迹是 HINT-ATAC 检测到的足迹
• RUNX1 motif 轨迹是 JASPAR 数据库中与 RUNX1 motif 匹配的足迹
• K562 ChIP-seq 轨迹是来自于 ENCODE 的 RUNX1 ChIP-seq，这表明足迹检测可以再现真实的 TF 结合。
• 右边的框中内容说明了 Shift 和 extend 方法。首先讲两端向外平移 s-bp，然后再向内延伸 2s-bp
• c：阐述了通过 ATAC-seq 的重建网络。TF 的存在可以通过上述方法检测到 motif 或 footprint 来表示。
• NFR: nucleosome-free region; TF: transcription factor; HMM: hidden Markov model

ATAC-seq 的通常的 Peak caller 可分为两大类：count-based 或者 shaped-based****。Count-base 的 Peak caller 使用不同的统计方法来比较候选区域和随机背景下的 reads 分布形状。MACS2 [63]、HOMER [64] 和 SICER/epic2 [65,66,67] 假设泊松分布，而 ZINBA[68] 假设负二项分布****零膨胀模型。F-seq [69] 和PeakDEck [70] 使用核密度估计来分析片段分布。SPP [71] 没有对片段分布的假设，而是使用一个滑动窗口来计算来自上下侧翼窗口片段计数的得分。有些工具，如 F-seq 和 ZINBA，并没有得到积极的维护，因此应该谨慎使用。将混合模型聚类应用于生物学重复时，JAMM 可以更准确地确定峰宽和边界 [72]。一般来说，基于计数的方法更容易解释和广泛使用。

shaped-based 的 Peak caller 目前不在 ATAC-seq 中使用，但是它们直接或间接地利用 reads 密度轮廓信息，并且被认为可以改进 ChIP-seq 中的 Peak caller [73]。PICS [74] 对片段位置进行建模，而不是计数并计算每个候选区域的富集分数。Polya Peak [75] 使用描述峰形的统计数据对峰进行排序。CLC [76] 从正峰和负峰中学习高斯滤波峰形。

目前，HMMRATAC 是唯一的专门为 ATAC-seq 的 Peak caller [61]。该方法采用 three-state 半监督隐马尔可夫模型 (HMM)，将基因组同时分割为高信号的开放染色质区、中等信号的核小体区和低信号的背景区。虽然HMMRATAC 的计算量更大，但它的性能优于MACS2 和 F-seq，并同时提供更多的核小体位置信息。

其他需要考虑的因素包括 Peak caller 是否导致了 Tn5 裂解偏差，以及它如何处理生物学重复。与 DNase-seq 类似，Tn5 的酶切会由于结合偏好而引入偏倚 [30, 31, 77]，这与 GC 含量有关，在 calling peak 时应进行调整 [22, 56]。生物学重复可以提高可重复性和减少假阳性 Peak。大多数工具都可以扩展到重复，方法是汇集原始 reads 数据，或者从单个样本中合并 Peak。重复也可以使用混合模型进行集成 [72]。

由这些工具生成的 Peak tracks 可以在 Fig. 3b 中显示出来。Count-base 的工具找到的 Peak 类似，但与 Shaped-based 的工具有很大的不同。此外，利用神经网络提取这些 Peak 的潜在序列特征，并显示出对已知 TF motif 的概括。这证实了 TFs 通过开放染色质在基因调控中发挥重要作用 [78, 79]。由于随着开放染色质的宽度 [32]的变化，对上述所有工具来说，参数的微调都是必不可少的 [9, 33]。将附近的窄峰缝合成宽峰的工具，如 MACS2、HOMER 和 SICER/epic2，也被认为可以提供更有意义的结果。但是，到目前为止，还没有针对 ATAC-seq 的 Peak caller 的全面的基准研究，作者建议使用积极支持的工具，如 MACS2 和 HOMER 来进行 Peak calling，如果计算资源足够的话，HMMRATAC 可以用于 ATAC-seq 的 Peak calling。

高级分析

Peaks

由于 ATAC-seq 本质上揭示了转录调控的多个方面，第三个主要步骤涉及四个不同层次的解释：Peak、motif、nucleosome 和 TF footprint。然而，只有少数工具是专门为 ATAC-seq 设计的。

差异 Peak 分析

目前，还没有专门针对 ATAC-seq 数据分析开发的差异 Peak 分析工具。一种直接的方法是找到候选区域(共同 Peak 或将基因组分 bin： consensus peaks or binned genome )，将这些区域的片段标准化并计数，并与其他条件进行统计比较 [80]。这可以通过手动或使用自动化工具来实现，比如 consensus peak 或基于滑动窗口的工具 (Fig. 4b)。

Fig. 4 Peak calling 和差异 Peak 分析工具的总结。
a：Peak caller 可以分为 count-based、shape-based、Markov model 三类。他们可以根据使用的统计方法或者模型进一步划分。
b：差异 Peak 分析工具可以被划分为基于 Peak 的差异分析和对全基因组划分滑动窗口的差异分析。
• 基于 Peak 的差异分析可以分为外部的 Peak caller 和 基于 RNA-seq DE 差异分析
• 基因滑动窗口的方法根据所使用的统计方法和模型进行划分
ZINB: zero-inflated negative binomial; HMM: hidden Markov model; DE: differential expression; NB: negative binomial

在基于共同 Peak 的工具中，HOMER、DBChIP [81] 和DiffBind [82] 依赖于 RNA-seq 差异 (DE) 分析包，如edgeR [83]、DESeq [84] 或 DESeq2 [85]。因此，它们都假设负二项（NB）分布，并且需要生物学重复以估计离散度。建议通过合并所有样本来 call 共同 Peak 以减少假阳性差异 Peak，这是 HOMER 的默认参数 [86]。但是，DBChIP 和 DiffBind 通过交集或并集操作生成共同峰 Peak。但是，相交操作会忽略样品或特定条件的 Peak，而并集操作通常会显示出较低的 P 值和更多的假阳性 [86]。

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

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

Peak 注释

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

Motifs

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

Motif 数据库和扫描

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

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

Motif 富集和活性分析

基于前面提到的 motif 搜索工具，可以得到 motif 在每个 Peak 区域的位置和频率，并与随机背景或其他条件进行比较。HOMER 使用超几何检验，而 MEME-AME [120] 使用秩和检验来比较 Peak 内的 motif 频率。MEME-CentriMo [121] 进一步确定了在 Peak 中心附近富集的 motif 。DAStk[62] 生成一个 MD 值 (motif displacement score) [122]。这是通过计算小窗口 (150 bp) 到大半径 (1500 bp) 的 motif 出现率来实现的。MD 值也可以通过 Z-检验 在不同情况下进行比较。这些方法采用不同的统计检验来比较在 Peak 区域和背景区域 motif 频率。

除了 overrepresentation test 外，每个假定的 TFBS 的可及性都被认为与 TF 活性相关联，并且可以通过片段计数来衡量。ChromVAR [56] 使用 Z-score 计算每个 motif 在多个条件下的可及性偏差，并根据已知的技术偏差进行调整（ GC bias，Peaks 中的平均可及性和 reads 比例 ）。 它是专门为scATAC-seq 数据设计的，具有大量可以被视为重复的细胞。然而，它在 bulk ATAC-seq 中的性能还没有得到评估。DiffTF 为所有 TFBSs 生成可及性倍数变化的分布，并根据每个 motif 的 GC 含量进行调整，然后与置换后的空背景进行比较，以评估其重要性 [123, 124]。总之，MEME-CentriMo 是一个广泛使用的 web 应用程序，它可以生成可视化报告，而 chromVAR 可以作为 scATAC-seq 的替代方案。

到目前为止所提到的所有工具都间接地从 Peak 区域内发现的 motif 来预测假定的 TFBSs。这种 TFBSs 可能包含大量的误报，并且可能是不完整的和混淆的。这是因为并不是所有的 TFs 都有相同的 motif，来自同一家族的 TFs 可以共享非常相似的 motif [125]。此外，预测的富集或活性变化可能具有微不足道的生物学意义，这妨碍了基于序列的 motif 分析结果的解释。

Footprints

解释 TF 调控的另一种方法是使用 footprint 。ATAC-seq 中的 footprint 指的是一个活跃的 TF 与 DNA 结合并阻止 Tn5 在结合位点内裂解的模式。这在开放的染色质区域留下了一个相对的消耗 (Fig. 3a) [47, 126, 127]。因此，活性结合 TFs 的足迹可以用来重建特定样本的调控网络。

然而，ATAC-seq 足迹分析存在障碍。首先，重要的是在预处理步骤中改变原始 reads，以解释 9-bp 的重复，以便准确地检测足迹 [9, 33]。其次，由于 Tn5 的结合偏好 [32, 128] 和瞬态 TF 结合的弱信号 [129],，足迹检测在实验和计算上都存在困难 [130]。DNase-seq 足迹的研究已经取得了很大的进展，除了酶的差异外，它也面临着相似的挑战。然而，在 ATAC-seq 上只测试了少量的足迹工具，并没有进行系统的评价 [48, 131, 132]。

足迹分析工具主要分为两类: De novo 和motif-centric。De novo方法根据典型足迹模式 (peak-dip-peak) 的特征，预测所有跨越 Peak 的足迹位置。然后这些假定的足迹位点被用来匹配已知的 motifs 或识别新的 motifs。相反，以 motif 为中心的方法需要先验 TFBSs 的输入，并使用监督或非监督方法将这些站点区分为有界或无界（Tabel 1）。

Tabel 1 足迹工具概述，包括软件类别、编程语言、算法或统计方法、DNase-seq 或 ATAC-seq 的偏差校正和输出统计。另外，最后一篇专栏文章举例说明了工具在 ATAC-seq 数据中的应用
FDR false discovery rate, HMM hidden Markov model, SVM support vector machine

De novo tools

对于 de novo 方法，重要的是数学上定义什么是 footprint 并从 Tn5 裂解偏差中去除 footprint 模式 [128, 134]。Boyle 等人 [135] 提出了一种隐马尔可夫模型，该模型在每个基上使用规范化和平滑的片段计数来检测不同的状态，如足迹、侧翼和背景。提示，HINT-BC ( DNase-seq 的偏差校正) 和最近的 HINT-ATAC 也使用 HMM，但只有 HINT-ATAC 对特定于 strand 的 Tn5 裂解偏差进行了校正 (Fig. 3b) [130, 133, 134]。Fig. 3b 显示了一个例子，其中 HINT-ATAC 在白血病样本中检测到的足迹也在带有 RUNX1 ChIP-seq 的 K562 细胞系中得到验证。由于这些基于 HMM 的方法需要使用人工标注的基因组区域进行监督训练，因此需要进一步评估它们在更大数据集中的通用性。Wellington 和 Wellington-bootstrap [136, 137] 比较了侧翼和候选足迹区域的 Tn5 切割数量，以找到局部最小值。Neph 的方法、Boyle 的方法、HINT 和 Wellington 没有考虑偏差校正，而 DNase2TF 和 HIN-BC 确实考虑了 DNase-seq 的偏差校正[47, 129]。参数调优是一个重要的考虑因素，它将影响最终的 calling。提出了一种基于 HINT 和 Wellington 的管道优化方案，该方案将 ChIP-seq 结合位点视为真阳性，并利用曲线下面积 (AUC) 分析对结果进行了评价 [48]。总之，目前只有 HINT-ATAC 处理 ATAC-seq 特定的偏差。

Motif-centric tools

Motif-centric 以 motif 为中心的方法侧重于先验的 TFBSs，并考虑与从头开始的方法相比，特定于 TF 的内存占用情况。面临的挑战是如何避免带有高质量 motif 的 TFs 的确定偏差。

基于从基因组区域提取的特征，非监督 motif-centric 方法将假定的 TFBSs 分为结合型和非结合型，例如：与 TSS 的距离、PWM 匹配得分、序列保守得分 [145, 146]，以及从序列 reads，例如，reads 数量和形状分布周围的假定 TFBSs [139,140,141, 147]。其中，CENTIPEDE 模型采用 multinomial model reads 分布，其性能对 TF 和 cell-type-specific 方式的参数敏感 [139,140,141, 147]，而 msCentipede 和 Romulus 解释了这些异构的足迹轮廓 [140, 141]。此外，msCentipde 可以对 Tn5 偏差进行建模，而 Romulus 可以改善低深度数据和低质量 motif 的性能。**PIQ **[147] 使用高斯对 reads 分布进行建模，并在提供重复时进一步增强了鲁棒性。非监督工具的准确性在很大程度上依赖于特征选择和构造，因此可以尝试使用特征工程和选择技术，例如 one-hot encoding、binning 和 clustering，来进一步提高性能。

相比之下，有监督的 motif-centric 工具需要高质量的 ChIP-seq 来将真正的 TFBSs 注释为训练数据。MILLIPEDE 和 BinDNase 都使用 logistic 回归 [142, 143]，而 DeFCoM 使用支持向量机 (SVM)， BPAC 使用随机森林分类器 [131, 144]。具体来说，BinDNase 分别为每个 TF 训练一个模型，以说明特定于 TF 的足迹模式。在 DeFCoM 中使用的支持向量机方法与 logistic 回归相比，对异常值具有更强的鲁棒性 [131]。此外，DeFCoM 在 ATAC-seq 数据上进行了测试，结果显示，与 DNase-seq 相比，DeFCoM 的性能略有下降，reads 增加了一倍。对于所有受监督的工具，由于不同的足迹模式，cross-TF/cell-type 验证的性能都会下降 [142]。这可能会妨碍它们在罕见细胞群或异种癌症样本中的应用。一组更大、更多样化的训练数据被证明可以改善足迹的表现 [144]，我们也希望整体学习是有益的，即训练多个学习者共同预测。此外，所有这些工具都是使用 DNase-seq 数据进行训练的，因此应该使用 ATAC-seq 数据进行再训练，以考虑不同数据的固有偏差。一般来说，由于 TF 和 cell 类型特定的足迹模式具有很大的可变性，因此对它们进行建模仍然很困难。

如果对整体 TF 足迹模式在不同条件之间的变化感兴趣，可以使用 BaGFoot [132]。在序列深度归一化和偏差校正后，计算所有 TF 的足迹深度和侧翼可及性。该方法对分析类型 ( DNase-seq 或 ATAC-seq )、Peak caller 和偏差校正方法都具有鲁棒性 [132]。

对 footprint 分析的评论

对于足迹分析有几点需要注意的地方。

• 首先，有监督的以 motif 为中心的足迹工具通常比无监督的同类工具和De novo 的方法表现得更好，但其通用性较差 [130, 131]。他们使用来自特定细胞类型的特定 TFs 的 ChIP-seq 和 DNase-seq 数据进行训练。因此，它们的内容上可能不是通用的，也不适用于 ATAC-seq。此外，训练数据并不总是可以从感兴趣的样本中获得，cross-TF/cell-type 的预测应该谨慎进行 [131, 144]。这些工具对 ATAC-seq 的通用性仍然需要广泛的评估。

• 其次，偏差校正在 DNase-seq 和 ATAC-seq 足迹检测中都很重要。最近，Tn5 偏好的 motif 已被确定，并显示混淆了一些 C2H2 锌指 TFs [128]。

• 第三，目前还没有一个通用的指南，以实现有效的足迹最小的 ATAC-seq 测序深度。尽管建议每个样本的 reads 量超过 2亿 reads，但 DeFCoM 被描述为在较少的测序 reads 情况下工作得相当好 [10, 48, 131]。随着深度的增加，由于结合亲和力和转化率的不同，TFs 和 细胞类型之间的足迹改善也不同 [131]。然而，饱和度分析需要为每个样本的测序深度的实验设计提供合理的建议。

• 第四，de novo 方法对于低质量和 novel motifs 仍然具有优势。尽管由于所选择的分析工具、参数设置和评价指标，不同研究对足迹方法的评价并不一致，但作者认为，由于 HINT-ATAC 具有特定于 ATAC-seq 的偏差校正，因此它可能是一个不错的选择[130, 131]。

• 此外，研究人员可以结合多种工具的结果来获得高度可靠的足迹。

  尽管如此，ATAC-seq 中的足迹分析对于理解 TF 调控和进一步重建细胞特异性调控网络是有用的，因此需要在特定的内容中进行广泛的软件比较和开发的基准测试。

核小体定位

核小体由一个组蛋白八聚体和大约 147 bp 的DNA复合物组成 (Fig. 3a)，通过改变染色质的可及性来影响 TF 结合 [2, 103, 148]。在标准的 ATAC-seq 库中，较长的片段对应于核小体相关区域 (Fig. 3a) [9]。已经开发了一些分析工具来检测这些核小体片段富集的区域。然而，Schep 等人指出，ATAC-seq 中的核小体检测比 MNase-seq 中的数据更困难，这是由于 ATAC-seq 中的 reads 信号覆盖超出了开放的染色质区域 [149]。

针对 MNase-seq 开发的软件，如 DANPOS2、PuFFIN、iNPS、NucTools 等，可以对核小体相关片段进行过滤后应用于 ATAC-seq 数据 [149,150,151,152,153]，而 NucleoATAC 和 HMMRATAC 则是专门针对 ATAC-seq 开发的。NucleoATAC 通过将位置信号与 V 形图交叉关联，为每个碱基设计一个信号评分，其表现优于 DANPOS2。V 型图是一个点图，用于显示片段大小和中点位置，并在物种间保守[149, 154, 155]。对信号分值进行归一化和平滑处理，用对数似然法求得局部极大值。如前所述，HMMRATAC 可以同时检测开放染色质和核小体相关区域(Fig. 3b) [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 进行整合

因为开放的染色质是大多数 TFs 结合的前提条件，所以 ATAC-seq Peak 通常与 TF ChIP-seq Peak 重叠，但通常更宽。因此，TF ChIP-seq 和 ATAC-seq 可以在同一个实验系统中相互验证彼此的质量和可靠性 [157]。TF ChIP-seq 中独特的 Peak 可以指示先锋 TFs 与封闭的染色质结合，然后招募染色质重组或其他 TFs 并启动转录 [98, 103]。基于假设的 TFBS 的分析，如 motif 富集和足迹检测，可以通过合并真实的 TF ChIP-seq Peak 来减少假阳性，从而进一步改进 [54]。ATAC-seq 也可与组蛋白修饰 ChIP-seq 整合，发现其与活性染色质标记(H3K4me3、H3K4me1、H3K27ac 等)呈正相关，与非活性染色质标记 (H3K27me3) 负相关[9, 157, 158]。综上所述，ChIP-seq 和 ATAC-seq 的整合有助于理解 TF 和组蛋白促进染色质可及性的改变。作者认为 ATAC-seq 将成为特异 TF ChIP-seq 之前的先锋检测方法，因为该方法简单，而且样本要求更低。

与 RNA-seq 进行整合

研究人员还对通过 RNA-seq 定性或定量地将染色质可及性变化与基因表达变化联系起来。直观地，研究人员可以发现 DE 基因是否在相应的 TSS 周围也有显着差异的染色质可及性 [159]。此外，可以推断 DE 基因受开放染色质中与特定的 motif 或 footprint 相关的 TF 调控。在单细胞水平上，Litzenburger 等人试图结合 scRNA-seq 和 scATAC-seq 来识别当 GATA 结合位点可及性改变时其表达变化的靶基因 [160]。Cao 等人使用 LASSO 回归模型识别导致靶基因表达改变的远端 Peak [161]。结合 scATAC-seq 和 scRNA-seq 的耦合聚类方法可以提高亚群检测的准确性 [162]。整合 ATAC-seq 和 RNA-seq 有助于破译基因调控和细胞异质性。

重建调控网络

虽然 ATAC-seq 可以同时检测到数百个 TF motif 或位点，但通过将 footprint/motif 与下游基因连接，可以重建细胞特异性调控网络 (Fig. 3c) [104, 163]。然而，以往在 DNase-seq 方面的研究仅限于启动子区域，仅研究 TF-TF 的调节 [104]。启动子内的 Peak 只占所有 ATAC-seq Peak 的一小部分，而大多数位于远端增强子内，降低了推断调控网络的能力 [9]。增强子在线性基因组中可能非常遥远，但在空间上接近(3D)它们的目标基因。这就导致了直接预测增强子目标基因的困难。许多研究认为远端 Peak 是增强子，并将其与最接近的基因联系起来，类似于启动子分析 [164,165,166]。Pliner 等人利用 scATAC-seq 提出了Cicero，它可以准确地再现co-accessible Peak，并将增强子和启动子链接到相同的靶基因。通过正交试验验证了该方法的有效性 [167]。虽然它已被证明在 scATAC-seq中有效，但尚不清楚该方法是否适用于样本量小得多的 bulk ATAC-seq。然而，Cicero 是使用 ATAC-seq 将远端增强子与基因调控联系起来的先驱。

虽然单独使用 ATAC-seq 可以重建无定向的 TF-基因调控网络，但当整合 RNA-seq时，定向调控可以进一步推断为激活或抑制。Duren 等人提出了一种利用配对基因表达和染色质可及性 (PECA) 数据预测靶基因表达与转录因子表达、染色质重构因子表达和染色质可及性的关系模型 [168]。Miraldi 等人使用ATAC-seq 衍生的二元 TF-基因相互作用作为先前的网络，进一步完善从 RNA-seq 数据推断出的调控网络 [166]。Berest 等人根据整个基因组 TFBSs 上 TF 表达与可及性的相关性，将 TFs 分为激活因子或抑制因子 [124]，并假设可及性与组蛋白修饰相似，与激活因子的 TF 表达呈正相关，与抑制因子呈负相关 [124, 169]。这种方法只允许以全局方式进行分类。

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

ATAC-seq 数据分析流程

越来越需要集成流程来处理 ATAC-seq 数据。已经开发了一些工具，但是通过将前面讨论的工具拼接在一起，对下游分析有不同的重点。

例如，esATAC [170] 和 CIPHER [171] 侧重于 Peak 注释，而图形用户界面(GUI)工具 GUAVA [172] 则侧重于差异 Peak 检测和功能注释。ATAC2GRN [48] 是另一个专门针对足迹优化的流程。

这些流程将为具有会点点****编程技能的研究人员探索 ATAC-seq 数据提供一个有用和方便的入口。然而，这些流程的一个普遍问题是缺少参数调优的灵活性。大多数参数都是根据经验硬编码的，因为它们的组合随着工具的数量呈指数级增长，这使得在任何给定的上下文中都很难修改流程。总的来说，具有可视化和用户界面的流程将更适合于非程序员来研究数据。

单细胞 ATAC-seq

利用 microfluidic、nano-well 和 combinatorial indexing 技术，scATAC-seq 现在能够以较低的成本测量数千个细胞的染色质可及性 [33,34,35]。每个碱基上的染色质可及性是二元的（binary ），scATAC-seq 数据将是稀疏的，因为在二倍体生物中，只有两个 DNA 拷贝。这是分析 scATAC-seq 数据的一个挑战。尽管列出了针对 bulk ATAC-seq 的分析，但针对单细胞的另一个重要分析是聚类。Chen 等人最近对 scATAC-seq 中聚类方法的基准研究表明，SnapATAC、Cusanovich2018 和 cisTopic 优于其他方法 [23, 173,174,175]。具体处理稀疏 scATAC-seq 数据的这三种方法的特点是：基于窗口的基因组结合、可及性的特征二值化（binarization）、覆盖偏倚校正、主成分分析降维 [175]。本研究为未来 scATAC-seq 软件开发提供了有益的启示。

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

未来展望和结束语

近年来 ATAC-seq 得到了迅速的发展，已成为研究染色质可及性的一种选择方法。现在已经有了优化的方案，可以处理单个细胞、血液样本和冷冻组织，提高了信噪比 [26, 33,34,35, 178]。尽管在流程方面取得了进展，但生物信息学分析工具的进展缓慢，没有确定全面的分析流程。这给 ATAC-seq 结果的解释带来了当前和正在进行的障碍。

在这篇综述中，作者系统地讨论了 ATAC-seq 分析流程中的所有主要步骤，以供读者考虑，从原始测序解读开始，一直到生物学意义解释的终点。在这里，作者提供了一个可用工具的指南和建议的分析步骤，以考虑促进 ATAC-seq 数据的适当生物学解释。校准和 QC 步骤类似于 RNA-seq 和 ChIP-seq。至于 Peak calling，大多数 ChIP-seq y衍生工具都与 ATAC-seq 数据兼容。然而，全面的基准测试将有助于选择合适的工具，并指导特定于 ATAC-seq 的 Peak callers 的未来发展。越来越多的证据表明，现有工具的改进或参数化可以适用于 ATAC-seq 数据。

对于下游分析，差异 Peak 分析可以提供染色质可及性变化的概况。然而，这些变化可以从 reads 数目和 Peak 的形状中产生，并且可以通过基于计数或滑动窗口的方法检测到。这两种方法的性能仍然需要在 ATAC-seq 中进行进一步的评估，并且可能针对特定的环境。为了推断生物学功能和相关的 TF，Peak 注释和 motif 富集分析是初步了解的很好的初次洞察力分析。

Motifs 和 footprint 分别是调控事件的直接和间接指标。检测足迹的困难来自酶切偏差和瞬态 TFs 的瞬时信号。最近的文章没有使用数学公式来定义足迹，而是首次尝试采用监督学习的机器学习算法的快速发展 [131, 144]。此外，由于 ATAC-seq 数据固有的弱点，检测核小体仍然很困难，在这种情况下，通常 reads 的覆盖率低于 Peak。NucleoATAC 和 HMMRATAC 已经尝试过了; 但是，在这方面的检测方法仍然有很大的差距。

分析的另一个考虑是利用 ATAC-seq 数据或与多组学数据整合重建基因调控网络。这是特别吸引人的，因为ATAC-seq 可以处理低至 500 个细胞，并允许研究明确的亚群，特别是在发育生物学和临床样本。scATAC-seq 为研究异质细胞群体中的染色质生物学提供了另一种选择。

综上所述，ATAC-seq 是一种信息丰富的分析方法，在分析染色质状态、TF 足迹、核小体位置和调控网络重建等方面，对特定的生物信息学分析工具有着巨大的需求。作为一个起点，作者建议研究人员可以建立一个有效的工作流程，结合 FastQC、trimmomatic 和 BWA-MEM 进行预处理，MACS2 进行 Peak calleing。对于高级分析，作者建议使用 csaw 进行差异 Peak 分析，使用 MEME suite 进行 motif 检测和富集，使用 ChIPseeker 进行注释和可视化，使用 HMMRATAC 进行核小体检测，使用 HINT-ATAC 进行足迹分析。如果 RNA-seq 数据可用，可以使用 PECA 方法重建调控网络。但是，研究人员可以随时参考这篇综述以找到每一步的替代工具，我们建议根据实验系统的内容和收集到的数据来选择工具。

作者预计这篇综述将鼓励研究人员认识到 ATAC-seq 数据分析的复杂性和当前的主要障碍。在不久的将来，新的 ATAC-seq 专用工具和全面的基准研究将使更多的生物学问题得到回答。

从早上到现在翻译（谷歌）完这一篇文献收获很大，ATAC-seq 数据的分析并没有我们考虑的那么简单，如果要做的深的远不止我们认为的那点东西。应该是目前最全的 ATAC-seq 数据指南综述吧。

​