前言背景知识关于ATAC-seq发展现状关于染色质可及性测序的方法小结MNase-seqDNase-seqFAIRE-seqATAC-seq小结关于Tn5转座酶原理ChIP-seq中阴性对照的设置pioneer factor关于ATAC-seq数据分析shift-extend方法预测peaks的原理第一部分——pre-analysis比对前质控比对比对后质控小结第二部分——peak calling小结第三部分——高级分析PeaksPeak differential analysisPeak annotationMotifsMotif database and scanMotif enrichment and activity analysisFootprintsDe novo toolsMotif-centric tools对于footprint分析的评价Nucleosome positioning第四部分——多组学数据联合分析建立调控网络结构与ChIP-seq联合分析与RNA-seq联合分析建立调控网络未来展望和总结生信教程补充后记
今天分享一篇文献，主要介绍ATAC-seq中的一些主要问题。
文献标题：From reads to insight: a hitchhikers guide to atAC-seq data analysis
文献地址：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-1929-3
下载地址：https://genomebiology.biomedcentral.com/track/pdf/10.1186/s13059-020-1929-3
在阅读这篇文章的过程中，前面因为我自己跑过流程，所以看起来还是比较轻松，但是到了后面的高级分析部分，就有些难度了，很多之前模模糊糊的地方都要再去查找文献去找到答案，这里放一个我觉得收获许多背景知识的中文博士论文：
ATAC-seq数据分析软件开发及其在肥胖诱导的慢性炎症研究中的应用，作者：左祖奇
因为在知网可以下载，但是没有账号你们可能还是没法下载，所以我把它也放到了百度网盘里，有需要的在公众号回复 “ATAC” 拿到我下好的pdf文件吧。
DNA序列包装成核小体→染色质→染色体
因为人体基因组是高度压缩状态，而转录和翻译都是需要在松散结构下的染色质情况下才可以进行，所以认为染色质的开放程度和基因转录密切相关。
关于ATAC-seq的原理，我之前写过推文介绍：
不过又来又学到了更多知识，只不过一直没有补充进入，今天发现了一个写的还不错的推文，在这里推荐给大家去学习相关的基础知识：
近些年来ATAC-seq技术的发展非常的迅速，从2013年greenleaf发表的第一篇关于ATAC-seq的文章后，这项技术迅速得到大家的喜爱：
包括在去年，greenleaf与10xGenomics公司和做，开发了10XscATAC-seq的测序方法，各种新技术层出不穷。
但是针对ATAC-seq的数据分析工具不多，因为ChIP-seq数据和ATAC-seq数据的相似性，目前主要使用的都是以前开发用于ChIP-seq的工具，默认为这两种数据分布结构是相似的，但是并没有人真正系统地去评估这两种数据分布。
今天分享的这篇综述主要是介绍在ATAC-seq分析过程中的一些思路和套路流程。主要包括4个方面：
pre-analysis：质控和比对
peak calling
高级分析：例如peak的注释，找motif，找核小体，找TF的印记（footprints）等方面
多组学数据联合分析
目前用于研究染色质可及性的方法主要有以下四种：MNase-seq、DNase-seq、FAIRE-seq和ATAC-seq：MNase-seq是通过对核小体保护的DNA片段测序，从而间接反映染色质可及性的方法，其他三种均为对检测染色质上的开放区域测序，直接反映染色质的可及性。
微球菌核酸酶（ Micrococcal nuclease， MNase）是来源于金黄色葡萄球菌分泌的一种核酸酶，同时具备核酸内切酶和外切酶的活性。MNase优先对裸露的DNA或核小体之间起连接作用的DNA进行切割和消化。所以这种方法一般用于检测开放区域。
标准的 MNase-seq的流程主要用于对核小体片段（~150bp）或更长的片段进行测序。
脱氧核糖核酸酶I（ DNase I）是由人的基因 DNASEI编码的核酸内切酶，可以非特异性的对双链DNA进行切割。但是没有外切酶活性。
在基因组学和染色质的研究中DNase I敏感的位点被认为是开放的，可接近的染色质的特征。低浓度的 DNase i对基因组上非核小体占据的的开放区域进行切割，这些区域被称为DNase I敏感位点（ DNase I hypersensitive sites，DHSs）。
DNase-seq目前已成为检测染色质可及性的“金标准”。
DHSs中序列特异性的转录因子的结合也会阻止 DNase I对DNA的切割，从而可以在单个碱基水平观察到转录因子的占据情况，即转录因子的印迹分析（footprint）。转录因子的印迹分析已被用于发掘人类细胞中细胞特异性的转录因子结合的“基序”（motif），并揭示了转录因子结合与染色质结构、基因表达和细胞分化的相关性。
甲醛辅助的调控元件的分离（ Formaldehyde-Assisted Isolation of Regulatory Elements，FARE）是一种直接检测无核小体占据的DNA序列的方法。其原理是，缠绕有DNA的核小体和无核小体结合的DNA，在苯酚和氯仿中的溶解度不同，因而在苯酚和氯仿形成的两相液体中呈现差异分布：缠绕有DNA的核小体分布于两相界面处，而无核小体结合的DNA分布于亲水相中。
具体的实验流程包括：
使用甲醛对染色质进行交联以便巩固蛋白质-DNA的结合。
通过超声处理，打断染色质，形成DNA片段或者DNA-核小体复合物。
通过苯酚-氯仿抽提纯化无核小体结合的DNA片段。
建库测序
FAIRE直接富集了活化染色质的区域，同时无核小体占据的区域并没有被酶解。
该方法已被用于真核生物细胞全基因组范围内的：
开放染色体区域检测
核小体位置确定
转录因子的印迹描绘
ATAC-seq建库过程简单快捷，只需要两步操作，同时仅需要较少的细胞数目，而且可以在很高的分辨率下揭示染色质的结构。
ATAC-seq仅仅使用500到50000个细胞就可以实现与 DNase-seq使用百万数量级的细胞才能达到的灵敏度和特异性。
ATAC-seq目前亟待解决的问题是对测序数据分析，原有的分析方法不适用于ATAC-seq的数据分析或仅可以有限度的使用。
http://www.bioon.com.cn/sub/showarticle.asp?newsid=72268
Tn5转座子是一种细菌转座子，最早由E. coli中发现，是一段含有若干抗性基因和编码转座酶基因的DNA片段。
其中IS50R和IS50L的序列高度同源，只是IS50L的一个碱基存在突变。
IS50具有19bp的倒置末端（外末端outside end，OE和内末端inside end，IE），两末端倒置有7个碱基不同。此倒置末端是转座酶（Tnp）的作用位点。
IS50L和IS50R均含有编码转座酶（TnP）以及转座阻遏蛋白（lnh）的基因，但由于IS50L中的碱基突变，造成翻译提前终止，所以仅有IS50R可以产生正常的有活性的TnP和lnh。
两个转座酶（Tnp）分子结合到Tn5转座子的OE末端，形成两个Tnp-OE复合体，随后两个复合体通过Tnp的C末端相互作用进行联会，形成一个Tn5转座复合体，此时Tnp产生切割DNA的活性。
随后Tnp利用切割活性，经过一系列化学反应切除供体DNA，离开供体链。
结合到靶DNA上时，Tn5转座复合体识别并攻击靶序列（Target site），将转座子插入到靶序列中，粘性末端通过DNA聚合酶、连接酶作用进行填补，两端形成9bp正向重复序列。整个转座过程完成了基因从原始DNA被剪切之后粘贴在另一受体DNA的过程，实现了基因的“跳跃”。
（解释1）
（解释2）
根据之前的报道，Tn5转座酶以同源二聚体的形式结合到DNA上，在两个Tn5分子间隔着9-bp的DNA序列。根据这个情况，每个Tn5同源二聚体的结合事件会产生2个「Insertions」，中间隔着9bp。因此，真实的"开放"位置的中心在Tn5二聚体的正中间，而不是Tn5的插入位置。为了尽可能的还原真实情况，我们对Tn5的「Insertions」进行了校正，即正链的插入结果往右移动4bp(+4 bp), 负链的插入结果往左偏移5bp(-5 bp)。
推荐看一篇推文：
简单归纳后要点如下：
为什么需要设置阴性对照：
因为超声破碎过程中DNA的断裂不均一，尤其是一些开放染色质区域在超声样本中优先累积，未经过IP的样本超声破碎后会产生数量不小的peaks。
可以有去除背景噪音的作用（排除因本身表达水平高或一些非特异性结合所造成的假阳性peaks）。
还可以根据Input中的靶序列的含量以及染色质沉淀中的靶序列的含量，按照取样比例换算出ChIP的效率（如果用同一引物进行PCR，ChIP组和input组亮度差不多，说明ChIP效率高，样本中所有的目的基因片段都被ChIP下来了）。
如何设置对照：
input对照：少了加抗体的步骤，在交联和超声裂解后，并没有加入抗体，跳过了免疫沉淀过程，直接将这些DNA进行测序。
Mock-ip对照：使用与目标蛋白无关的非目标抗体（IgG或者标签）进行“模拟”的IP，为了防止抗体的非特异性结合。
pioneer factor，先驱转录因子，是一种特殊的TF，可以结合在核小体DNA上，直接介导染色质可及性的改变。
有几点是我之前没有太注意到的，这里标注下。
ATAC-seq数据中包括了开放染色体区域检测（call peaks），核小体位置的检测以及转录因子的印迹（footprints）的分析。
由于真核细胞的转录起始位点被报道处于开放状态，因而可将TSS信号的强度作为检测 DNase-seq和ATAC-seq信噪比的标准，同时，在全基因组范围内，DNase-seq和ATAC-seq的所获得的片段长度分布，应可以见明显的核小体占据“峰”。
由于其建库过程中，可能引入线粒体DNA，因此需要检测其中线粒体DNA的比例。
染色质可及性分析的首要目的是寻找到基因组上的开放区域。所谓信号峰搜寻（call peaks），就是在全基因组范围内找出测序读长累积形成的脉冲峰的位置及信号强度。这些峰的位置代表了基因组上的开放区域，峰的髙度或面积代表了该区域的开放程度。同时，由于这些区域与基因的表达调控密切相关，他们与各自附近的基因在基因组上的位置关系和功能注释也值得关注。
信号峰的搜寻一般会产生存储有信号峰位置信息的BED格式的文件。对信峰上信号强度统计，可以获得每个开放区域的可及性。在具体的生物学比较分析中，研究者可通过比较不同组别之间信号峰的强度差异，或对差异信号进行聚类分析，筛选出感兴趣的开放区域。
染色质上的开放位点意味着没有核小体的占据，这些区域里往往包含大量的转录因子结合位点。转录因子与DNA的结合，也对结合区域提供了保护，避免了被 DNase或Tn5酶的剪切。转录因子在DNA上的结合区域很短，一般为8-30bp长，相比于结合区域周围，转录因子结合区域酶切的信号强度往往更弱，统计结合区域及附近区域遂个碱基上的酶切信号强度可以看到结合区域呈现明显的凹陷，该凹陷指示该转录因子在该区域的确发生了结合。染色质可及性分析的目的之一就是找到在开放区域上富集有哪些转录因子的结合位点，以及描述这些位点上转录因子的结合情况。
根据算法的不同，转录因子富集分析的方法主要有两大类：
一类是根据已知的转录因子在基因组上结合位点的偏好位点信息，在提供的DNA序列上进行比对，进而获得转录因子的富集情况，这类分析只能找到已知的转录因子
一类是对提供的DNA序列进行 de novo分析，这类分析有可能发现新的转录因子结合区域
对已知转录因子的搜寻，依赖于已知的转录因子结合的位点信息，这些住息来源于前人对转录因子所做的ChP-seq分析，且往往以位置权重矩阵position weight matrices，pwMs）的形式存储。专门存储转录因子PwMS的数据库主要有 MatBase， JASPAR， TRANSFAC和 UniPROBE。
文中提到了一个shift-extend方法来预测peaks，于是我找了很多资料，最后终于找到了一个不错的解释：
ChIP-seq 分析------原理：https://www.jianshu.com/p/dc493cb7b1b3
对于一个DNA序列来说（有正负链的），它mapping的位置正负链都有的，对这些reads位置进行统计画图可以看到一个红色的peak，一个蓝色的peak。这两个peak说明的是一个事情，就是这个地方有富集。最后对这两个peak进行merge，最后变成了一个富集区域。灰色的peak!
所谓的shift-extend，就是把PE测序片段进行延伸，然后这样就可以直接得到灰色的覆盖区域最多，peaks也就最高了。
通过看图示，发现具体过程如下：
先将片段向外移动s个单位，然后再向内延伸2s个单位。
这个是所有测序技术都需要进行的QC流程，主要是看看接头有没有去除干净，GC比例是否合格，测序质量情况如何等等，可以使用linux平台的工具如：FastQC
测序文件在5'开始时和3'端结束前，测序质量可能会有一个大幅降低，这个是可以接受的，原因和测序中酶活性以及机器设计原理相关。
一般在我们运行完去除接头这些比对前质控操作后，一般会再运行一次FastQC来查看质控效果。
然后就对通过质控fastq文件进行mapping。一般使用BWA或者Bowtie2工具。
对于比对完成的bam文件我们可以使用samtools或者经典的Picard工具来分析下比对后的情况，文献这里建议：
unique mapping rate > 80% 属于合格，会比较好。
对于哺乳动物来说，ATAC-seq结果中应该至少检测到 50 million的reads用于分析开放区和做差异分析。找到200 million的reads用于分析TF 的footprints
做完比对后，我们仍然需要做质控，从unique mapping reads/rates，duplicated read percentages，fragment size distribution等等方面去评估。
如果遇到下面的情况，则reads需要被去除：
测序数据质量太差
PE测序中无法找到overlap区域
线粒体基因（因为线粒体基因都是可接近的，所以ATAC-seq数据中会有较多线粒体的序列，这部分需要需要我们去除）
ENCODE数据库中记录的需要被去除的区域
PCR重复序列
上述的质控可以提高数据的准确性，降低假阳性率。
但是，上述的可能并不是ATAC-seq质控的特异指标，因为我们在做RNA-seq或者ChIP-seq数据分析中也会做这些质控，下面有几个ATAC-seq质控的特异指标：
片段长度分布图如下：
这个图如何看可以看我之前写的推文：
对于nucleosome-free regions（NFR）片段 < 100 bp，这些片段理论上是富集在不同基因的TSS区域的。
对于核小体结合区域片段，在TSS区域应该存在一个低谷，而在TSS两侧翼区则应该存在富集现象。
如下：
上述这种评估可以使用ATACseqQC工具进行评估。
还有一个比较容易忽视的细节，也是比较重要的，根据一个朋友的经验，他提醒我关于reads需要移动的一个tips——"+4和-5"规则：
正链+4bp
负链-5bp
作者提供了他们自己分析所使用到的工具pipeline：
FastQC➔ trimmomatic➔BWA-MEM➔ATACseqQC
目前存在的一些用于peak calling的工具以及他们背后用到的统计学分布总结：
一般call peaks都是使用MACS2工具，在ENCODE官网上的ATAC-seq的pipeline上也推荐使用MACS2来call peaks。
目前存在的call peaks工具99%都是和ChIP-seq分析混用的，只有一个工具是专门针对ATAC-seq数据而开发的——HMMRATAC。
我们知道，在做ChIP-seq时，需要有正常input control对照。但是在ATAC-seq中一般不设置input control对照。
对于ATAC-seq的PE数据，在经过比对后，得到的范围涉及到了NFR和核小体结合区：
而对于所谓的开放区，其实是来自NFR区域的比对结果的，或者使用一种生信方法——shift-extend。
目前call peaks的工具大致分为2大类：①基于计数 ②基于分析形状
基于计数的call peaks使用不同的统计方法来比较某个特定区域内的reads分布和随机情况下的reads分布形状。如MACS2、HOMER、SICER/epic2都是假设是泊松分布；而ZINBA则假定是零膨胀负二项分布；等等。
因为F-seq和ZINBA 并不是经常有人维护，所以作者不建议使用。
总的来说，基于计数的call peaks使用更多，更容易解释。
基于分析形状的call peaks不常使用。
HMMRATAC是唯一一个专门针对ATAC-seq的call peaks工具。优点在于：结果比MACS2和Fseq等工具找到更好，并且可以同时提供给我们核小体的位置信息。缺点在于：计算量非常大，用到很多机器学习方面的算法，如三态半监督隐马尔可夫模型（一听就很厉害，让人不想去看，哈哈哈哈哈哈）
关于实验中的设置生物学重复问题：设置生物学重复可以减低假阳性，同时提高可重复性。大多数工具在使用时都可以通过参数设置处理生物学重复。
放一个真实数据的ATAC-seq图：
分成3大部分：HMM based工具、count-based工具以及shape-based工具。
RUNX1 motif track：是从JASPAR公共数据库里得到的RUNX1 footprint结果
K562 ChIP-seq track ：是从ENCODE公共数据库里得到的RUNX1的ChIP-seq结果。（相当于是一个标准答案，因为是直接用ChIP拉下来的片段）
count-based工具之间结果相似，但是和shape-based工具之间结果有出入。
这些数据证明，可以通过peaks的结果找到已知转录因子的motifs。反过来，也说明这些转录因子与染色质的可及性密切相关。
目前没有工作去比较ATAC-seq中peak calling的表现性能，作者推荐使用MACS2和HOMER这种工具来peak calling，如果服务器足够强大，就推荐使用HMMRATAC去call peaks。
目前存在的一些Peak differential analysis的工具总结：
目前没有专门针对ATAC-seq数据找差异peaks的工具。
目前存在的工具一般分成2种：
consensus peak-based工具：如HOMER，DBChIP，DiffBind。他们内部其实主要依赖于差异基因分析的统计学原理，例如我们常用的edgeR、DESeq2等。我们很清楚这些工具的前提是他们假设数据分布是负二项分布（NB）。
consensus peak是指：不同生物学样本重复得到的peaks进行合并后，找到的一些所有重复样本中都存在的peaks。这样可以减少假阳性结果。HOMER默认会用将所有样本的peaks混合pool在一起后算出consensus peak。而DBChIP，DiffBind则通过在不同样本之间取交集来得到consensus peak。
Sliding window-based工具：使用这种方法的话，无须去生成consensus peak，他们会评估所有全基因组上每个bin区域，这样当然会有更多的假阳性情况，于是需要做FDR进行校正p值。
独立工具有：PePr和DiffReps，他们使用负二项分布检验、G检验或卡方检验。ChIPDiff则应用HMM去计算2个临近窗口的相关性。
还有些工具，如csaw，则是依赖于其他的DE分析的R包edgeR。
Sliding window-based工具因为是对整个基因组范围进行分析，所以是一种unbiased的方法，但是这也提示，使用这种方法进行分析时，需要一个很严格的过滤标准！
目前绝大多数的研究都证实ATAC-seq数据中reads的分布符合泊松分布，这和RNA-seq数据的分布是一致的。
shape-based的差异peaks分析工具没有专门针对ATAC-seq数据的，但是因为shape-based的方法可以用到同一个数据的2个维度——reads+分布形状，所以作者认为这种方法应该会给我们提供更多的信息。并且作者推荐使用csaw，因为这个工具核心是依赖于edgeR，这样结果更好解释。
使用工具：HOMER，ChIPseeker和ChIPpeakAnno使用最多。
一般对于peaks的注释，都是找离peaks距离最近的gene（内含子或外显子）或者调节元件（promoter，5′ UTR，3′ UTR等）。
简单来说，就是exon与intron是互斥的，exon包含UTR和CDS。具体的关系可以看我之前分析CCDS文件的系列过程：
一般来说生信分析的软件，对于peaks注释后的可视化展示最经典的如下饼图：
代码的话其实直接用Y叔的R包就可以了，非常的简单，有空了把相关画图代码找找放上来吧~
得到这些注释结果后，可以通过GO/KEGG/Reactome等数据库进行富集分析，其实代码也非常的简单。日后用到再整理吧~
所谓的Motifs就是那些可以结合TF的DNA序列，而TF结合的位置称为TFBS（TF binding sites）。TF如果想要对转录进行调控，就必须和DNA上的顺式作用元件结合，而TF结合的前提通常来说是这段DNA序列是可接近的，也就是ATAC-seq中可以测到。不过也有少数TF可以和那些非开放区域进行结合。
TF调节转录的机制：1）竞争组蛋⽩或⾮组蛋⽩；2）co-factor互作
具体介绍TF调节转录机制的综述推荐：
文章地址：https://www.nature.com/articles/s41576-018-0089-8
下载地址：https://www.nature.com/articles/s41576-018-0089-8.pdf
和其他的生信分析思路类似，想要探索Motifs的信息，我们需要先有一个Motifs的数据库，例如：
JASPAR（包括多个物种，可以通过API进行访问、存在Bioconductor的R包）
CIS-BP和TRANSFAC（真核生物）
HOCOMOCO（人和鼠的数据）
RegulonDB（大肠杆菌数据库）
HOMER（Linux）以及Bioconductor上的一些R包（TFBSTools和motifmatchr）都可以通过检索给定的DNA序列来判断TFBSs。
PWMScan可以直接在线使用，需要提供bowtie index后的文件。
MEME工具（之前使用过）包括几个组件组成，其中：
FIMO去搜索Motif
MAST合并多个Motif
MCAST推断调节模块
推荐使用的工具——MEME和PWMScan
前面的工具，找到了Motif的位置和发生频率，接下来就看这些Motif在peaks里的频率了。HOMER用到的统计学原理是超几何检验，这个比较好理解，不懂的可以去看我之前写过关于超几何检验的推文。MEME-AME则用到的是秩和检验 。DAStk通过计算MD分数来判断。当然还有更多其他的方法，不过这些方法都是通过应用不同的统计学方法来比较Motif在peaks中的频率，从而得到真正的Motif吧。
通过计数fragments的读数，可以得到TFBS的可及性，这个和TF的活性相关。而ChromVAR工具就是针对scATAC-seq数据而设计的，但是是否可以 应用于bulk ATAC-seq数据目前没有研究。DiffTF则针对所有TFBS计算一个可及性改变FC
这里提到的所有工具都是用来间接预测peaks区域内的TFBSs。但是这里找到的TFBSs可能有一些是错误的。因为目前并不是所有的TF都有明确的Motif序列，而且，来自相同家族的TF可以有共同的Motif结合序列。
我们解析TF调节也可以用Footprints。所谓的Footprints是指：激活的TF结合的DNA序列，这段序列因为和TF结合而不受Tn5酶切处理。
使用Footprints进行分析存在的几个要点：
前面提到关于read的移动，需要正链的插入结果往右移动4bp(+4 bp), 负链的插入结果往左偏移5bp(-5 bp)。
因为Tn5酶具有偏好性，所以对于一些短暂结合的TF，Footprints的检测存在困难。在过去用DNase-seq时这种困难也存在。
分析Footprints的工具主要分成2大类：
de novo：根据Footprints的典型特征，预测所有peaks区域的Footprints，得到的结果和已知的Motif去做匹配，少数无法匹配到的则为新发现的Motif。
motif-centric：需要我们提供一个TFBSs的信息文件，并通过机器学习的思想去区分这些属于结合状态和非结合状态。
这种方法的一个重点在于要用数学方法上去定义什么是一个footprint，并且尽量降低由于Tn5酶切偏好性引起的footprint噪音。
这里列举HINT-ATAC工具，使用隐马尔科夫模型（HMM），同时矫正了Tn5酶切偏好性：
HINT-ATAC找到的footprint，同样在K562细胞系中，用RUNX1的ChIP-seq结果中得到了验证。
对于使用HMM的工具，本质上都是需要监督学习的，所以也就是说，需要我们手动去注释一些基因区域，因此这类工具在更大范围内的使用问题仍需测评。
利用非监督学习的方法进行聚类，基于一系列参数如：距离TSS距离，PWM分数，reads分布，reads数目等等将可能的TFBSs分成结合状态和非结合状态。CENTIPEDE工具对于参数的变化比较敏感；msCentipde可以提高低深度和低质量数据的表现；PIQ在有生物重复的情况下可以提高结果的鲁棒性。
利用监督学习的方法进行聚类，基于高质量的ChIP-seq数据来注释真正的TFBSs。MILLIPEDE和BinDNase使用逻辑回归，DeFCoM使用支持向量机SVM，BPAC使用随机森林进行鉴定。
一般来说，监督学习工具会比非监督学习工具和de novo工具效果更好，但是其通用性就稍逊一筹。
偏差校正在DNase-seq和ATAC-seq足迹检测中都很重要。
能够有效实现足迹分析的ATAC-seq最小测序深度是多少，目前没有标准。虽然建议每个样品read数超过2亿，但有报道称DeFCoM对于更少的测序read数也能有不粗的表现。
对于低质量和之前未报道过的motif，de novo方法仍然具有优势。
作者认为HINT-ATAC可以是一个不错的选择，因为它具有ATAC-seq特异性的偏好校正。
同样的，研究人员可以结合多种工具的结果来获得高度可靠的足迹。
在ATAC-seq数据中，较长的reads片段对应着开放区中缠绕核小体的DNA片段。有许多工具用来分析检测这些缠绕核小体的DNA序列，但是根据研究证明，由于ATAC-seq数据中这些区域的覆盖深度较浅，所以相比与MNase-seq数据来说，分析更加困难。
针对MNase-seq开发的软件如DANPOS2，PuFFIN，INPS，和NucTools，可以在ATAC-seq数据过滤得到核小体相关片段后使用，而NucleoATAC和HMMRATAC是专为ATAC-seq开发的。
所有这些工具都具有典型ATAC-seq实验的相同潜在缺点，即染色质开放区之外的覆盖率较低。期待未来开发用于ATAC-seq的生物信息学方法，以更有效和精确地捕获核小体的占位。目前作者认为HMMRATAC和NucleoATAC是用于ATAC-seq核小体检测的两个有用且特异性的工具。
由于开放染色质是大多数TF结合的前提条件，但是不是所有开放染色质都有TF的结合，因此ATAC-seq峰通常与TF ChIP-seq峰重叠，但通常更宽。联合TF ChIP-seq和ATAC-seq可以相互验证彼此的质量和可靠性。
ChIP-seq中存在TF的峰，而在ATAC-seq中不存在，可能指示了先驱转录因子（pioneer factor），它结合到封闭染色质，然后招募染色质重塑因子或其他转录因子并起始转录。
ATAC-seq也可以与标记组蛋白修饰的ChIP-seq联合分析，验证与活跃染色质标记（如H3K4me3的，H3K4me1，H3K27ac等）正相关，与不活跃的染色质标记（如H3K27me3）负相关 。
由于ATAC-seq实验方法的\简便性和样品需求较少，因此可以在做ChIP-seq实验之前，把ATAC-seq当成一种预实验。
我们可以验证差异基因在各自的TSS周围是否也具有明显的染色质可及性差异，从而从染色质可及性的角度验证自己的理论假说。
可以推定差异基因受到开放染色质中特定TF的调控。
关于增强子的作用机理：
可以看到，启动子promoter一般在target gene的上游，而Enhancer的话可以在gene的上/下游，或者在更远的位置。发挥作用时，只需增强子和TF结合后，通过增强TF的活性，促进Promoter的转录活性。
增强子在线性基因组中可能非常遥远，但在空间上接近其目标基因。这导致增强子的直接靶基因难以预测，因为很多研究都是直接把远端增强子联系到最近的基因上，而非真正的target gene。
对于scATAC-seq，Pliner等人推出了Cicero，可将增强子和启动子联系到同一靶基因。尽管已证明Cicero可以用于scATAC-seq，但是没有证据证明Cicero否适用于样本量小的bulk ATAC-seq数据。
ATAC-seq近年来发展迅速，在实验protocol取得了较大的进展，但生物信息学分析工具的进展缓慢，没有成熟的分析pipeline。
在整个分析过程中，比对到参考基因组和质控步骤与RNA-seq和ChIP-seq中类似。至于call peak，大多数ChIP-seq的工具都与ATAC-seq数据兼容，ATAC-seq特异性的call peak工具较少。
对于下游分析，peak差异分析可以概述染色质可及性的变化。为了推断生物学功能和相关的TF，peak注释和motif富集分析是初步了解的首选。
motif和footprint分别是调控事件的直接和间接指标。检测footprint的困难来自酶切偏倚和TF结合DNA时间短引起的信号微弱。
由于ATAC-seq数据固有的弱点（峰以外的区域read覆盖率很低），核小体检测仍然很困难。
作者建议的分析流程：
用FastQC，trimmomatic和BWA-MEM进行预分析
用MACS2进行peak calling
使用csaw进行peak差异分析
使用MEME进行motif检测和富集
使用ChIPseeker进行注释和可视化
使用HMMRATAC进行核小体检测
使用HINT-ATAC进行足迹分析
如果同时有RNA-seq数据，则可以使用PECA方法重建调控网络。
这篇综述内容还是很多的，也学到了许多关于ATAC-seq的新知识。
另外也找到一个生信分析的系列教程：
http://qiubio.com/new/book/chapter-05/#第四章-chip-seq数据分析chapter-4-chip-seq-data-analysis
有需求的同学可以自行学习~
另外放上我当初学习ChIP-seq的代码教程吧，因为当初是第一次学习，所以很多就会相对比较细致了，但是也可能会存在一些问题，大家自行判断学习：
近来开学在即，本着不信谣不传谣的党中央精神，我在家默默看完了这篇文献。希望大家也要在家坚持：不信谣不传谣~
联系客服
