ATAC-seq 入门（一）：先设计比较，再看质量控制

这是一份面向 bulk、双端测序（paired-end）ATAC-seq 的学习教程，不适用于把单细胞 ATAC-seq 按同一流程直接处理。文中的样本名称和问题场景均为教学假设，不代表本站作者完成过相应实验，也没有真实样本的分析结果。

本系列采用已核验的 nf-core/atacseq 2.1.2 作为流程学习入口。第一篇不急着运行峰识别，而是回答一个更早的问题：手里的数据，是否真的能够支持想做的比较？

一、先把“开放”说准确

ATAC-seq 利用转座酶对相对可及的染色质区域进行标记，随后从测序片段的分布观察染色质可及性。它提供的是一类测量信号，不是所有调控结论的直接答案。原始方法论文是理解实验原理的合适起点，而不是一份可以无条件照搬到所有组织的参数表。原始 ATAC-seq 方法论文

例如，某个区域在处理组中信号增高，可以暂时表述为“该区域的相对可及性发生变化”。此时还不能直接写成“某转录因子结合增强”，更不能凭最近的基因名称写成“这个基因导致了表型”。测量、统计比较、机制解释，是三层不同的问题。

开始前，建议先写四句话：比较的细胞或组织是什么；处理条件是什么；独立样本来自哪里；准备控制哪些已知混杂因素。如果这四句话写不清，后面的软件即使正常结束，也可能只得到一个漂亮但难以解释的结果目录。

二、生物学重复不是多几个 FASTQ

假设准备比较 CONTROL 与 TREATMENT，各有三个独立培养并单独制备文库的样本。这是每组三个生物学重复。如果同一个文库被分到两个测序 lane，两个 lane 增加的是技术测量，不是独立样本数。把它们当成两个重复，会夸大统计上的信息量。

在固定的 2.1.2 版本里，同一实验组使用相同 sample 标签，以 replicate 区分生物学重复；同一生物学样本的重复测序使用相同的两个标签。务必按该版本的输入语义填写，不要凭“每行就是一个独立样本”的直觉处理。nf-core 2.1.2 输入说明

下面是结构示例，路径只是占位符，必须换成实际可访问的文件：

1
2
3
4
5
6
7
sample,fastq_1,fastq_2,replicate
CONTROL,/data/C1_R1.fastq.gz,/data/C1_R2.fastq.gz,1
CONTROL,/data/C2_R1.fastq.gz,/data/C2_R2.fastq.gz,2
CONTROL,/data/C3_R1.fastq.gz,/data/C3_R2.fastq.gz,3
TREATMENT,/data/T1_R1.fastq.gz,/data/T1_R2.fastq.gz,1
TREATMENT,/data/T2_R1.fastq.gz,/data/T2_R2.fastq.gz,2
TREATMENT,/data/T3_R1.fastq.gz,/data/T3_R2.fastq.gz,3

流程输入表之外，还需要独立的实验元数据表。建议保存 sample_id、condition、batch、donor、library_id、lane、protocol、read_length，并注明哪些字段未知。sample_id 可以采用流程输出的 CONTROL_REP1 等标识，但不能让一个样本在计数表、实验记录和质控表里出现三个不一致的名字。

三个重复只是这个教学设计的选择，不是任何项目都足够的承诺。样本异质性、预期效应和研究目标都会影响所需的独立样本数。ENCODE 的标准强调生物学重复与重复间一致性；这些要求比“把深度补足就行”更接近比较问题的核心。ENCODE ATAC-seq 数据标准

三、批次问题必须在比较前暴露

考虑两个假设设计。设计甲把所有对照在第一天建库、所有处理在第二天建库；设计乙在两天里都安排对照和处理。甲的处理条件与建库批次完全重合。发现组间差异后，无法仅凭现有数据判断来自处理还是建库日。把 batch 加进公式，不能凭空补出缺失的比较信息。

设计乙至少提供了同批次内的组间比较机会，但依然需要记录供体、细胞比例、建库步骤等因素。批次交叉不是免检证书。尤其是 bulk 数据，细胞组成变化也可能改变可及性图谱，不能将这种混合信号自动解释为同一种细胞内部的调控变化。

下面的标准库练习只使用内存中的虚构元数据，不读取测序文件，也不需要安装生信软件。它检查样本标识、重复数与批次交叉情况：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from collections import Counter, defaultdict

rows = [
    ("CONTROL_REP1", "control", "day1"),
    ("CONTROL_REP2", "control", "day2"),
    ("CONTROL_REP3", "control", "day1"),
    ("TREATMENT_REP1", "treatment", "day1"),
    ("TREATMENT_REP2", "treatment", "day2"),
    ("TREATMENT_REP3", "treatment", "day2"),
]
ids = [sid for sid, _, _ in rows]
assert len(ids) == len(set(ids)), "Duplicate biological sample ID"
counts = Counter(condition for _, condition, _ in rows)
assert set(counts) == {"control", "treatment"}, "Both comparison conditions are required"
assert all(n >= 2 for n in counts.values()), "No replicated comparison"
by_batch = defaultdict(set)
for _, condition, batch in rows:
    by_batch[batch].add(condition)
