入门ATAC-seq，你需要知道什么？
我们有一个《明码标价》专栏，汇聚了几百个生信工程师提供各种各样的ngs组学数据分析服务：
今天，我们对基因组的了解依然知之甚少。然而正因如此，我们的工作每天都充满了乐趣，因为这是一条到处都充满了惊喜的发现之旅。这是一张非常奇妙的卫星地图，里面孕育着关于人类的几乎所有的秘密。我希望能够推动中国在这个领域的发展，在这张地图上做出新的发现并标上我们自己的图标。
——颉伟
为什么要研究基因组的开放区域？清华大学生命学院研究的表观遗传的颉伟研究员的一段话可能对你有所启发。克里克和沃森早在1953年就发现了DNA双螺旋结构，人类早在2001就将基因组的所有碱基都测完了，但科学家逐渐发现这张地图还不够完整，我们需要新的工具来寻找这张地图上的宝藏。
为什么需要ATAC-seq？
DNA是中心法则中的首要环节，其ATCG四个碱基的排列组合决定了你是猩猩还是人类，而由受精卵发育而来的个体的每个细胞都拥有同样的DNA序列，但是为何同样DNA序列的细胞的表型会如何不同，为何肝细胞是肝细胞，神经细胞是神经细胞？是什么造成了他们生产蛋白不同，决定蛋白生成的RNA不同呢？原因可以用表观遗传来解释。
在DNA到转录成RNA的过程不单只是转录这么简单，它的得到转录许可的路上有许多的关卡需要突破，这包括：染色质可及性，DNA修饰，组蛋白修饰等等。即是说在细胞里的DNA并非所有都可以被表达出来，不同细胞的关卡不一样，造成了不同的表型。
而这其中，染色质可及性即DNA上开放的区域，与细胞属性密切相关，因为其对基因的激活和抑制尤为重要。在许多真核生物（例如人类）中，基因组在核小体（染色质)的帮助下紧密排列和组织。核小体是由8个组蛋白组成的复合物，每个核小体占了约147bp的DNA。当DNA正在被转录成RNA的时候，DNA将从核小体复合物松开。许多因素，如染色质结构、核小体的位置和组蛋白修饰，在染色质的组织和可及性中起着重要作用。
而ATAC-seq就是用来检测染色质可及性的的测序方法，也是一种确定基因表达调控机制的方法。该方法可以帮助识别启动子区域、潜在的增强子或抑制子。启动子是靠近转录起始点(TSS)的DNA区域。它包含转录因子的结合位点，这些转录因子将招募RNA聚合酶。增强子是可位于启动子下游或上游1Mb的DNA区域。当转录因子与增强子结合，并与启动子区域接触时，该基因的转录增加。相反，抑制子会减少或抑制基因的表达。
对于染色质可及性的检测方法一开始并不是ATAC-seq，还有Faire-Seq和DNase-Seq等，但由于ATAC-seq比等相似的技术更容易、更快和更少地需要细胞。
为了找到开放染色质区，基因组被TN5转座酶处理。在ATAC-Seq中，修饰后的TN5将与NextEra接头相对应的DNA序列插入到基因组的开放区域，同时，DNA被转座酶的活性剪切。
这也就决定了ATAC-seq与Chip-seq call出来的peak代表的意义不同。Chip-seq peak是被目的蛋白结合拉下来的DNA，一般只有一个峰，而ATAC-seq是被Tn5转座酶切开、没有被组蛋白结合、染色质开放的DNA位点，如果是TF结合的区域，一般会有一个山谷般的存在。ChIP-seq和ATAC-seq在TF或者Tn5结合区域都会形成一个双峰的reads结合模式，但在判断peak的时候，会有不同的标准。chip-seq是由于TF一起沉淀下来的DNA fragment一般会大于TF的结合区域，read的位置并不是真实的TF结合的位置，需要向内shift；而ATAC-seq也需要shift，但一般是往两边shift。
得到DNA片段后，可以为测序准备建库了，这包括用完整的NextEra接头和纯化、PCR扩增等步骤。基于上述原因，ATAC-Seq推荐使用双端配对的方法。
了解完了什么是ATAC-seq过后，我们来看下ATAC-seq分析的具体用途。
ATAC-seq能分析得到什么？
ATAC-seq可用于：
生成表观基因组图谱 得到在不同组织或不同条件下对应可及性区域。 得到核小体位置 鉴定重要转录因子 生成转录因子结合区域的特征(footprinting)
组合分析：
ATAC-seq+RNA-seq： 一般来说，RNA-seq会优先于ATAC-seq先测，但是差异基因富集出来的基因通路只是一种相关性。想要分析出其中是谁调控了目的基因，就可以通过ATAC-seq来做motif分析，寻找潜在的调控因子，然后再进行后续的实验验证或者chip-seq的验证。
+HiC： 对于一些想了解染色质高级结构对生命行为的作用的时候，通常会需要用到ATAC-seq等技术，因为Hi-C分析得到高级结构compartmentA/B、TADs、Loops等信息，通常只是相关性，但通过ATAC-seq，你可以获得promoter、enhancer等信息，更能知道高级结构是如何影响启动子、增强子从而影响基因表达的。
+组蛋白修饰： ATAC-seq可以预测出一个位点的开放程度以及可能有某种转录因子的结合，但是它不能知道这个因子是促进基因表达，还是激抑制基因表达的，有时候只通过基因层面上的鉴定来判断转录因子对基因的促进或者是不够的，它只是一种相关性。而这个时候如果你能提供像H3K27ac这类激活型组蛋白、H3K27me3这类抑制型组蛋白将能使你的数据更加的solid。
我们从国内较早一批研究iPSCs的学者，如裴端卿的工作可以看到，他们在17年和18年分别在CSC上报导了两篇关于成纤维细胞被诱导成iPSCs机制的研究，虽然两篇诱导的手段不同，但是可以看到，在解析iPSCs重编程中的染色质可及性的时候，不仅用到了ATAC-seq来描述细胞的身份转变，还通过H3K27ac来指征该区域的激活。其中一篇还通过调控成纤维细胞关键基因启动子区去乙酰化修饰，来达到了促进重编程的进程。
--scATAC-seq+scRNA-seq RNA-seq现在逐渐不能满足研究者的胃口，单细胞RNA-seq越来越常见，而ATAC-seq也开始出现了scATAC-seq来检测每个细胞的染色体开放情况。而有甚者，现在更前沿的技术可以一个细胞里同时进行RNA-seq和ATAC-seq，并且是单细胞水平的检测。SHARE-seq，能够实现在单细胞中同时高质量，高通量的检测基因表达和染色质可及性。该技术可以使用染色质潜力算法（chromatin potential），用ATAC和RNA的差异来预测细胞的对变化的方向。相对于以往仅依赖于RNA的预测手段，染色质潜力能够大大提前预测的时间。
数据上游分析
我想现在我已经说服你相信ATAC-seq是有用的，不仅能为你的工作增彩，最主要的是能帮你发现新的东西。那接下来，就一起学习下入门ATAC-seq分析所要掌握的知识点吧。
感谢下列文章，给出的宝贵流程和参考：
ATAC-seq data analysis: from FASTQ to peaks
ATAC-seq Data Standards and Processing Pipeline in ENCODE
ATAC-seq数据分析实战
Harvard FAS Informatics - ATAC-seq Guidelines
第一篇文章是我学习ATAC-seq的首选文章。它非常耐心细致地讲解了从raw data到ATAC-seq 的peak数据的分析流程，非常建议读一读。
第二篇文章是ENCODE官网的文章，比较晦涩一些，但里边不仅有标准的流程，还囊括了公认的质控指标，以及中间输入和输出文件的解释。
第三篇是Jimmy老师讲解的手把手ATAC-seq教程，他在生信知识分享这块我非常钦佩，也是国内做的知名度比较高的。
在分享之前的我们先看看大致流程是怎么样的？要理解这个流程其实不难，简单地将ATAC-seq理解成一个黑箱，你输入什么，然后得到什么？网上对于基本的流程讲的都很充分，这里简单地放个流程图。
1. 碱基质控与mapping
从公司得到fq文件后，初始的步骤其实与RNA-seq大差不差，都是得到bam文件。我一般就是走fastqc--trim_golare--bowtie2的流程。但ATAC-seq的mapping 记得带上这个参数--very-sensitive -X 2000。
2. 序列筛选
有几点是ATAC-seq需要特别注意的，序列筛选还是对bam文件或者sam文件进行处理。
去除线粒体和叶绿体基因，由于这些基因上没有组蛋白的结合，通常不在研究的范围，而其的存在会影响整体的分布，需要去掉。 去除Encode上规定的blacklist区域，有些区域是已公布的信号异常区域，需要去除。 X、Y染色体区域，这个因人而异吧，不想研究性别的也可以去除。 去重，Picard里的MarkDuplicates是常用的软件，一般会先标记上重复的reads。
ref: 讨厌又迷人的reads去重复
通过对这些区域的标记与识别，我们可以用bedtools以及samtools的常用工具将这些区域去掉。那么去重之后，就是我们称之为合适的比对序列了，可以进行后续的质检分析。
3. Shifting reads
由于Tn5酶是以二聚体的方面结合到染色体上的，其跨度大致是9bp，在第一篇ATAC-seq出来的时候，作者就考虑到了这个问题，在分析的时候，需要回补这个9个bp的碱基差。具体做法就是将正链正向移动4bp，将负链负向移动5个bp。我一般用alignmentSieve 一步到位。注意，不做reads shift 对单碱基分辨高的分析会有影响，例如TF motif footprinting。
4. Peak calling
Peak calling做的一件事情就是在统计学上判断真实的peak，因为Tn5在染色体上结合是个概率事件，如w何判断这个位置的reads足够为一个peak，这就需要用到统计检测，用的最多也是用来做chip-seq peak calling的软件macs。注意：你需要设置好你认为是Peak的标准。
p-value。 建立模型的方式，这会影响macs是如何判断peak的标准，--nomodel的意思是让其不要建立双峰模型来使两个“相邻”的峰shift成一个峰，而是向外shift(也就是在nomodel后要加上--shift -75 --extsize 150的参数)这个在前面提到过。
ref:Chip-seq处理流程
5. 峰文件bw的生成与组间标准化
网上许多教程没有讲的一件事情，如果你想要做不同组间样品的ATAC信号的比较，那么做normalize是必要的，注意，除了并不是bam to bw时单个样本自身将ATAC信号进行RPKM标准化外，还应进行组间标准化，即用haystack中的haystack_hotspots来做这件事情，输入的是bw文件；你可以选择输出峰的bin size，越小，峰会越shark。同时也要注意排除blacklist；这个包是要调用python3.6的。
6. 流程文件类型
除此之外，你还需要了解一下分析的大体流程以及中间产生的文件。我们从最大的一个人类肿瘤的ATAC分析的文章中来学习一下，它将TCGA中的23种肿瘤的ATAC都测了一遍，是一个很好resource。
ref: The chromatin accessibility landscape of primary human cancers
ATAC-seq上游处理之后，有两个最基本的文件类型——bed、bw.（chip-seq也是。）上述两者文件都可以导入到IGV展示peak的情况，但它们的功能有些不同。
Bed文件是Callpeak过后的峰位置文件，Bed文件每行至少包括chrom，chromStart，chromEnd三列；另外还可以添加额外的9列，这些列的顺序是固定的。但通常bed文件是没有表头的。我们一般用bed文件来定义特定的峰区域，就像你选择特定基因一样。 另外用macs call peak后会有peak文件具体来说有些许不同，narrowPeak其实就是峰的bed格式的文件，可以直接用， summit_bed指的是峰的中心的那个点。
Ref:++BED 文件格式++
bw文件是用于方便可视化peak的文件，因为上游处理完的bam文件通常都较大，不方便于快速展示，而将其转变成bw(bigwig)或者wig就会方便的多，而bigWig文件的显示性能又相较wig文件快得多，故bw是更常用的。而相较于bed文件相说，它不只提供了peak的位置，还有peak的高低。
质量控制
每步的质控对于得到正确的分析结果至关重要。Encode推荐了几个重要的指标。
image
比对率
通常要求在95%以上，但80%也是可接受的。
插入片段长度统计
插入片段长度是重要的评估实验好坏的指标，统计出的插入片段长度应该符合实验预期的长度。
FRiP(Fraction of reads in peaks)
peaks中的reads与总reads的比例，即文库中结合位点片段占背景reads的比例，可理解为“信噪比”，也是样本富集效果的评价指标，可一定程度上反应富集效果。通常要求大于0.3，大于0.2也是可以接受的。
库的复杂程度
它可以用这几个参数来描述：NRF、PCB1、PBC2，与reads 结合的独特性有关。通常要求，NRF>0.9, PBC1>0.9, and PBC2>3。
重复性鉴定
可以从两个方面出发，一个是Bam文件的重复性，一个是peak的重复性。其中，我喜欢用deeptools时的plotCorrelation；peak的话，用IDR做。如果生物学重复之间bam文件的相关性还可以，那我一般会用samtools merge 将其merge到一起，也可以用IDR选取重复样本中重叠的peak来做后续的分析。
常规下游分析
相信，你对ATAC-seq上游已经有了一个基本的认识，那么在质控没有问题后就可以开始下游分析了，总的来说，下游分析其实就是围绕着bw、bed文件展开的。
TSS等位点的peak plot展示
输入你想看的peak区域（bed文件），以及特定样本的bw文件，你就可以画出你的这些样本在这些区域的peak打开情况的热图。这里的需要用的到命令有两个：computeMatrix用来构建矩阵、plotHeatmap或plotProfile把矩阵展示出来。
这里有一些细节需要注意一下：
这两个命令的是包deeptools-3.2.1里的命令，所以你需要下载安装这个包。又因为这个包需要调用python，也需要安装python-3.6。 在画图的时候，你需要注意，要去掉一下不表达的峰区域，使图更加好看。 定义好你需要看区域中心，比如TSS位点、峰的中心，展示一个区域等等，这会影响你用computerMatrix的一些参数；你也需要定义好你想看的范围。
ref:快速使用 Deeptools 对 ChIP-seq 数据画图！
Motif 分析
homer中findMotifsGenome.pl可以做这个事情，注意其是用perl编写的软件，因此要安装perl-5.24，再用的时候记得调用一下。需要注意的是其有输入前要先将narrowPeak转变成特定格式的tmp文件才能读入。
ref: ATAC-seq分析实操生信技能树健明教程
Peak 对应的基因注释
这个一般将bed输入Great网站做就好，不同于传统的费舍尔精确检验，功能富集分析的p值是基于二项分布的计算得到的。需要注意的是，你输入的peak最好是有名字的，就是说，bed文件需要加一列名字列，不然的话，就难以找到基因与peak的对应关系。你还需要注意你关注的peak对应的基因的范围，一般看远端的peak的话，TSS上游5kb，下游1kb就可以。有时候国内的网有时候不太友好，可能要翻一下。
ref: 使用GREAT对peak进行功能注释
Peak 对应的区域注释
用homer中的annotatePeaks.pl做就好了，输入目的区域的bed文件即可。
Peak 的操作
对peak的上述操作，网上应该都讲解的比较详尽，但是如果你只想研究某些区域，想对这些区域进行选取的时候怎么办呢？这个时候操作bed文件的bedtools就是一个很重要的工具，它可以实现对peak进行取交集、取并集等等操作。
ref：Bedtools使用教程
进阶分析
除一些常规的分析，生信人员经常也会接到一些个性化的分析要求。这里也列举的常见的几种分析方式和策略。
peak 分类
2019年发表上CSC上的文章里展示了peak的方式：A图寻找差异peak的分类方式看着是k-means聚类来做的，实则不是，k-means是随机按三个群来划分peak的，无法做到如此。而具体按两组样本之间的差异peak分类，可以通过bam to count思路来做，用bedtools multicov后再结合deseq2来做；B图是将peak进行promoter及enhancer的分类。
PCA
主成分分析是看样本分布的情况。一个是可以用deeptools里的plotPCA做，还可以将bam转变成count，再用deseq2做。选定特定基因对应的区域画图
根据某基因list对应的区域来做图，它与直接将peak进行分类不同，它是先将基因分好类，再找对应的peak区。通常要用Great来反向找到基因对应的peak区域。而在进行基因或者peak分类时，组别在三者及以上时，可以采用差异基因的两两组合或者差异peak的两两组合，来实现peak的独特分类方式。如图，我将三个样本基因表达模型分成了大致三个类型。后续可以将这些基因对应的peak的热图画出来，以观察其是有无ATAC信号不同的分布特点。
五百君这次的ATAC-seq知识到这里就分享完了，希望对你的ATAC-seq学习有所帮助。
如果你觉得图文知识不过瘾，在b站还有视频教学哦！