拿到ATAC-seq数据的那一刻我才意识到自己之前把这事儿想简单了。当时我手里的FASTQ测序报告很漂亮Q30超过90%看起来是个标准的好数据。结果到了比对这一步alignment rate只有55%我一度以为是参考基因组下载错了。后来一层层排查才发现是建库用的Nextera接头没有切干净read2末端全是接头序列bowtie2硬生生把接头当基因组匹配把大量reads丢弃掉了。这种数据表象很好、往下跑全是坑的经历我相信很多跑过ATAC-seq的人都有同感。这篇文章是我做ATAC-seq数据分析全流程的学习记录和系统总结覆盖从原始FASTQ到最终差异peak、motif分析的完整链路。不只给命令还会解释每个参数为什么这么设、每个QC指标到底在衡量什么、以及我在实际项目中踩过的坑和排查思路。适合刚开始接触ATAC-seq数据分析的生信同学也适合做实验想自己看懂分析结果的朋友。1. 先搞明白ATAC-seq到底在测什么Tn5机制决定了下游分析的每个细节1.1 Tn5转座酶不是切一刀那么简单ATAC-seq的核心试剂是高活性Tn5转座酶。它同时干了两件事一是识别并切割开放染色质区域的DNA二是把测序接头直接连到切口两端。这就是为什么ATAC-seq建库速度极快、起始量极低——打断和加接头在一步内完成没有传统建库里独立的片段化末端修复A尾连接接头流程。这个机制带来一个直接后果Tn5的切割偏好决定了你看到的信号分布。它优先进入核小体间开放的区域核小体占据的区域基本不被切割。所以reads聚集的地方就代表染色质开放、有调控潜力的区域。这也是为什么ATAC-seq信号能同时反映转录因子结合位点、增强子、启动子等调控元件的开放状态。很多人把ATAC-seq当成ChIP-seq来跑这个类比帮了倒忙。ChIP-seq是抗体富集后再测序信号的分布由抗体特异性决定ATAC-seq没有富集步骤信号就是全基因组范围内Tn5切割事件的投影背景和噪声的逻辑完全不同。比如ChIP-seq里不太需要担心的线粒体reads在ATAC-seq里可以占到全部有效reads的20%-50%因为线粒体基因组小而且高度开放Tn5特别爱切它。1.2 Tn5的9 bp交错切口与read偏移校正这是ATAC-seq数据分析里最容易被忽略、却直接影响peak质量和motif定位准确性的细节。Tn5转座酶在DNA上切割时两条链的切口位置不是平齐的而是错开9 bp。测序得到的read实际上是从这个交错切口的两个端点分别延伸出来的。反映到比对结果上就是正链reads的起始位置和真实Tn5切割位点之间差了4 bp负链reads的起始位置差了5 bp。如果不做校正直接堆叠信号你会看到每个真实Tn5结合位点两侧各出现一个小峰形成双峰形貌MACS2在这种情况下call出来的peak中心位置会系统性偏移后续做motif分析时motif的中心往往对不准peak summit富集显著性也会被稀释。实际操作中通常把正链reads向右移动4 bp、负链reads向左移动5 bp让两条链的信号对齐到同一个cut site。这一步在deepTools里有现成参数alignmentSieve --ATACshift后面我会给出具体命令。做过这个校正和没做过在TSS区域的信号图上差异非常明显——校正后的TSS峰更尖、更深。1.3 插入片段长度分布是免费的质检报告ATAC-seq的文库由不同长度的DNA片段组成无核小体区域产生的片段一般在50 bp-150 bp单核小体片段在180 bp-250 bp双核小体大约350 bp-450 bp以此类推。把比对上的双端reads按照插入片段长度画直方图正常情况下能看到明显的约200 bp周期性——这就是经典的核小体阶梯。这个特征本身就是重要的质控指标如果建库成功、Tn5切割正常阶梯必然存在如果片段长度分布杂乱无章说明建库失败或者样本已经降解。同时片段长度分布还告诉你ATAC-seq的双端reads插入片段跨度很大。比对时如果参数限制得太死比如把最大片段设成500 bp虽然能保留下大部分无核小体信号但多核小体片段会被丢掉设成2000 bp是比较稳妥的上限这也是为什么后面对比参数里会专门提到-X 2000。2. 从FASTQ到BAM我把每一步都在干什么拆开讲透2.1 原始数据质控与Nextera接头切除拿到FASTQ的第一件事永远是FastQC这没什么好说的。跑完后用MultiQC把多个样本汇总到一起看省时间。这里要强调的不是怎么跑而是看什么除了常规的Q20/Q30、GC含量、碱基分布一定要留意Overrepresented sequences那一栏。如果里面出现大段的接头序列或者重复序列说明接头没有切干净或者文库复杂度有问题。ATAC-seq建库用的是Tn5转座酶配套的Nextera体系接头序列里有一段特征性的19 bp序列CTGTCTCTTATACACATCT。这是Tn5的mosaic end序列在read的3端一旦出现就说明测序长度已经穿透插入片段、进入到了接头区域。这部分序列不切掉后续比对时会被当成基因组序列去匹配轻则浪费比对深度重则造成大量错配和比对失败。我用cutadapt切接头双端reads同时处理cutadapt \ -a CTGTCTCTTATACACATCT -A CTGTCTCTTATACACATCT \ -m 30 -j 8 \ -o sample_trim_R1.fq.gz \ -p sample_trim_R2.fq.gz \ sample_raw_R1.fq.gz sample_raw_R2.fq.gz解释一下参数-a对应read1的3端接头-A对应read2的3端接头Nextera体系两条reads的接头序列一样-m 30表示切除后长度不足30 bp的reads直接丢弃这一步非常关键因为ATAC-seq有大量短片段接完接头后如果两个接头在中间碰头会产生全长几乎都是接头的reads这种reads没有任何分析价值-j 8是并行线程数。切完接头后再跑一次FastQC确认接头残留消失。这里有个经验如果没有切干净就直接比对alignment rate可能从95%掉到60%以下而很多人会误以为是样本质量问题去重新提RNA或者重新建库其实纯粹是数据处理环节的疏忽。2.2 比对参数的选择为什么用bowtie2而不是bwaATAC-seq比对最常见的选择是bowtie2主要原因是它对双端reads的局部比对控制更精细而且原版Buenrostro流程用的就是bowtie2实测表现稳定。bwa-mem也能用但在处理ATAC-seq这种大量短插入片段时bowtie2的--very-sensitive模式通常能得到更高的比对率。我的完整比对命令bowtie2 -p 16 --very-sensitive -X 2000 --no-mixed --no-discordant \ -x /path/to/GRCh38/genome_index \ -1 sample_trim_R1.fq.gz \ -2 sample_trim_R2.fq.gz \ -S sample.sam 2 sample_bowtie.log samtools view -bS - 16 sample.sam -o sample.raw.bam samtools sort - 16 -O BAM -o sample.sorted.bam sample.raw.bam samtools index sample.sorted.bam--very-sensitive是bowtie2的预设组合参数等效于开启更严格的seed匹配和更细的延伸策略慢一些但值得。-X 2000定义了双端reads允许的最大插入片段长度前面说过设成2000是为了保留多核小体级别的长片段。--no-mixed和--no-discordant的意思是丢弃那些只有一端能比上、或者两端比上但方向/距离明显不合理的reads。对ATAC-seq来说基因组上真正的Tn5插入应该都是正常的双端比对这两个参数能帮你去掉大量噪声。比对完成后先看sample_bowtie.log里的整体比对率。ATAC-seq的正常区间一般在85%-95%具体取决于样本质量和参考基因组的匹配程度。如果低于80%先别急着往下走后面第5章我会专门讲排查思路。2.3 去除线粒体reads、重复reads和Tn5偏移校正比对得到的BAM接下来要过三关顺序不能乱。第一关是去线粒体reads。samtools idxstats可以快速看到每条染色体的reads数统计。如果chrM的占比超过30%不用慌这在ATAC-seq里很正常线粒体基因组小、开放程度高Tn5就是爱切。问题在于这些reads会干扰后续的peak calling尤其会在线粒体区域堆积出极高的假信号所以必须去掉。samtools view - 8 -h sample.sorted.bam | \ awk !/^/ $3 ! chrM || /^/ | \ samtools view -b - 8 -o sample.nuclear.bam第二关是去PCR重复。建库过程中PCR扩增会让同一个片段被复制成多个完全相同的reads这在统计上会被当成多次独立的Tn5切割事件人为抬高信号。用picard处理picard MarkDuplicates \ Isample.nuclear.bam \ Osample.dedup.bam \ Msample.dedup.metrics.txt \ REMOVE_DUPLICATEStrue这里有一个ATAC-seq特有的争论Tn5切割是随机的不同细胞来源的片段如果恰好切割位置完全相同在数据上看起来也是重复readsMarkDuplicates无法区分PCR重复和天然相同都会去掉。对于起始量高的细胞系样本这个损失可以接受对于起始量极低的单细胞或者稀有样本去重可能导致真实信号大量丢失。我的做法是常规样本去重超低起始量的样本保留重复但记录一下dup比例在MACS2里用--keep-dup all并在结论里说明。关键是同一个项目里所有样本的处理方式必须一致否则比较分析会引入系统性偏差。第三关是Tn5偏移校正。这一步把正链read起点向右移4 bp、负链read起点向左移5 bpalignmentSieve --ATACshift --numberOfProcessors 16 \ -b sample.dedup.bam \ -o sample.shifted.bam samtools index sample.shifted.bamalignmentSieve是deepTools里的工具--ATACshift参数专门为ATAC-seq设计。它会在BAM层面做偏移校正同时保留end motif信息。校正完的这个BAM既是后续peak calling的输入也是生成bigWig可视化的输入。最后生成bigWig用于IGV可视化注意这一步要用校正后的BAMbamCoverage -b sample.shifted.bam -o sample.bw \ --normalizeUsing CPM --binSize 10 -p 16用CPM归一化是为了让不同样本间的信号量纲一致--binSize 10表示每10 bp一个信号bin兼顾分辨率和文件大小。多组学比较时这个归一化方式要提前统一不然画出来的热图差异可能只是测序深度的差异。3. Peak calling与质量评估这些指标比peak数量更重要3.1 MACS2的两种跑法理解比抄命令重要MACS2是ATAC-seq peak calling的事实标准但参数用法因数据而异我见过太多人因为抄错参数跑出一堆奇怪的peak。核心问题在于ATAC-seq的reads是Tn5切割后产生的短片段没有ChIP-seq那种抗体富集带来的固定片段长度所以MACS2默认的建模逻辑不完全适用。有两种主流跑法。第一种是Shifted BAM配BAMPE模式macs2 callpeak \ -t sample.shifted.bam \ -f BAMPE \ -g hs \ -n sample \ -q 0.05 \ --keep-dup all \ --call-summits-f BAMPE让MACS2直接利用双端reads定义的插入片段信息每个片段的两个5端就是Tn5的两个切割点这也正是shifted BAM的意义所在——校正后片段的5端精确对应Tn5结合位点。这种跑法保留了真实片段长度信息适合对peak位置精度要求高的场景。第二种是老流程里更常见的ChIP模拟跑法macs2 callpeak \ -t sample.dedup.bam \ -f BAM \ -g hs \ -n sample \ -q 0.05 \ --nomodel \ --shift -100 \ --extsize 200 \ --keep-dup all \ --call-summits--shift -100 --extsize 200的含义是把每条read当作一个200 bp的测序片段来处理先用--shift -100把reads的5端向5方向移动100 bp再向3方向覆盖200 bp。这样单端reads也能近似出双端片段的覆盖效果。这个魔法数字的历史包袱很重最初是为了在单端测序时代模拟核小体占位现在双端测序已经很普及我个人更推荐第一种BAMPE跑法它不会强行把所有片段压成200 bp。两个参数必须注意--keep-dup要和前面的去重策略一致。如果你已经用picard去了重就填all让MACS2不再做任何重复过滤如果你没去重填auto让MACS2自己评估重复率但两种策略不要混用。-g hs是指定有效基因组大小人基因组常用2.7e9小鼠用1.87e9这个值影响背景模型的估计设错了peak数量会明显异常。跑完后同一套参数处理所有样本得到narrowPeak文件。3.2 blacklist过滤这是最容易被跳过的步骤即便MACS2结果正常也不能直接拿去注释。基因组里有大量区域——着丝粒、端粒、rDNA重复区、卫星重复等——在没有任何真实染色质开放信号的情况下也会因为序列比对歧义和Tn5偏好堆积出高信号。ENCODE维护了一套blacklist文件按基因组版本区分覆盖这些假信号区域。bedtools intersect -v \ -a sample_peaks.narrowPeak \ -b blacklist_hg38.bed sample_peaks_clean.narrowPeak过滤后在blacklist里检查一下被去掉的peak比例。正常情况下应该低于2%-3%。如果超过5%说明数据的背景噪声偏高或者比对质量存在问题不要急着掩盖它回到比对和去重环节排查。3.3 QC指标TSS富集、FRiP、PBC、NSC/RSC怎么看peak call完了不等于分析做完了。在进入下游之前我会把所有样本的QC指标系统地过一遍。下面这几个指标是我最看重的附上通常的参考阈值具体项目里以ENCODE指南和同批次数据横向比较为准。指标衡量什么参考阈值常用工具TSS富集启动子区开放信号的集中程度 8好的样本往往10-20deepTools / ATACseqQCFRiPFraction of Reads in Peaks 0.2bedtools 手算PBC文库复杂度PCR重复水平 0.9picard / preseqNSC / RSC信噪比与交叉相关强度NSC 1.05RSC 1phantompeakqualtools插入片段阶梯核小体周期性mono-nuc峰清晰bamPEFragmentSizeTSS富集是最直观的指标之一。用deepTools计算并画图computeMatrix reference-point \ -S sample.bw \ -R hg38_tss.bed \ -b 2000 -a 2000 \ --referencePoint TSS \ -p 16 -o sample_TSS_matrix.gz plotProfile -m sample_TSS_matrix.gz \ -o sample_TSS_profile.pdf \ --plotTitle TSS enrichment正常的ATAC-seq数据TSS中心位置会有一个尖锐的峰两侧信号低。峰越高越尖说明开放信号越聚焦在启动子区域如果TSS富集低于5或者TSS区域根本看不出峰那么后续所有结论都要打问号。顺便说一句hg38_tss.bed在UCSC和ENCODE都能下载注意版本和GTF注释文件要配套。FRiP的计算逻辑很简单就是比对上的reads里有多大比例落在peak区域。我的统计方式total$(samtools view -c sample.dedup.bam) in_peaks$(bedtools intersect -u -a sample.dedup.bam -b sample_peaks_clean.narrowPeak | wc -l) awk -v t$total -v p$in_peaks BEGIN{print FRiP , p/t}FRiP在ATAC-seq里通常比ChIP-seq高因为可及性信号是弥散分布在全基因组调控区域的不是只有少数几个富集位点。低于0.2时我会先怀疑是不是peak calling参数过严再怀疑数据质量问题。NSC和RSC来自phantompeakqualtools原理是基于read交叉相关的信噪比分析。跑一遍能拿到一串指标数值含义是NSC大于1.05说明信号对比背景可区分RSC大于1说明存在明显的核小体或转录因子尺度的富集峰。这两个指标在比对后、cutadapt后、去重后各阶段都可以检测用来判断哪一步操作改善了数据质量。3.4 peak注释位置分布不是终点得到干净的narrowPeak后注释是常规操作。R里的ChIPseeker是我用得最多的一步就搞定library(ChIPseeker) library(TxDb.Hsapiens.UCSC.hg38.knownGene) peak - readPeakFile(sample_peaks_clean.narrowPeak) annot - annotatePeak(peak, TxDb TxDb.Hsapiens.UCSC.hg38.knownGene) plotAnnoPie(annot)注释结果的解读要小心一个现象如果样本里大部分peak都落在promoter区域确实符合某种细胞状态但也可能是因为建库时细胞固定过度导致只有启动子区域开放。反过来如果大量peak落在异染色质/基因间区需要检查blacklist是否过滤干净。把peak位置分布和已知生物学背景放在一起看而不是机械地数注释比例。4. 多组比较与motif分析从差异peak到生物学结论4.1 构建一致peak集与差异可及性分析如果你的实验有多个条件想知道不同条件下染色质可及性哪里变了第一步不是直接比较每个样本的peak而是构建一个统一的一致peak集。原因很简单不同样本的peak边界会有微小差异直接比较各自的peak列表会被边界噪声主导。我常用的策略是同一条件的重复样本先取重叠peak要求重叠比例超过50%再把不同条件合并成一个总的一致peak集。复杂的方案是用IDRIrreproducible Discovery Rate筛选两个重复之间的高置信peak但IDR对重复数和质量问题很敏感日常项目里用bedtools intersect取交集就够。得到一致peak集后统计每个样本在这个peak集上的reads计数。R里用DiffBind把样本表-计数-归一化-差异分析串起来非常省事library(DiffBind) dba - dba(sampleSheet samples.csv) dba - dba.count(dba, bUseSummarizeOverlaps TRUE) dba - dba.normalize(dba) dba - dba.contrast(dba, categories DBA_CONDITION) dba - dba.analyze(dba, method DBA_DESEQ2) res - dba.report(dba, th 0.05, bUsePval FALSE)这里有两个经常出错的地方。一是样本的测序深度差异如果处理组平均深度是对照组的两倍不做归一化直接做差异分析会得到大量假阳性peakDiffBind内置的归一化能缓解问题但设计实验时最好还是让各组的测序深度尽量接近。二是生物学重复的数量少于3个重复时DESeq2的离散度估计很不稳定结果只能当筛选不能当定论。如果你不想用DiffBind也可以用featureCounts从一致peak集统计原始count矩阵再用DESeq2走标准流程结果大同小异。4.2 motif分析的背景控制比命令更重要的坑差异peak拿到手大多数人下一步就是跑motif富集。这里最大的坑不是工具不会用而是背景选择错误导致结果全是无效motif。HOMER是最常用的motif富集工具基本命令findMotifsGenome.pl \ sample_up_peaks.bed hg38 \ homer_up_out \ -size 200 -len 6,8,10,12 -p 8-size 200决定取peak中心周围多大的序列范围这个建议匹配peak的宽度量级-len指定motif长度范围6-12 bp覆盖了绝大多数转录因子结合位点的宽度。但假如你的差异peak大量集中在启动子区域而HOMER默认使用的是全基因组随机背景你富集出来的很可能是Sp1、NFY这类泛启动子motif它们只是启动子区域的普遍特征和你的实验条件没有直接关系。更聪明的做法是直接拿下调peak当上调peak的背景做两组对照findMotifsGenome.pl \ sample_up_peaks.bed hg38 \ homer_up_vs_down \ -bg sample_down_peaks.bed \ -size 200 -len 6,8,10,12 -p 8这样富集出来的motif是上调可及性区域vs下调可及性区域之间真正有区分度的结合位点而不是启动子通用特征。我在实际项目里用这个策略后结论贴合生物学预期很多。另一个常见问题是peak数量太少比如只有几百个时motif富集统计功效不足很容易出假阳性。可以先看看差异peak的数量是否足够通常要至少1000个以上再做motif富集才有意义。4.3 footprint分析深度不够不要强行做转录因子足迹分析footprinting是基于ATAC-seq数据的进阶玩法。原理是转录因子结合在DNA上时Tn5切不到被蛋白质占据的区域所以在转录因子结合位点内部会出现信号低谷两侧信号较高形成脚印。这个分析能给出转录因子结合的证据但它对数据质量的要求非常苛刻。如果样本的可用reads数低于5000万人基因组footprint基本不可信。即便深度足够Tn5自身的切割偏好也会制造大量假脚印需要用TOBIAS里的ATACorrect步骤先校正Tn5偏好TOBIAS ATACorrect \ --bam sample.dedup.bam \ --genome hg38.fa \ --peaks sample_peaks_clean.bed \ --outdir atacorrect_out TOBIAS FootprintScores \ --signal atacorrect_out/sample_corrected.bw \ --regions sample_peaks_clean.bed \ --output sample_footprints.bw我在低深度样本上试过一次footprint结果得到的足迹在随机区域也大量出现后来加了深度门槛只对测序量达标的样本做footprint结论才稳定下来。拿不到高质量数据时不要硬凹footprint图。4.4 与基因表达整合分析的思路ATAC-seq的最终目的是解释调控机制所以把可及性变化和转录组变化整合起来是常见需求。简单做法是用GREAT对差异peak做功能富集看它们可能影响哪些通路进阶做法是给peak匹配最近基因或共定位的基因再做表达相关性分析。但这里要克制一点远距离调控关系增强子-基因在ATAC-seq里只能靠相关性推断不能证明因果。即使某个增强子在某个基因附近且两者信号正相关也不代表这个增强子真的调控这个基因。如果要进一步收敛可以考虑叠加Hi-C数据或者做CRISPR验证这就超出纯ATAC-seq数据分析的范畴了。把可及性变化当作候选调控区域列表来使用是更稳妥的定位。5. 我实际踩过的坑比对率低、peak数量异常、重复性差5.1 比对率只有55%一次完整的排查链路开头提到的那批数据alignment rate只有55%我的排查过程很有代表性。第一步回看bowtie2日志。系统记录显示总reads 4000万比对上的只有2200万未比对reads中大量是mate mate比对失败或者discordant。这说明不是参考基因组索引问题而是reads自身有问题。第二步重看cutadapt报告。我发现跑trim时作业脚本里路径写错了程序直接跳过了cutadapt这一步reads根本没有切除接头。第三步手动截取几万条未比对reads用FastQC看序列组成结果Overrepresented sequences栏赫然出现大段Nextera接头序列。到这里基本坐实接头污染导致read2在末尾的几十个碱基全是接头bowtie2在局部比对时无法跨越这段非基因组序列大量reads因此判为未比对。第四步重新跑了cutadapt把接头切干净后再比对alignment rate回到94%。这个过程让我养成了一个习惯比对前一定核对trim日志比对后顺手用FastQC抽查几个样本确认接头已清除。自动化的流程如果缺乏中间确认环节一个小疏漏会影响后面所有步骤。5.2 peak数量太多或太少先定位再调参我曾经处理过一个细胞系样本MACS2跑出18万个peak明显不对。逐一排查下来原因有三个叠加blacklist没过滤贡献约2万个假peak、没有去除线粒体reads在chrM上堆积出大量假峰、PCR重复没有去峰被重复reads人为抬高。反过来另一个极端是peak只有几百个。那次是数据本身的可用reads太少全样本比对成功后去掉chrM、去重后剩下不到800万readsTSS富集只有2这种样本无论怎么调MACS2参数都救不回来。遇到这种情况我的建议是老实回到实验环节提高细胞活力、增加有效reads测序量而不是在分析上折腾。用一张表总结我现在的排查逻辑现象优先级排查项验证方式比对率80%接头切除、参考基因组版本、样本物种FastQC、bowtie日志peak数量异常高blacklist、chrM、重复未去、q阈值过松分层检查各环节peak数变化peak数量异常低有效reads不足、MAPQ过滤过狠、MACS2参数idxstats、TSS富集重复样本相关性低批次效应、测序深度差异、样本来源差异PCA、pairwise相关热图5.3 重复性与批次效应别到最后一步才想起来多组学项目的失败往往不在单个样本的质量而在样本间的一致性。两个生物学重复之间信号相关性低于0.85就该警惕了。我常用的检查方法是deepTools的multiBigwigSummary加plotCorrelationmultiBigwigSummary bins \ -b ctrl_rep1.bw ctrl_rep2.bw treat_rep1.bw treat_rep2.bw \ -o all_samples_corr.npz plotCorrelation \ -in all_samples_corr.npz \ --corMethod pearson \ --whatToPlot heatmap \ --plotNumbers -o all_samples_corr.pdf做完相关热图之后再对peak计数矩阵做一次PCA看样本是否按照生物条件聚在一起而不是按照测序批次聚在一起。如果PCA显示同条件下两个重复分别落在两端那就是批次效应后续差异分析结果的可靠性要大打折扣。处理批次效应的方式有很多RUVSeq、ComBat-seq等但最好的办法还是设计实验时把同一条件下重复样本分散到不同批次从源头规避。最后再分享一个我自己的习惯所有中间文件必须保留处理参数全部写进脚本而不是靠记忆每个样本的cutadapt日志、bowtie日志、MarkDuplicates的dup rate、MACS2的peak数量汇总成一个总表。这样当某个样本在最后一步出现异常时我能在一分钟内回溯到具体是哪个环节出了问题。ATAC-seq数据分析不是一个黑盒但如果你自己都不记录它就真的成了黑盒。