ATAC-seq数据分析全流程：从Tn5酶原理到peak calling实战

1. 先从整体上拆解ATAC-seq：我们到底在测什么

1.1 一个Tn5酶就是一整套“切割+标记”系统

ATAC-seq全称是Assay for Transposase-Accessible Chromatin with high-throughput sequencing，核心主角是Tn5转座酶。理解这个实验的关键，是先理解Tn5这个酶在干什么。它不像普通限制性内切酶那样只负责切断DNA，它是“切割”和“标记”一起完成的：识别开放染色质区域后，切断DNA的同时，把自己携带的测序接头直接连到断口两端。这一步做完，你的文库其实已经带着P5和P7接头了，后面只需要补几个PCR循环就能上机测序。

这也是我最初理解ATAC-seq时绕的一个弯：它跟ChIP-seq那种“先把DNA切碎，再在两端补接头”的逻辑不同，ATAC-seq在建库这一步就已经完成了片段化与接头连接的耦合。这个机制决定了后续数据分析里很多现象：比如插入片段长度分布天然呈现核小体周期性的振荡（约200bp一个周期），比如Tn5偏好性会导致某些区域覆盖度异常高，比如数据里会有大量线粒体reads——因为线粒体基因组没有核小体包裹，染色质状态是完全“开放”的，很容易被Tn5切碎。

熟悉了这套底层逻辑，你再去看数据就不会觉得有些现象是“异常的”，而是“本来就应该这样”。这也是为什么我建议不要只背流程，先把Tn5的工作机制搞明白，后面所有参数调整都围绕“这个酶到底切了哪些地方”来理解。

1.2 ATAC-seq能回答什么问题，以及它和ChIP-seq的区别

ATAC-seq解决的问题是“染色质开放程度”。染色质开放性对应的是调控元件的活跃状态：启动子、增强子、绝缘子等调控区域，一般会暴露出可供转录因子结合的开放染色质。检测这些区域的开放状态，能帮你回答这样几类问题：

• 不同细胞类型或处理条件下的开放染色质差异在哪里，哪个增强子被激活或沉默；
• 转录因子结合位点富集在哪些区域，推测某个TF在目标区域是否参与调控；
• 核小体定位与基因表达的关系，尤其是启动子区域核小体占位变化；
• 与RNA-seq联合分析，寻找开放区域与差异表达基因的关联。

很多人会把ATAC-seq和ChIP-seq放在一起问到底选哪个。我的看法是，它们的定位完全不同。ChIP-seq测的是“某个特定蛋白结合在哪些位置”，它依赖抗体质量，而且一次只能看一个转录因子或组蛋白修饰。ATAC-seq测的是“所有开放区域”，没有抗体依赖，一个样本就能拿到全基因组尺度的调控信息。代价是它看不到具体的转录因子结合，只能通过motif分析去推测哪些TF可能结合在开放的peak区域。

所以实际研究里，我更倾向于把ATAC-seq当作“先遣侦查兵”：先扫一遍全基因组的开放情况，缩小候选调控位点范围，再用ChIP-seq验证特定TF或组蛋白修饰是否真的结合在那里。两者互补而不是互斥。

1.3 分析流程总览与环境准备

一套完整的ATAC-seq数据分析流程，我把它分成五个阶段：质控与预处理、比对与过滤、peak calling、下游功能分析、可视化与结果解读。这里我直接给出一张我在实际项目里使用的流程总览，你可以照着搭：

关于硬件，我之前在实验室里用16核32线程、128GB内存的服务器跑过不少ATAC-seq数据。一个人类样本的FASTQ（大约5000万到8000万条reads），比对到参考基因组用Bowtie2大概耗时15到25分钟，MACS2 peak calling只需几分钟。如果你们实验室只有笔记本，也不是不能跑，但建议用预处理后的数据，或者直接下载公共数据（如ENCODE、GEO）来学流程。

环境方面，我个人强烈推荐用conda管理生物信息工具：

conda create -n atacseq -c bioconda -c conda-forge python=3.9 \ fastqc fastp fastp multiqc bowtie2 bwa samtools picard bedtools \ macs2 deeptools homer

实测下来，conda安装工具能省去很多编译依赖的麻烦。不过也要注意，conda的版本有时候不是最新的，比如MACS2在conda里大概率是2.2.x，够用就好，不必追求大版本更新。

