整合全基因组测序控制DNA剂量提升ATAC-seq峰识别准确性

《NAR Genomics and Bioinformatics》：Controlling for DNA dosage with whole-genome sequencing improves ATAC-seq peak calling

　　

ATAC-seq（transposase-accessible chromatin using sequencing，转座酶可及染色质测序）是一种用于跨细胞类型绘制染色质可及性并鉴定顺式调控元件（cis-regulatory elements, CREs）的强大方法。研究人员揭示，ATAC-seq峰识别（peak calling）受到DNA剂量局部偏倚的影响，该偏倚源于拷贝数变异（copy number variations, CNVs）以及DNA复制时序（DNA replication timing, RT）。作为FANTOM联盟工作的一部分，研究人员提出通过整合全基因组测序（whole-genome sequencing, WGS）数据来增强基于模型的ChIP-seq分析（model-based analysis of ChIP-seq, MACS）流程，类似于ChIP-seq分析中的输入对照（input controls）。该调整增加了ATAC峰的数量和宽度，并改善其与转录起始位点（transcription start sites, TSSs）的邻近性。这些ATAC峰表现出经典的CRE表观遗传标记，并富集与性状和疾病相关的遗传变异。此外，与未使用WGS对照所识别的峰相比，使用WGS所识别的峰数量增强了与基因表达水平的相关性。总体而言，这些结果表明，整合WGS作为对照可显著提高ATAC-seq峰识别的准确性。关键的是，即使低深度WGS数据也足以改善峰识别，使该方法具有成本效益和可扩展性。研究人员将该方法实现为开源Nextflow流程，以提高可及性和可重复性。通过挑战均匀基因组可见性（uniform genomic visibility）假设，该方法对其他基于DNA测序的技术也具有广泛意义。

论文解读：整合全基因组测序控制DNA剂量改善ATAC-seq峰识别

研究背景与问题
ATAC-seq（transposase-accessible chromatin using sequencing，转座酶可及染色质测序）可全基因组探测染色质可及性，所需细胞少，适用于稀缺原代样本、异质组织及单细胞分析。通过选择性分析开放染色质，ATAC-seq可系统注释顺式调控元件（cis-regulatory elements, CREs），这些区域可招募转录因子（transcription factors, TFs）并承载活跃转录起始。CREs通常定义为ATAC-seq峰。MACS（model-based analysis of ChIP-seq，基于模型的ChIP-seq分析）是最常用峰识别程序，原本用于ChIP-seq，用于ATAC-seq需参数调整。MACS以滑动窗口识别相对基因组背景显著富集的区域，并计算期望读数λ。由于局部富集受多因素影响，MACS估计λ_local=max(λ_BG, λ_1k, λ_10k)，其中λ_BG为全基因组常数，λ_x来自以候选区域为中心的x kb窗口。ChIP-seq常有匹配输入对照，局部λ_x从对照估计；ATAC-seq常无输入对照，MACS直接用ATAC-seq读数覆盖估计λ_x。但ATAC-seq读数可因可测序DNA量局部变化而偏倚，如拷贝数变异（copy number variations, CNVs）。此外，DNA复制时序（DNA replication timing, RT）也影响DNA剂量：DNA按精确时空程序复制，具局部、细胞类型特异和进化保守性，因此分裂细胞中早期复制区域比晚期复制区域更丰富、更易被测序。直接从ATAC-seq信号计算λ可能掩盖RT引起的DNA剂量局部变异。因此，研究人员旨在评估将WGS（whole-genome sequencing，全基因组测序）衍生信息作为MACS流程对照以计算λ的影响，类似ChIP-seq中输入对照。研究人员观察到，整合WGS识别ATAC峰增加峰数和宽度，并推断更多靠近转录起始位点（transcription start sites, TSSs）的峰。这些峰具有经典CRE表观遗传标记，富集与特定性状和疾病相关遗传变异。研究人员还显示，使用WGS识别的ATAC峰数量与基因表达相关性优于无对照峰。结果表明，将WGS数据纳入MACS流程通过更好考虑局部DNA剂量效应，提升ATAC-seq峰识别准确性。该方法作为开源Nextflow流程发布，挑战所有位点具有相同可见性的假设，对广泛DNA测序方法有重要意义。该研究由Christophe Vroland、Mohammad Salma、Artemy Zhigulev、Pelin Sahlen、Eric Soler、Laurent Bréhélin和Charles-Henri Lecellier等完成。