for batch, conditions in sorted(by_batch.items()):
    print(batch, sorted(conditions))
if all(len(conditions) == 1 for conditions in by_batch.values()):
    raise ValueError("Condition may be completely confounded with batch")

把所有对照改成 day1、所有处理改成 day2，观察程序为什么停止。这个简化检查只能提示一种风险，不能代替第三篇中的设计矩阵满秩检查，更不能替代实验设计审查。

四、参考文件是一套坐标系统

至少需要匹配的基因组 FASTA 与基因注释，同时准备适用的 blacklist 和用于 TSS 统计的注释。不要把“人类参考”当成足够明确的版本说明：GRCh37 与 GRCh38 的坐标不能直接混用；同一组装的不同补丁、附加 contig 和染色体命名也需要核对。

建议把参考包写成一张清单：组装名称、FASTA 来源和校验值、GTF 发布版本、blacklist 的版本和适用组装、TSS 文件的生成规则、线粒体 contig 名称。chr1 与 1 不一致，可能让区域相交悄悄变成空集；线粒体叫 chrM 还是 MT，应从实际 FASTA 中确认，而不是凭记忆填。

blacklist 的作用是标记容易产生异常富集或不可靠解释的区域，不是删除“不符合预期”的结果。TSS 文件则定义统计中心的位置：基因、转录本的取舍以及重复 TSS 的处理，都可能影响富集曲线。换了参考定义，旧报告里的阈值就未必能原样比较。流程提供相关参考参数，但参数存在不等于参考文件已经匹配。nf-core 2.1.2 参考参数

五、质控是一条证据链，不是一个总分

先看原始 reads 的质量分布、接头和 R1/R2 完整性，再看比对与过滤损失，最后看信号结构和重复关系。高比对率不必然代表好文库：大量线粒体片段、重复片段或者背景片段也可以成功比对。只报告最终百分比，会隐藏信息在什么步骤丢失。

建议分层记录以下内容，并为每个数字写清输入文件和分母：

检查层 需要观察的内容 不应直接推出的结论
原始数据 接头、质量、读段数、配对是否完整 reads 多就有统计功效
比对与过滤 核基因组映射、重复、线粒体、有效片段 比对率高就能进入差异分析
信号结构 TSS 曲线、片段长度分布、峰中富集 单个指标过线就是合格
重复关系 相似度、异常样本、批次聚类、峰一致性 聚类符合预期就没有混杂

FRiP 需要特别小心。“峰内 reads / 可用 mapped reads”与“至少命中一个峰的独立片段 / 可用独立片段”不是同一个统计口径。还要说明峰集合来自单样本还是统一集合、是否去重复、是否去线粒体、采用何种重叠规则。若两个报告的分母不同，不能简单用数值大小评价文库优劣。TSS enrichment 也依赖所用注释、窗口和背景归一化定义。ENCODE 质控术语

ENCODE 给出了特定数据标准下的深度、FRiP、TSS 及重复一致性要求，可以作为参照，而不是所有样本的通用裁决器。组织、协议、参考注释及分析口径不同，阈值的适用性需要另行评估；既不能把边界值当万能通行证，也不能任意降低标准来保住想要的结果。

六、假设异常样本，怎样处理才有依据

设想一个处理组重复线粒体比例偏高、TSS 曲线平坦，同时与同组另外两个样本不一致。第一步不是删掉它，而是查配对关系、参考命名、实验记录和过滤日志。可以先排除文件错配、文库标识错置或参数差异。若问题仍在，记录哪些证据支持其质量不足，再按事先约定的规则作决定。

不能因为它使差异峰变少而删除，也不能只因 PCA 上距离远就断言它失败。如果剔除后某组不再有可用生物学重复，或者剩余样本令条件与批次完全混杂，就应停止该差异比较。还能保留描述性质控和探索图，但不要发布伪装成可靠推断的显著峰列表。

七、这一篇应留下什么

第一篇的输出不是峰，而是可审查的分析起点：输入文件清单、独立样本元数据、参考包清单、质控口径表和预先记录的继续／暂停条件。也应写下尚不知道的事项，例如供体关系未确认、某样本建库日期缺失，而不是用猜测填满表格。

进入下一篇前，自问：我能解释生物学重复与 lane 的区别吗？每个样本的 R1/R2 是否对应？参考文件是否采用同一坐标？各组是否有足够的可用独立样本？批次是否允许区分条件效应？FRiP 和 TSS 的定义是否已记录？如果任一关键问题未解决，先补证据，再运行更复杂的统计。

下一篇：从 FASTQ 到可审计的峰集合。最后一篇会把“统一区域、统一计数、明确模型”连起来：差异可及性不是比较峰数量。