2. 上游分析实操：从原始数据到可用的比对结果

2.1 数据质控：不能只看Q30，还要关注这几个关键指标

很多人拿到FASTQ的第一件事就是跑FastQC，然后看到Per base sequence quality是绿色就放心了。但实际上ATAC-seq数据有几个比Q30更值得关注的指标，我建议你在质控阶段就养成记录这些数值的习惯：

reads总条数和有效比对率。ATAC-seq文库的复杂度跟细胞起始量、Tn5酶用量、PCR循环数都有关。建库环节如果细胞量太少，或者Tn5浓度偏高，可能出现大量重复reads，有效信息量骤降。我经手过一批数据，一个样本看起来有6000万条reads，去掉重复后只剩下2800万，有效利用率不到50%，这直接影响后续peak calling的深度。所以拿到fastq的第一件事，我就建议同时记录原始reads数和经过比对、去重后的有效reads数，前后一对就知道文库质量如何。

线粒体reads占比。因为线粒体基因组是裸露的，Tn5会优先切割，导致线粒体reads在总reads中占比很高。人类细胞系的ATAC-seq数据，线粒体reads占比经常在20%到50%之间，高的甚至超过70%。这个数字本身不代表建库失败，但如果占比过高（比如超过80%），会影响核基因组区域的测序深度，这时候就得考虑在建库环节优化细胞核提取，或者在分析环节把线粒体reads直接过滤掉。

插入片段长度分布。这也是ATAC-seq特有的质控指标。从BAM文件里提取插入片段长度画出来，你会看到约200bp周期性振荡的模式：第一个峰在0到100bp附近（无核小体区域），后面的峰以约200bp为间隔递减（单核小体、双核小体、三核小体）。如果这个振荡模式完全消失，说明染色质结构被破坏或者Tn5处理过度了。

数据处理这一步，我会先用MultiQC把FastQC报告汇总在一起：

fastqc -t 16 *.fastq.gz -o fastqc_raw/ multiqc fastqc_raw/ -o multiqc_raw/

然后对reads进行去接头和低质量过滤。ATAC-seq的reads有些是短插入片段，测序时容易出现R2引物直接读到R1的接头，所以切接头这步很重要。我用fastp比较多：

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ -h sample_fastp.html \ --detect_adapter_for_pe \ -q 20 -u 30 \ -l 35 \ -c

参数说明一下：-q 20表示碱基质量低于20的位点会被修剪，-u 30控制的是如果一条reads上有30%以上的碱基质量低，整条reads丢弃，-l 35是过滤后最短长度阈值。我没有把过滤条件调得特别狠，因为ATAC-seq的reads本身偏短，太激进会损失很多有效信息。

2.2 比对工具选择：Bowtie2还是BWA-MEM

ATAC-seq比对的主流选择是Bowtie2，理由有几个：它对短reads的比对速度快，消耗内存低，而且比对结果里能保留MAPQ信息用于后续过滤。BWA-MEM在处理长reads或需要split alignment的情况下有优势，但ATAC-seq的双端reads通常只有50到150bp，Bowtie2完全够用，跑得还快。

我常用的比对命令如下：

bowtie2-build hg38.fa hg38 bowtie2 -p 16 --very-sensitive -x hg38 \ -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ -S sample.sam 2> sample_bowtie2.log samtools view -bS -@ 8 -o sample.raw.bam sample.sam samtools sort -@ 8 -o sample.sorted.bam sample.raw.bam samtools index sample.sorted.bam

--very-sensitive参数会让Bowtie2花费更多时间换取更高的比对灵敏度，这个参数对ATAC-seq数据我建议开着，因为开放染色质区域的reads通常较短，转座子偏好性还会带来不均匀的比对难度，灵敏度不够会导致部分真实信号丢失。

比对的参考基因组版本一定要和后续注释、peak注释保持一致。比如人类数据，我一般统一用GRCh38/hg38。如果你用的是UCSC的hg19或者Ensembl的GRCh37，后续用ChIPseeker注释的时候就要指定对应的TxDb，否则坐标不一致，peak注释就会错乱。这个坑我踩过一次，注释结果大范围偏移，最后重跑了一遍，白白浪费半天时间。