关键技术方法概括（不超过250字）
研究人员使用ENCODE K562（ATAC-seq ENCSR483RKN；WGS ENCSR711UNY）和MCF-7（ATAC-seq ENCSR422SUG；WGS SRR8652105/PRJNA523380）数据。CNVs用Control-FREEC计算；Repli-seq数据来自GSE148362。WGS经fastp质控、BWA-MEM比对至hg38、Samtools和Picard处理，生成局部读深BigWig。峰识别比较MACS3标准流程与整合WGS对照流程，采用伪重复、合并及--call-summits。利用支持向量机回归（support vector machine regressor, SVR）预测ATAC信号，Fisher精确检验、bedtools、RNA-seq TPM、H3K27ac ChIP-seq、CAGE、ReMap、STARR-seq和GWAS/ClinVar/gnomAD/phyloP等分析。

研究结果
CNV and replication timing are confounding factors of ATAC-seq
研究人员从ENCODE获取K562 ATAC-seq信号，以K562作为tier-1细胞系。通过WGS调用CNVs，利用已发表七组分Repli-seq数据确定RT，其中S1–S6对应S期渐进阶段，第七组分G2代表最晚复制区域。示例位点显示，ATAC-seq信号在拷贝数扩增区（超过两个DNA拷贝）和早期复制S1区更高，在晚期复制区降低。全基因组水平，ATAC-seq信号与拷贝数状态和RT相关。MACS无外部对照计算的λ参数与CNV和RT高度相关。虽然CNV此前已被认为是ATAC-seq峰识别混杂因素，但RT是额外混杂因素。研究人员建立SVR模型，以CNV和RT为预测变量，能以高准确度预测ATAC-seq信号，Spearman ρ约0.6，RT相关特征最重要。由于CNV和RT可估计且与WGS覆盖相关，研究人员评估在ATAC-seq峰识别中用WGS控制DNA剂量变异的益处。

Controlling for DNA dosage variations in ATAC-seq peak calling
研究人员修改MACS峰识别流程，从匹配WGS数据估计λ，类似ChIP-seq输入对照。标准流程与ENCODE略有差异，使用不同split seed生成伪重复并使用MACS3而非MACS2，但结果峰仅轻微差异。用WGS考虑局部DNA剂量变异改善MACS λ参数与CNV及无CNV区域RT的相关性。两种流程峰Poisson P值分布相当。ATAC-seq折叠富集（fold-enrichment, FE）以标准λ计算时，非CNV区域可比扩增区强；WGS提供更大对比，扩增区信号强于早期复制区。因此，WGS计算的λ更好捕获CNV和更精细RT。WGS λ在扩增区升高，在早期复制区降低，导致扩增区检测峰比例减少、早期复制区峰增加。两种流程峰summit在CNV边界附近分布相似。

Comparing ATAC-seq peaks called with and without WGS
WGS识别峰总数和summit数超过标准流程，但多数峰由两种流程共同发现。WGS识别150 050个峰区域、448 359个summit；标准识别117 397个峰区域、230 870个summit。WGS窄峰几乎为宽度两倍：平均/中位长度624/496 bp，标准376/334 bp。定义new峰（仅WGS对照发现，n=58 279）和common峰（两流程均发现，n=91 771）。WGS整合流程恢复大多数标准峰，说明多数显著峰保留。标准峰富集于注释TSS，new峰也位于其紧邻区域。整合WGS可检测正好位于基因TSS及附近的峰；标准只识别第一类，遗漏许多TSS附近峰。new峰ATAC信号略低于common峰，T_n5序列偏好相似。结论：标准MACS因缺匹配对照掩盖许多开放染色质区域。

Using WGS as control improves ATAC-seq peak calling
new峰显著富集Alu元件，其有助于增强子-启动子特异性。new峰含表达数量性状位点（expression quantitative trait loci, eQTLs）比例略高于common峰。在K562正交数据中，H3K27ac（组蛋白H3第27位赖氨酸乙酰化，活跃增强子标志）信号与new峰比common峰更一致；全基因组H3K27ac峰summit更常落在new峰内。ReMap TF ChIP-seq和CAGE（cap analysis of gene expression）峰summit也如此。STARR-seq（self-transcribing active regulatory region sequencing）测量增强子活性，new峰信号略低于common但高于非峰区。new峰富集GWAS（genome-wide association study）表型性状和ClinVar疾病相关致病变异。使用gnomAD和phyloP，new与common峰间个体变异和种间保守性相当。更宽WGS峰（>2 kb）在所有测试基因组注释中富集，长度归一化至1 kb后仍如此，证实生物学相关性。结论：WGS对照改善ATAC-seq峰识别。

