ATAC-seq差异分析方法怎么选？
文章标题：Comparison of differential accessibility analysis strategies for ATAC-seq data
发表期刊：Scientific Reports
影响因子：3.998 (2019)
研究背景（冗长内容可跳过，直奔结果）
哺乳动物基因组中的基因调控涉及到不同类型的调控原件，诸如启动子、增强子和绝缘子。据估计在人类和小鼠基因组中有超过两百万调控原件，并且这些调控原件涉及不同的表观修饰过程，可以调控细胞类型特异和发育特定阶段的基因表达。活性调控原件必须维持在一个活性状态，这样可以和不同的转录因子结合来激活或者沉默靶基因。
ATAC-seq（assay for transposase-accessible chromatin followed by sequencing），是近几年开发出来用于检测染色质可及性的技术。相较于其他技术诸如DNase-seq、Mnase-seq和FAIRE-seq，ATAC-seq在不同的组织和细胞类型上实验操作相对简单。不仅如此，ATAC-seq实验可以允许超低细胞数，甚至低到单细胞水平。这些优势使得ATAC-seq成为各大基因组学学会（包括ENCODE、TCGA、PsychENCODE、IHEC和TaRGET II）最广泛使用的来定义开放染色质的技术。
在ATAC-seq中用来确定开放染色质区域的peak-calling分析通常是由ChIP-seq分析改良而来的。然而，ATAC-seq和ChIP-seq有着根本的差异，即ATAC-seq实验没有control或者input样本。不仅如此，peak caller，像是macs2，通过判断局部环境和基因组背景来判断开放染色质区域。在peak calling之后，多个样本的开放染色质区域先被合并，然后通过估计开放染色质区域上reads的差异来进行差异分析。由于我们通常假定ATAC-seq数据的差异分析与RNA-seq很相似，所以诸如edgeR和DESeq2等用来分析差异基因表达的软件被广泛的用来进行ATAC-seq数据的差异分析。
大多数开放染色质区域在两个条件下是相同的，只有一小部分开放染色质区域有可以被鉴定出的显著差异。开放染色质区域上的reads分布遵循一定的分布，即负二项分布。但是我们需要对于这些广泛使用的原本用来计算基因差异表达的工具用在ATAC-seq数据上的敏感性和特异性的综合比较，这样可以给用户提供指导。
在这项研究中，比较了4种广泛使用的软件包（DESeq、DESeq2、edgeR和limma）和两种经典的统计学方法（Wilcoxon秩和检验和t检验）以检测其敏感性和特异性；分析了不同样本量和测序深度对分析方法性能的影响；研究了不同P值和fold-change对判定差异可及区域（DAR）的影响；还研究了通过移除不需要的变量方法来矫正批次效应能极大提高ATAC-seq分析的敏感性。
结果
一、6种差异分析方法在模拟ATAC-seq数据上性能的比较
A. 作者比较了RNA-seq数据和ATAC-seq数据的信号分布。
数据来源是公开数据，同一个样本的RNA-seq数据和ATAC-seq数据。发现RNA-seq数据和ATAC-seq数据的信号分布模式不太一样。
B. 作者基于真实信号分布构建的模拟数据集。
设计了在1,5,10CPM信号强度上（分别代表低、中、高信号OCR）分别有10%，20%，50%和100%平均差的比较组。（图展示的是5CPM信号强度）
C. 作者在3种条件下应用6种方法的结果：
- 1. Wilcoxon秩和检验和t检验在所有信号水平的敏感性都比较低；
- 2. limma在所有信号水平都有最高的敏感性；
- 3. 在5和10 CPM组中，edgeR的性能都堪比limma，然后是DESeq和DESeq2。
D.作者对6种方法计算了假阳性率。
- 1. 所有方法的假阳性率能小于5%，特别是Wilcoxon秩和检验，在所有条件中都找不到任何假阳性结果；
- 2. 相比于limma和edgeR，DESeq和DESeq2对假阳性率有更好的控制。DESeq2有最低的FPR暗示其在所有方法中有最好的特异性。
E. ROC曲线分析显示，在高信号组（5,10CPM）中DESeq2，edgeR和limma有高敏感性和低FPR，但在低信号组（1CPM）中，edgeR和DESeq2敏感性下降严重。
二、差异分析中样本数和测序深度对敏感性的影响
A. 作者将每一边样本重复数范围从2扩大到20，从而研究样本重复数对每个分析方法敏感性和特异性评价的影响。
- 1. 随着样本数量增加，除了DESeq2之外，所有方法的敏感性都在增加；
- 2. 相比于limma，Wilcoxon秩和检验和t检验，另外三种基于负二项分布的方法在1 CPM组（20%平均差）有较低的敏感性。而前三种方法在重复数达到15之后显示出几乎100%的敏感性；
- 3. Wilcoxon秩和检验在重复数少于6时几乎不能确定真正类（True positive，TP），而t检验需要至少4重复；
- 4. 在重复数有限的情况下，limma在1 CPM（50%平均差）、edgeR在10CPM（20%和50%平均差）的情况下有最好敏感性。
B. 作者研究了提升样本重复数对假阳性的影响。
- 1. 所有的方法都将FPR控制在5%以下；
- 2. 由于TP的召回率由于样本数而提升，DESeq、DESeq2和edgeR在1CPM组中的FPR保持在0左右，而TP偏低；
- 3. 在两重复情况下，edgeR在1CPM（10%，20%，50%和100%平均差）能有25%的TP，而在其他两组中edgeR能有40%的TP而FPR只上升了3%；
- 4. DESeq和DESeq2对FPR控制的最好，但TP的召回率偏低。
C. edgeR、DESeq2和DESeq都是基于负二项分布的，测序深度对他们的结果影响较大。Wilcoxon秩和检验、t检验和limma则不会。作者模拟了3种测序深度（10M、20M和30M有效reads数）在6重复条件下来检测6种方法的性能。
- 1. 测序深度的增加能极大优化edgeR、DESeq2和DESeq的敏感性，特别是在1CPM组中；
- 2. Wilcoxon秩和检验、t检验和limma在20%和50%平均差的条件下，测序深度从10M提升到20M时，敏感性也增加了。这可能是由于测序深度带来的更准确的CPM计算引起的。
三、在真实ATAC-seq数据上用DESeq2，limma和edgeR进行差异分析的性能评价
模拟数据显示DESeq和DESeq2有更好的特异性，edgeR有更好的敏感性。
A. 作者使用发表过的小鼠肝和肾的ATAC-seq数据来比较分析方法。
B. 作者比较了用edgeR、limma和DESeq2得出的组织特异的DAR。（DESeq2 and limma：padj < 0.01；edgeR：FDR < 0.01）显示出3种方法在确定组织特异的DAR有比较高的一致性。
C. 作者用down-sampling法来评价三种方法在真实数据上的性能。作者将在每一步中新发现而在原来6×6比较中没出现的DAR称作潜在假正类。
- 1. edgeR在低重复数上有最高的敏感性；
- 2. DESeq2在所有down-sampling步骤中出现的潜在假正类控制的最好。
D&E. 作者还分析了q值和fold-change阈值设定对分析结果的影响。
F. 基因表达及其pomoter区对应DAR的相关性分析。
四、通过移除不需要的变量的方法来改进ATAC-seq差异分析
批次效应去除方法——RUVSeq（去除不需要的变化，R包）。
由于这个包经常使用在RNA-seq数据上，作者测试了RUVSeq包在两个ATAC-seq数据差异分析上的性能。
结论
- 当需要高敏感性时或者样本数有限时，作者推荐使用edgeR来找DAR。
- 当在大样本数中需要特异性时，作者推荐DESeq2来找DAR。
- 为了确保充分的敏感性，作者推荐每个条件至少三个重复。
- 改善差异分析的敏感性需要充分的测序深度，特别是再低信号区域，如enhancers。
- 用PCA去检查样本分布是很重要的。
- RUV策略能修正数据的批次效应，从而改进敏感性。
文章转自微信公众号：嘉因生物