比对完成后，先别急着过滤，看一眼比对报告的指标。我用samtools flagstat：

samtools flagstat sample.sorted.bam > sample.flagstat.txt

重点关注比对率（mapped ratio）。一个合格的人类ATAC-seq样本，比对率通常在90%以上。如果低于70%，先别急着往下跑，去排查建库问题或者参考基因组是否选择正确。我遇到过一种情况是样本来自大鼠细胞，但比对参考基因组用成了人类，比对率只有可怜的20%左右，这种情况下后续分析毫无意义。

2.3 过滤顺序有讲究：线粒体、重复、黑名单

比对完成之后，BAM文件里的reads并不都能直接用于peak calling，需要按照一定顺序过滤。这个顺序我建议固定下来，因为它会影响你对每个步骤结果的理解：

第一步：过滤线粒体reads。很多教程把这一步放在去重之后，但我的经验是放在最前面更合理——线粒体reads占比太高时，提前过滤能显著减小后续处理的数据量，跑得更快。操作就是用samtools把MT染色体上的reads剔除：

samtools view -b -@ 8 -h sample.sorted.bam \ -o sample.nomt.bam \ -U sample.mt.bam \ --exclude-chr MT samtools index sample.nomt.bam samtools flagstat sample.nomt.bam

这里用了-U选项把线粒体reads单独存一个文件，方便随时统计线粒体占比。人类参考基因组的线粒体染色体名是"MT"，在小鼠里是同样是"MT"，但有些老版本参考基因组的命名可能是"M"，注意确认一下。

第二步：去除重复reads。Tn5酶有一个特性，它倾向于在相同的插入位点重复切割，形成大量PCR重复。另外建库时PCR扩增也会造成重复。我一般用Picard的MarkDuplicates：

picard MarkDuplicates \ I=sample.nomt.bam \ O=sample.nomt.dedup.bam \ M=sample.markdup.metrics.txt \ REMOVE_DUPLICATES=true \ VALIDATION_STRINGENCY=LENIENT samtools index sample.nomt.dedup.bam

REMOVE_DUPLICATES=true是直接删除重复reads，有些教程会建议设成false然后只标记不移除，给下游filter一个选择空间。我个人的习惯是直接移除，因为ATAC-seq的peak calling是基于覆盖度的，重复reads会严重扭曲开放区域的覆盖度信息，保留它们相当于给高覆盖区域不断加权重。

不过有一点要注意，不要用samtools rmdup，这个工具是旧时代的产物，不能正确处理双端测序数据，会导致大量信息丢失。

第三步：过滤低质量比对的reads。我习惯用MAPQ阈值来过滤，对于Bowtie2比对结果，MAPQ低于10的reads不保留，这些reads通常是比对到多位置或者比对质量存疑的：

samtools view -b -@ 8 -q 10 sample.nomt.dedup.bam \ -o sample.nomt.dedup.q10.bam samtools index sample.nomt.dedup.q10.bam

第四步：过滤黑名单区域。ENCODE项目提供了多种物种的blacklist区域（包括常见的高信号假区域、着丝粒、端粒等），这些区域的reads即使比对成功也属于噪音。从ENCODE官网下载对应参考基因组的bed文件，用bedtools过滤：

bedtools subtract -A \ -a sample.nomt.dedup.q10.bam \ -b hg38.blacklist.bed \ > sample.final.bam samtools index sample.final.bam

注意bedtools subtract输出的是sam格式，需要手动转成bam并排序索引。如果照搬这一条命令，建议后面再加一步：

samtools view -bS sample.final.bam | samtools sort -O BAM -o sample.final.sorted.bam

2.4 别忘了做核小体信号的质量评估

过滤完之后，先别急着跑MACS2，我强烈建议先做一步核小体信号评估。这一步能直观反映你的ATAC-seq文库质量，而且非常快。

最直接的方式是用deeptools或者samtools提取插入片段长度分布。我从BAM文件里提取Tn5插入位点，实际上每个Tn5切割事件会产生两条reads，它们的前端位置相差4bp。分析时一般把reads比对结果的5'端朝正向链+4bp位置，反向链-5bp位置作为Tn5的插入中心，然后生成插入位点信号：

