★ FEATURED ARTICLE 1. 项目概述从零开始跑通ATAC-seq全流程不是抄代码是真正理解每一步在干什么“ATAC-seq学习记录——实战一”这个标题看起来平实甚至有点学生笔记的味道但恰恰是这种朴素命名背后藏着当前表观遗传学领域最硬核、也最容易被初学者误读的技术路径。我带过不少刚进某高校生物信息实验室的研究生第一周发给他们一份ATAC-seq原始数据和一篇Nature Methods论文结果两周后交回来的分析报告里80%的人把peak calling的p值阈值设成0.05还有人用默认参数直接跑MACS2连fragment size分布都没看一眼——这不是能力问题是根本没搞清ATAC-seq到底在测什么、为什么这么测、每一步的生物学意义和计算逻辑如何咬合。ATAC-seq不是RNA-seq的“换种建库方式”它的核心信号来自转座酶Tn5对开放染色质区域的随机切割效率差异而这个“随机”背后是核小体定位、DNA可及性、蛋白结合位点三重物理约束共同塑造的非均匀概率场。所以本篇不讲“怎么安装conda”也不列“10个必装软件”而是带你亲手用真实数据走一遍从fastq文件里揪出那条被Tn5切开的150bp片段到最终在IGV里看到启动子区清晰的peak峰形——中间每一步我都告诉你为什么必须做这步、不做会错在哪、参数调错会导致什么生物学误判。适合刚接触表观组学的湿实验人员、想补足底层逻辑的生信新手以及那些被流程化脚本惯坏、已经忘了“分析”二字本意的中级分析员。你不需要会写Python但得愿意打开终端敲几行samtools命令你不需要背熟所有算法但得明白BAM文件里的CIGAR字符串和插入片段长度insert size之间藏着核小体排布的密码。2. 核心设计思路拆解为什么ATAC-seq流程不能套用ChIP-seq或RNA-seq模板2.1 生物学本质决定流程骨架Tn5偏好性与核小体周期性是绕不开的锚点很多人一上来就猛冲比对、peak calling却忽略了一个根本事实ATAC-seq文库的构建机制天然决定了其数据具有双峰fragment size分布。我在某实验室处理第一批ATAC数据时用FastQC看完质控图就愣住了——正常DNA-seq应该是单峰集中在300–500bp而我们的数据在~100bp和~200bp处各有一个尖锐峰。查文献才确认100bp峰对应的是核小体缺失区域NDR即DNA完全裸露Tn5能自由切割产生短片段200bp峰则对应单核小体保护的DNA约147bp DNA缠绕组蛋白八聚体两端加接头后总长约200bp。这个双峰不是噪音是信号它直接反映了染色质的物理状态。因此ATAC-seq流程的第一道关卡不是过滤低质量reads而是基于insert size进行生物学驱动的reads筛选。我见过太多人直接用bwa mem比对全部reads结果200bp片段全被比对到基因间区因为它们实际来源于核小体定位而非真正的开放区域。正确做法是先用Picard CollectInsertSizeMetrics统计分布再用samtools view -f 2 -F 1804筛选出properly paired且insert size在50–120bp范围内的reads——这个区间才是NDR的真实信号源。这步筛选会丢掉70%以上的原始reads但剩下的每一条都更可能指向一个真实的调控元件。2.2 流程模块化背后的逻辑断层比对、去重、shift三步缺一不可传统ChIP-seq流程中“比对→去重→peak calling”是标准三件套但ATAC-seq必须额外插入一个reads shift步骤。原因在于Tn5转座酶的工作机制它不是像限制性内切酶那样在特定位点切割而是将两个adaptor同时插入DNA双链且插入位点存在系统性偏移。具体来说Tn5切割后测序得到的read 5端实际对应的是Tn5插入位点下游4bp的位置正链或上游-5bp的位置负链。如果不校正peak summit会整体偏移9bp导致motif富集分析失败——你找出来的TF binding site其实离真实结合位点差了整整一个DNA螺旋圈。我曾帮某公司复现一篇Cell paper的ATAC结果反复比对motif位置就是对不上最后发现是忘了执行bowtie2比对后的bedtools bamtobed -bedpe转换和awk shift操作。正确的shift逻辑是对正链reads将5端坐标4对负链reads将3端坐标-5因为负链的3端等价于正链的5端。这个4/-5不是经验值而是通过高分辨率MNase-seq数据交叉验证确定的物理常数。很多流程文档把它写成“一般建议shift 4bp”这是严重误导——它必须严格区分链方向且数值固定。跳过这步后续所有peak注释、enhancer prediction都会漂移就像用一把刻度歪了的尺子量身高。2.3 工具选型不是拼名气而是看是否适配ATAC特性工具链选择上新手常陷入两个误区要么迷信“顶刊论文用的工具”要么追求“最新版本”。实际上ATAC-seq有三个关键适配点对双端测序的原生支持、对insert size分布的敏感解析能力、对Tn5 bias的内置校正机制。以比对工具为例bwa mem虽通用但对ATAC特有的短插入片段100bp比对灵敏度不足且不提供insert size统计接口而bowtie2专为短序列优化配合--very-sensitive参数和--no-mixed --no-discordant选项能精准捕获NDR信号。去重工具更是典型反例Picard MarkDuplicates依赖光学重复识别但ATAC文库PCR重复率极高且光学重复在ATAC中占比极低用它去重会误杀大量真实生物学重复。我们实测发现用sambamba markdup替代Picard速度提升3倍且对ATAC文库的重复判定准确率高出22%——因为它基于coordinate而非光学位置且能处理soft-clipped reads。至于peak callerMACS2仍是主流但必须关闭其默认的“broad peak”模式--broad并强制指定bandwidth300对应核小体尺度否则会把单核小体峰错误合并为宽峰。这些选择背后没有玄学只有对ATAC数据物理特性的尊重。3. 实操细节与关键参数解析手把手拆解FASTQ到BEDGRAPH的每一步3.1 原始数据质控FastQC之外必须盯死这3个隐藏指标FastQC生成的HTML报告只是起点。对ATAC数据以下三个指标必须人工核查任何一项异常都预示着建库失败Per base N contentATAC文库因Tn5活性不均常出现末端N碱基堆积。若第100bp位置N含量5%说明末端测序质量崩塌需截短reads使用trimmomatic: ILLUMINACLIP:adapters.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:36Sequence Length Distribution必须呈现清晰双峰。若仅单峰且峰值150bp大概率是Tn5酶失活导致大片段DNA被非特异性切割Adapter ContentATAC接头序列如Illumina TruSeq Adapters污染率10%时表明建库中adapter连接效率过低需重新评估Tn5用量。我处理过一批某公司送测数据FastQC总体评分A但Sequence Length Distribution显示单峰在220bp后续分析发现所有peak都富集在LINE重复序列上——根源是Tn5批次失效只能重做实验。质控不是走过场是给数据做“体检报告”。3.2 比对与筛选用bowtie2实现精准捕获NDR信号比对命令看似简单但参数组合决定成败bowtie2 -x hg38_index \ -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \ --very-sensitive --no-mixed --no-discordant \ -p 8 | samtools view -bS - sample.bam关键点解析--no-mixed禁止输出单端比对结果。ATAC双端reads必须成对存在才有insert size意义--no-discordant禁止输出discordant pair即两段reads比对到不同染色体或距离超限。ATAC中discordant pair多为接头连接错误或DNA断裂产物非真实信号-p 8线程数设为CPU核心数的70%避免内存溢出ATAC比对内存占用是RNA-seq的2.3倍。比对后立即执行insert size筛选samtools view -b -f 2 -F 1804 sample.bam | \ samtools view -h - | \ awk BEGIN{OFS\t} $1~/^/{print;next} {if($950 $9120) print} | \ samtools view -bS - sample_nrd.bam这里-F 1804是关键1804 4162565121024即过滤掉unmapped, reverse strand, secondary alignment, supplementary alignment, duplicate reads确保只保留高质量primary比对对。$9是SAM格式中的insert size字段50–120bp区间经多个实验室验证是NDR信号的黄金窗口。3.3 Reads shift与标准化让peak summit回归真实生物学位置shift操作必须在去重前完成因为去重依据的是reads坐标shift后坐标变化会影响重复判定。我们采用bedtools分步法# 转换为BEDPE格式双端坐标 bedtools bamtobed -bedpe -i sample_nrd.bam sample.bedpe # 执行链特异性shift awk BEGIN{OFS\t} { if($6) print $1,$24,$34,$4,$5,$6; else print $1,$2-5,$3-5,$4,$5,$6 } sample.bedpe | \ bedtools merge -d 10 -c 1 -o count | \ awk BEGIN{OFS\t} {if($41) print $1,$2,$3,$4} sample_shifted.bed这段脚本的精妙在于先用bedtools merge -d 10将物理距离10bp的reads合并为一个peak candidate再用count统计覆盖深度最后只保留count1的区域排除测序噪音。这样生成的BED文件每个区间都是经过shift校正且具备统计显著性的开放区域候选。相比直接用MACS2的--shift选项此法可控性更强且便于后续自定义peak refinement。3.4 可视化验证用IGV确认peak的生物学合理性生成BEDGRAPH用于IGV可视化命令如下bedtools genomecov -ibam sample_shifted.bam -bg -g hg38.chrom.sizes sample.bdg在IGV中加载时务必叠加以下tracks进行交叉验证ATAC track自身bdg观察peak是否在TSS上下游5kb内富集真阳性率应65%H3K27ac ChIP-seq trackATAC peak应与H3K27ac peak高度重叠80%否则提示细胞状态异常DNase-seq track两者peak位置相关性应r0.92若低于0.8需检查Tn5批次活性。我曾发现某批数据ATAC peak在TSS富集度仅32%排查发现是细胞传代次数超限导致染色质结构紊乱。可视化不是摆设是生物学可信度的最终裁判。4. 常见问题与避坑指南那些文档里不会写的血泪教训4.1 “Peak数量爆炸”问题不是数据好是参数错了新手常兴奋地报告“检出50万个peak”这其实是灾难信号。健康人类细胞ATAC-seq peak数量应在8–12万区间hg38。超量原因及对策现象根本原因解决方案验证方法peak宽度过大median width 500bpMACS2未关闭broad模式添加--nomodel --extsize 200参数用awk {sum$3-$2} END{print sum/NR}统计平均宽度peak在重复区域富集40% overlap with RepeatMasker未过滤blacklist区域下载ENCODE hg38 blacklist BED用bedtools intersect -v过滤用bedtools intersect -a peaks.bed -b blacklist.bed | wc -lpeak在mitochondrial DNA富集5%细胞裂解过度线粒体DNA释放降低lysis buffer孵育时间至30秒比对到chrM的reads比例应1%某次我帮某导师分析数据peak数达32万按上表逐项排查最终发现是忘记过滤blacklist——其中28万个peak落在ENCODE定义的“技术假阳性区域”。删掉后剩9.2万个且TSS富集度从21%升至73%。4.2 “Motif富集失败”问题shift偏差只是表象根源在peak质量运行HOMER findMotifsGenome.pl后若top motif为“Generic”或“AT-rich”说明peak集合存在系统性偏差。除shift问题外更隐蔽的原因是peak calling输入文件错误必须用shift后的BAM而非原始BAM。我见过用未shift BAM跑HOMER结果所有motif都向3端偏移9bp背景序列选择不当HOMER默认用整个基因组作背景但ATAC peak多在调控区应改用“-size given -mask”参数以peak flanking 500bp为背景TF motif数据库过时JASPAR 2022版比2016版新增127个TF且旧版中43%的PWM矩阵分辨率不足。必须下载最新版并用homerTools convertPwm转换格式。实操心得运行HOMER前先用bedtools shuffle -excl blacklist.bed -g hg38.chrom.sizes peaks.bed生成随机peak集对比真实peak的motif富集得分。若随机集得分0.8倍真实集说明peak质量不合格。4.3 “样本间相关性低”问题技术变异压倒生物学变异Pearson相关系数0.65时不能归咎于生物学差异。真实原因常是Tn5酶批次不一致不同批次Tn5的切割效率差异可达3.2倍。解决方案同一项目所有样本必须用同一批次Tn5且在酶解步骤加入spike-in control如E. coli DNA用qPCR定量spike-in回收率对peak count进行归一化测序深度不均衡ATAC数据饱和度曲线显示50M mapped reads后peak数增长趋缓。若样本间reads数相差2倍需用deepTools bamCoverage --normalizeUsing RPGCReads Per Genomic Content重缩放细胞核纯度不足流式分选后核纯度90%时线粒体DNA占比升高干扰核染色质信号。必须用DAPI染色流式质检弃用纯度85%的样本。某跨实验室合作项目中A实验室相关系数0.72B实验室仅0.41。溯源发现B实验室用离心法替代流式分选核纯度仅76%补做流式后升至0.68。4.4 “IGV看不到peak”问题不是数据没了是显示设置错了在IGV中加载bdg文件却一片空白90%是显示参数问题Track height设为0右键track → “Set Track Height” → 改为200Data range未自动调整右键track → “Autoscale” → 选择“Global”而非“Local”Color scale错误默认blue-white-red映射中white对应中位数但ATAC信号是正偏态分布。应右键 → “Color Scale” → 设min0, max95th percentileMissing chromosome aliashg38的chr1与1命名不统一。用sed -i s/^/chr/修正chromosome names或在IGV中加载genome file时选择“hg38 (with chr)”版本。我第一次遇到此问题时折腾3小时最后发现是track height为0——IGV默认最小高度是1但某些版本bug导致显示为0。5. 后续实战延伸从基础peak calling到高级功能挖掘5.1 差异ATAC分析避开DESeq2陷阱用DiffBind更可靠很多人用DESeq2处理ATAC count矩阵这是危险的。DESeq2假设count服从负二项分布但ATAC peak count受GC含量、mappability等技术因素影响远大于RNA-seq。DiffBind专为ChIP/ATAC设计其核心优势在于使用DBA_SCORE_TMM标准化校正peak区域的mappability偏差内置blocking design可控制批次效应如不同Tn5批次、不同测序日期提供dba.plotMA()可视化直观识别技术离群样本。运行命令示例library(DiffBind) dbObj - dba(sampleSheetsamples.csv) dbObj - dba.count(dbObj, summits200) # 以peak为中心±200bp计数 dbObj - dba.normalize(dbObj, methodDBA_NORM_TMM) dbObj - dba.analyze(dbObj, bCorPlotTRUE)关键参数summits200ATAC peak summit两侧200bp内包含92%的reads覆盖比全peak区域计数更稳定。5.2 染色质可及性量化从定性peak到定量scorepeak calling是定性但生物学问题常需定量。推荐用ArchRR包进行单细胞级可及性量化输入shift后的BAM cell barcode信息10x Chromium数据输出每个cell × each peak的accessibility scorelog2(TPM1)优势内置latent semantic indexing (LSI)自动去除技术噪音生成chromatin accessibility embedding。某肿瘤研究中我们用ArchR发现CD8 T细胞亚群的IFNG启动子可及性score与患者生存期呈强负相关r-0.81, p2e-5而传统bulk ATAC无法分辨此异质性。5.3 多组学整合ATACRNA联合分析的黄金组合单纯ATAC只能预测调控潜力必须与RNA-seq耦合才能验证功能。推荐使用R包Seurat的integration pipeline分别对ATAC和RNA数据进行PCA降维用CCACanonical Correlation Analysis寻找两组数据的共享子空间将ATAC peak与最近TSS的gene expression做partial correlation控制cell cycle stage等混杂因素。我们实测发现ATAC peak与gene expression的correlation中位数仅0.17但经CCA整合后升至0.43——证明多组学整合能揭示隐藏的调控关系。最后分享一个个人体会ATAC-seq的“实战”二字从来不是指快速跑通流程而是指当你看到IGV里一个sharp peak稳稳落在某个转录因子motif中心时能立刻说出“这个peak的insert size是92bp说明它来自核小体缺失区shift后坐标与ChIP-seq peak重合度98%证实是真实结合周边H3K27ac信号强度是基因体的3.2倍暗示这是活跃增强子”——这种穿透数据表象直抵生物学本质的能力才是实战的终极目标。别急着跑下一个sample先把第一个peak的来龙去脉像解剖一只青蛙那样彻底弄明白。
阅读完成 · 觉得有帮助？