Integrating WGS as control reveals cell-specific ATAC-seq peaks
比较K562和MCF-7峰识别，new峰比common峰更细胞特异：K562 new峰仅41.6%在MCF-7检测到，common峰为63.4%。基因周围开放染色质区域数量与表达水平强相关。研究人员计算RNA-seq基因表达与WGS整合流程识别的基因起始位点2 kb窗口内ATAC峰数量的相关性，类似但略高于标准流程。进一步测试WGS峰是否更好预测差异基因表达：计算K562与MCF-7间基因起始位点附近WGS ATAC峰比值，与RNA-seq计数倍数变化（K562/MCF-7）相关，WGS峰相关性超过标准峰。K562特异峰关联基因在K562中表达略高于MCF-7，MCF-7特异峰反之。K562特异峰比MCF-7特异峰更靠近最近基因起始。GO富集：以K562特异峰基因为前景、MCF-7为背景，富集免疫反应调节；反向富集细胞运动和迁移调节。结论：WGS新识别峰对基因表达和细胞特异性有信息。

Low-depth WGS is sufficient to enhance ATAC-seq peak calling
为确定所需测序深度，研究人员对WGS数据下采样，发现深度5足以恢复先前分析所用峰。该深度表明将WGS对照样本纳入ATAC-seq分析成本有效，值得常规实施。

讨论与结论
讨论部分总结：研究人员认为，ATAC-seq信号不仅与CNV强相关，也与RT强相关。既往研究记录ATAC-seq与RT关联，早期复制区可及性更高、晚期复制区较低，但其潜在混杂效应未被考虑。为阐明DNA剂量变异对ATAC-seq数据分析影响，研究人员检查其对MACS峰识别流程背景估计的影响，并整合WGS作为对照，类似常规ChIP-seq输入样本。整合WGS对照后，结果与无匹配输入数据的标准MACS明显偏离：峰更多、更宽、更常近TSS，并改善差异基因表达预测，可能影响基于增强子的基因调控网络。因此，考虑DNA剂量可揭示原本未检测到的调控区域。正交数据支持新峰的生物学相关性。研究人员主张改进ATAC-seq峰识别策略，明确整合DNA剂量信息。但不排除其他或复合对照模型同样或更有效，例如将GC含量和Repli-seq RT谱与WGS整合可能进一步改善背景归一化。RT显著影响峰识别，尤其在高增殖细胞系如K562和MCF-7中。RT是时间和细胞特异特征，比较异质或癌症来源细胞类型染色质可及性时需考虑复制动态。DNA剂量与染色质可及性关系也支持从ATAC-seq数据推断CNV。尽管bulk和单细胞ATAC-seq在输入需求、数据稀疏性和分析范围差异大，多数单细胞流程仍依赖聚合读数谱和基于MACS的共识峰识别，因此本文所述混杂效应可能延伸至单细胞ATAC-seq。类似从scATAC-seq推断CNV，RT也可从该类数据恢复。若低覆盖WGS足够，未来方法可整合匹配单细胞WGS数据集，提供本文bulk WGS对照的直接对应。ATAC-seq结果提示RT和CNV可能是广泛基于DNA测序技术的普遍混杂因素。RADICL-seq相互作用计数也与CNV和RT高度相关，即便考虑RNA特异校正后仍可能显著混杂。Hi-C归一化如迭代校正通常假设位点可见性相等，但CNV时扩增区更易被捕获，低拷贝区可能代表性不足；RT也可能引入偏倚。鉴于Hi-C接触频率与RT强相关，可能需要改进Hi-C归一化以更好考虑复制相关偏倚。解决这些挑战将提高数据准确性并深化对DNA复制、染色质组织和拓扑关联域形成机制联系的理解。

研究结论翻译：总体而言，这些结果表明，整合WGS作为对照可显著提高ATAC-seq峰识别的准确性。关键的是，即使低深度WGS数据也足以改善峰识别，使该方法具有成本效益和可扩展性。通过挑战均匀基因组可见性假设，该方法对其他基于DNA测序的技术也具有广泛意义。