samtools view sample.final.sorted.bam | \ awk -F'\t' '{ if ($9 > 0) { start = $2 + 4; end = $2 + $9 - 5; } else { start = $2 + $9 + 4; end = $2 - 5; } if (end >= start) { print $1"\t"start"\t"end; } }' > sample_tn5.bed bedtools sort -i sample_tn5.bed | \ bedtools merge -d 100 -c 1 -o count > sample_tn5_merged.bed

然后你去IGV里看看合并后的峰信号，如果TSS附近有明显的信号富集，说明文库质量没问题。更定量一点的做法是直接计算Fragments per Thousand Transcripts（FTT）或者TSS富集分数，从ENCODE官网下载TSSbed文件，然后计算reads在TSS区域的富集倍数。正常ATAC-seq样本的TSS富集分数应该在5到10之间，低于4说明文库质量堪忧。

3. Peak calling与下游核心分析：把信号转成生物学结论

3.1 MACS2参数怎么调才靠谱

Peak calling这一步，我用得最多的还是MACS2。虽然也试过Genrich和SPRING，但MACS2胜在参数直观、结果稳定、社区案例多。ATAC-seq的MACS2调用和ChIP-seq有一个关键区别：ATAC-seq没有严格意义上的“control”样本，它通常用Tn5酶切割背景（比如等量基因组的裸露DNA）或者直接不做对照。我做细胞系样本时，一般就直接跑单样本模式。

我的常用命令如下：

macs2 callpeak -t sample.final.sorted.bam \ -n sample \ -f BAMPE \ -g hs \ -q 0.05 \ --shift -100 \ --extsize 200 \ --nomodel \ -B --SPMR \ --keep-dup auto

这里几个参数我分别说一下为什么这么设：

-f BAMPE是告诉MACS2输入的是双端比对结果，它会把每个reads对当作一个DNA片段来处理，而不是单纯把单端reads延伸成固定长度。这个参数对ATAC-seq很重要，因为片段长度信息本身包含核小体周期模式，不应该被忽略。

--shift -100 --extsize 200 --nomodel这三个参数配合使用，是ATAC-seq单样本模式下的常见做法。--nomodel跳过MACS2自带的模型构建步骤（因为没有control样本，模型构建容易失败），然后手动把reads的5'端向3'方向做定点延伸：每个Tn5插入位点前后各延伸100bp，实际上就是每个插入位点变成200bp的信号窗口。这个做法本质上还原了“以Tn5插入点为中心的高斯分布信号”，比默认行为更适合ATAC-seq数据。

--keep-dup auto表示在peak calling过程中，对于重复reads，会自动根据局部覆盖度决定保留多少。因为前面已经物理去掉了绝大多数重复reads，这里保留一些重复是为了避免过度剪切信号，这个参数可以保持默认。

-B --SPMR是输出bedGraph和每百万条reads的标准化信号，方便后续deeptools做可视化。

跑完之后，你的输出会有_peaks.narrowPeak、_peaks.broadPeak、_summits.bed几个核心文件。ATAC-seq默认用narrowPeak即可，因为开放染色质区域的peak通常比较尖锐。如果你的数据来自某些特殊细胞类型，信号很大很宽，也可以考虑用broadPeak看看，但绝大多数场景下narrowPeak已经足够。

关于-g参数，它是估算有效基因组大小（effective genome size）。人类的推荐值是hs（约2.7e9），小鼠是mm（约1.87e9）。不要用基因组实际长度去算，因为有很多重复区域和黑名单区域不可比。

3.2 Peak数量是多少才算正常

Peak calling跑完，很多人第一个问题是：我这个样本出了多少个peak才算正常？这里我给出一个经验范围，仅供参考：

人类细胞系的ATAC-seq样本，在中等深度（约3000万有效非重复reads）情况下，MACS2默认参数下通常能检测到5万到12万个narrowPeak。小鼠细胞系数量略少，大概在3万到8万之间。如果peak数量低于1万个，大概率是数据质量不过关或者细胞类型特别特殊；如果peak数量超过20万个，可能是质量不好导致信号太弥散，也可能是细胞状态高度开放。

当然这些数值不是铁律。我见过神经元样本因为整体染色质高度开放，peak数量接近20万；也见过某些沉默期的细胞系peak数量只有几千个。关键是看重复样本之间peak数量是否一致、peak是否富集在TSS和增强子标记附近，以及FRiP值（Fraction of Reads in Peaks，定位到peak区域的reads比例）。FRiP是一个非常重要的质控指标，计算方式很简单：

bedtools intersect -a sample.final.sorted.bam -b sample_peaks.narrowPeak -wa -u | wc -l total_reads=$(samtools view -c sample.final.sorted.bam)

用peak内reads数除以总数，得到FRiP。合格的人类ATAC-seq样本，FRiP普遍在0.3以上，即至少30%的有效reads落在peak区域内。低于0.2就要怀疑peak calling参数或者数据质量了。

3.3 Peak注释与差异分析：别只看“落在哪个基因上”

拿到peak列表后，很多人习惯性地去注释“peak落在哪个基因的启动子区域”，然后就结束分析了。这种思路太粗糙。Peak注释是一个需要结合基因结构、功能元件、转录方向来综合判断的过程，不能只看一个最近基因。

我推荐用ChIPseeker这个R包做注释：

library(ChIPseeker) library(TxDb.Hsapiens.UCSC.hg38.knownGene) library(org.Hs.eg.db) library(clusterProfiler) txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene peak_annot <- annotatePeak( "sample_peaks.narrowPeak", TxDb = txdb, annoDb = "org.Hs.eg.db", tssRegion = c(-1000, 1000), level = "gene" )

注释结果会告诉你每个peak属于启动子区（Promoter，通常分布在TSS上下游1kb到3kb）、5' UTR、3' UTR、外显子、内含子还是基因间区。我拿到注释后习惯先画一张饼图看区域分布，正常情况下一半以上的peak应该落在启动子和远端的基因间区（增强子候选），如果大量peak落在外显子或重复区域，说明数据或者peak calling可能有问题。

需要特别警惕的是“peak离最近基因很远不代表没意义”。远端调控元件（enhancer）经常距离基因几十kb甚至几Mb，但因为染色质三维折叠，它们依然能调控基因表达。所以注释时不要只盯着邻近基因，建议结合H3K27ac或H3K4me1的数据，或者用HiChIP/interaction data来判断远距离调控连接。如果没有这类数据，至少要做一个“peak与差异表达基因关联分析”，取最近基因只是起点，不是终点。

差异peak分析是很多ATAC-seq项目的主要目标：比较两个条件下开放染色质有哪些变化。主流做法是先用DiffBind整理样本的peak矩阵，再用edgeR或DESeq2做差异检验。我的流程如下：

library(DiffBind) sample_sheet <- data.frame( SampleID = c("Ctrl1", "Ctrl2", "Treat1", "Treat2"), Condition = c("Ctrl", "Ctrl", "Treat", "Treat"), Replicate = c(1,2,1,2), bamReads = c("Ctrl1.final.bam", "Ctrl2.final.bam", "Treat1.final.bam", "Treat2.final.bam"), Peaks = c("Ctrl1_peaks.narrowPeak", "Ctrl2_peaks.narrowPeak", "Treat1_peaks.narrowPeak", "Treat2_peaks.narrowPeak") ) dba <- dba(sampleSheet = sample_sheet) 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)

这里有几点实操心得：

第一，DiffBind默认的count窗口以overlap为基础，但ATAC-seq的peak在不同样本间位置会有微小漂移，我建议bUseSummarizeOverlaps = TRUE，这样可以避免边界抖动带来的计数误差。

第二，差异分析前的标准化非常关键。DiffBind默认会用总reads数做TMM标准化，但如果样本间文库复杂度差异很大，我建议额外设置bFullLibrarySize = TRUE，用有效reads总数而不是原始reads数来标准化。这样做能避免线粒体reads占比不同带来的偏差。

第三，差异peak分析至少要两个生物学重复，没有重复就不应该做统计推断。如果只有单样本，建议只看峰的有无、强度变化趋势，别写p值。

3.4 Motif分析与可视化：推测哪些转录因子在干活

开放染色质区域之所以开放，通常是为了让转录因子结合到顺式元件上。ATAC-seq不能直接告诉你是哪个TF结合在peak上，只能通过序列motif来推测。最常用的工具是HOMER：

findMotifsGenome.pl sample_peaks.narrowPeak hg38 motif_output -size 200 -len 8,10,12

HOMER会自动对peak中心区域做de novo motif发现，同时再用已知motif数据库扫描。输出结果里有knownResults.html、homerMotifs.all.motifs等文件。我拿到结果后会优先看排名靠前的known motif的z-score和p-value，z-score超过10的一般比较可靠。

不过motif分析容易出假阳性：一个GC-rich的motif可能在很多peak中都富集，但不代表那个TF真的结合了。我建议把motif富集结果和RNA-seq或蛋白表达数据联合分析，比如motif富集到CTCF，但你的细胞里CTCF根本不表达，那大概率是背景噪音。另一种常用策略是用ATAC-seq信号在特定TF的已知ChIP-seq peak处做富集验证。

可视化这块，我推荐deeptools的computeMatrix和plotHeatmap，在TSS或者peak中心做信号profile：

computeMatrix reference-point \ -S sample.bw \ -R sample_peaks.sorted.bed \ --referencePoint center \ -b 2000 -a 2000 \ -bs 50 \ -p 16 \ -o matrix_sample.gz plotHeatmap -m matrix_sample.gz -out sample_heatmap.pdf \ --sortRegions descend \ --clusterUsingSamples 1

这里input的sample.bw需要先从bam转换成bigwig，我一般用bamCoverage：

bamCoverage -b sample.final.sorted.bam -o sample.bw \ --normalizeUsing CPM --binSize 25 --smoothLength 75 -p 16

之所以用CPM标准化，是为了让不同样本的bigwig可以直接比较。绘图时看到一个清晰的以peak中心为原点、两侧信号递减的模式，说明peak中间是真正的Tn5富集区域，这个结果放文章里也好看。

4. 实战中的质量评估与常见问题排查

4.1 建立一套自己的样本验收标准

分析跑多了就会发现，与其在最后面对一堆峰值找问题，不如一开始就建立一套样本验收标准，把不合格的样本尽早拦截下来。我在实际项目里使用这套验收清单，你可以直接抄来用：

这套标准不是死的。比如线粒体占比这项，我做过一批细胞分选后的样本，线粒体占比普遍超过60%，但核基因组数据质量很好，最终分析结果也可靠。所以遇到个别指标不合格时，要结合多个指标一起判断，别因为一两个数值异常就直接否掉整个样本。

另外，我强烈建议在正式分析前跑一两个样本作为预实验，确认流程能跑通、peak数量合理、差异方向符合预期，再批量跑其他样本。这样可以避免后期发现系统性错误，导致所有样本都要重跑。

4.2 常见报错和排查思路实录

我把自己这些年跑ATAC-seq流程踩过的坑整理成了一份速查表，都是实操中确实遇到的问题：

报错1：Bowtie2比对率极低，只有20%-40%

这个容易让人慌，但先别急着怀疑数据。我遇到过的原因主要有三个：参考基因组物种选错了（比如把小鼠样本比对到人类基因组）、FASTQ文件前后端标识错乱（R1和R2调换）、reads含有大量接头未切除。排查顺序：先用head看看FASTQ的序列长度和碱基组成，再检查参考基因组版本，最后用Kraken2或者比对到rRNA序列统计污染情况。

报错2：MACS2报错“can not find uniquely mapped paired-end reads”

多半是BAM文件中没有配对信息，或者所有reads都被前面的过滤步骤去掉了。检查方式是：

samtools view -c -f 2 sample.final.sorted.bam samtools view -c -f 12 sample.final.sorted.bam

第一个统计正确配对的reads数，第二个统计未配对的reads数。如果未配对reads数占总数的比例很高，回头看排序和过滤步骤是不是弄乱了配对关系，或者bedtools subtract输出后没重新排序。

报错3：peak calling结果出现超大量peak，动辄几十万个

绝大多数情况是Tn5过切导致信号弥散，或者样本之间相互污染。先看插入片段分布是否还呈现核小体周期，再看TSS富集分数是否下降。如果这些指标没问题，可以尝试把MACS2的q-value阈值从0.05收紧到0.01，或者尝试用Genrich重新calling一次做对比。

报错4：两个生物学重复的peak交集不足50%

先别急着怀疑差异分析，先看是否存在批次效应：两个重复是不是在同一次建库中完成？线粒体reads占比有没有明显差异？总reads数差异大吗？我建议用deeptools的plotCorrelation检查重复相关性，并且画PCA图观察样本聚类情况。如果重复确实分得很开，考虑用limma的removeBatchEffect做批次校正，但前提是你要记录好批次信息。

报错5：TSS富集分数只有2左右

说明文库中真正的开放染色质信号很少，噪声占主导。这个情况在建库环节很难挽救，分析环节也不要强行用MACS2出peak。我遇到过TSS富集分数只有1.8的样本，强行跑出来的peak基本全是假阳性，最终我选择放弃该样本。在我的经验里，与其在后期花大量时间校正一个质量差的样本，不如回头优化建库流程更实际。

4.3 从分析流程走向研究结论：几条实用的解读框架

最后聊一下，分析所有的数据和peak之后，怎么把结果转化成生物学结论。我总结了几条常用的解读框架，供参考：

第一，结合motif分析和基因表达做因果推断。如果发现某处理条件下的peak显著增加，且这些peak中富集到某个TF的motif，而这个TF的靶基因又刚好在RNA-seq中显著上调，就可以形成“TF在条件A下结合增强子，激活下游基因”的假设。这个假设虽然还需要ChIP-qPCR验证，但在前期筛选阶段非常高效。

第二，分析差异peak在基因组分布的偏好。如果差异peak大量落在增强子区域，提示处理条件主要影响了远端调控元件；如果大量差异peak落在启动子区域，则说明转录起始层面的调控更突出。这两种情况对应的后续验证方案不一样，前者更适合做3C/HiChIP，后者更适合做启动子报告基因实验。

第三，警惕线粒体reads污染导致的假差异。我在项目中遇到过一个问题：处理组比对照组线粒体reads占比显著更高，导致核基因组reads变少、peak信号整体降低，最终计算出的差异peak几乎全是“下降”。这是个典型的系统误差，处理组根本不是染色质关闭，只是测序深度被线粒体reads稀释了。排查方法很简单：查看各样本的线粒体reads占比，如果组间差异很大，就要考虑用线粒体占比作为协变量，或者在分析前对样本做下采样统一有效reads数。

第四，不要忽略插入片段长度选项。很多下游分析会区分短片段（无核小体，<100bp）和长片段（单核小体，180-247bp）。短片段对应TF结合位点等开放区域，长片段对应核小体周围区域。在做TF footprint分析或者核小体定位分析时，分别对两类片段提取信号，往往能得到比混合分析更清晰的结果。我在做转录因子结合分析时，都会单独提短片段做一轮MACS2，这样得到的peak更尖锐，更符合TF的直接结合特征。

写在最后的几点经验

ATAC-seq这套流程，本身没有什么不可逾越的难点，但想跑出可信、可复现、对生物学问题有用的结果，细节决定成败。我自己最大的体会是：做生物信息分析，不能只关心“命令能不能跑”，更要时时回头想“这个结果在生物学上是否合理”。每一批数据拿到手，先花时间做质量评估，确认数据靠谱，再往下分析；跑完peak calling之后，先别急着做差异和注释，先去IGV里肉眼看看几个TSS和增强子区域的信号，确认peak不是随机噪音堆出来的。这样虽然多花半小时，但能避免后续所有分析建立在不可靠结果上。

另一个实用的建议是：把你的流程固定在版本控制里。我自己的ATAC-seq流程脚本都放进Git仓库，每次修改都标记版本号。这样几个月后回看某个项目时，能知道做差异分析时用的是哪个MACS2版本、哪些参数，结果遇到审稿人质疑时也能迅速复现。数据分析和湿实验一样，都需要可复现性。

最后再分享一个小技巧：如果你在做细胞类型比较或者发育时间序列的ATAC-seq，建议在正式分析前先做一次全局主成分分析（PCA）。PCA能在不引入任何生物学假设的情况下，帮你快速发现离群样本、批次效应和组间分离趋势。我几乎每个ATAC-seq项目都会先跑这一步，至少在正式差异分析前做到心里有数。分析做完之后，也别忘了把不同的可视化结果一起打包存档，到了写文章或者做补充材料的时候，你会感谢那时候留存了高质量图片的自己。