上一章得到峰集合，并不意味着差异分析已经完成。“对照没有叫到峰，处理叫到了峰”是峰识别结果的比较；“在考虑深度、批次与重复变异后，这个区域的信号是否有统计支持的变化”才是差异可及性的问题。前者可能由测序深度或阈值差异造成，不能代替后者。

本文继续讨论 bulk paired-end ATAC-seq，沿用两组各三个独立生物学重复、非配对的学习设计。唯一可以直接运行的小练习使用人工数据，只检查结构和计数约定；真实 BAM 计数及 DESeq2 部分是资料准备完成后的模板，本文没有运行真实测序分析，不提供杜撰的差异区域数量、p 值或图。

一、先过 QC 门，再问差异

统计软件不会因为样本被错标、两组使用不同过滤标准而自动阻止分析。开始之前，要回看比对、重复率、线粒体比例、片段长度结构、TSS 富集、FRiP 与重复一致性，并确认每项指标使用的输入和单位。阈值需要结合协议、组织与实验背景；不能把一个固定数当作所有文库的通行证。

QC 失败、身份存疑或有效生物学重复不足时，先调查，必要时只保留描述性结果。不要为了让 PCA 分组漂亮而删样本；排除应有独立的技术证据、记录，并重新检查排除后的设计。把同一个样本的 lane、R1/R2 或 pooled BAM 当作新增重复，不能增加独立实验单位。

本例的三个重复是教学约定，不是充分功效保证。真实设计还要根据变异、目标效应和研究目的考虑样本量；只有一个来源每组一次时，不应套模板产生看似正式的显著性结论。DiffBind 官方教程将计数、质量检查、设计与归一化分开讨论，也适用于理解 ATAC 数据的这些边界。

二、每个样本必须站在同一张区域地图上

比较矩阵的行应来自一个固定的候选区域集合，列对应所有保留的独立样本。先根据明确规则建立 consensus regions，再在这些相同坐标上对所有样本计数。某个样本没有在该处叫到峰，也要在同一区域读取它的计数，而不是将该格直接填成零或缺失。

沿用 nf-core/atacseq 2.1.2 的 Bowtie2、narrow peaks 学习路线，可在 bowtie2/merged_library/macs2/narrow_peak/consensus/ 查看统一区域及相关计数材料。但这里的流程输出是起点，不是对所有实验设计的完整统计承诺。固定版本输出说明

要保存候选集合的生成规则、原始峰来源、重复支持要求和最终 BED 校验值。合并区间会改变宽度和计数机会；改成固定宽度的 summit 窗口则是另一个分析定义，需要重新计数并记录。不能让处理组用自己的峰宽、对照组用另一套峰宽，再把两个表按最近位置拼起来。

参考 FASTA、BAM header、BED、blacklist 和后续注释必须属于同一组装并使用一致的 contig 名称。chr1 与 1 不匹配不是小排版问题。BED 的零起点半开区间与 SAF 的坐标约定也不能混用；转换时保留 region ID，并用已知小区间验证边界。blacklist、线粒体与低质量比对过滤应采用同一规则，记录删除了什么；不能只在某一组移除可疑区域。

三、一条双端片段不能被 R1 和 R2 算成两次

本章选择片段级计数：同一有效 read pair 在一个区域内贡献一次，不将两端各算一个观察。不要将切点数、覆盖碱基总和、CPM、bigWig 信号或峰的富集分数混称为原始片段计数。沿用上一章的片段主线，不再额外给 BAM 叠加 Tn5 shift；另做切点或 footprint 分析时使用单独派生文件及明确校正记录。

“计数为整数”只是必要检查，不能证明计数单位正确。一个长片段可能命中多个区域；应在分析前决定是排除歧义分配、唯一分配，还是采用明确记录的多区域策略。本章的严格教学约定是：排除歧义分配，保留的片段只归一个 region，从而不会跨行重复累计。真实方法可以不同，但所有样本必须一致，不能临时为某个区域选择最有利的规则。

特别注意固定流程配置中的 -O --fracOverlap 0.2：-O 允许多区域分配，--fracOverlap 是重叠比例要求，不等于分数权重开关。两者不能证明一个片段只出现于一行；是否产生分数计数还要核对其他参数与实际输出。不要盲目四舍五入小数，也不要把整数多重分配误称为唯一分配。2.1.2 计数配置源码、featureCounts 作者说明

下面是准备完成后才可运行的重新计数模板，不是复刻上面的流程默认值。它面向已过滤、索引完整且参考一致的双端 BAM；使用支持这些参数的 Subread 版本并记录 featureCounts -v。所有占位输入必须预先核验。-p --countReadPairs 明确片段计数，-B -C 排除指定不完整或异常配对，未加 -O 或 --fraction，歧义多区域分配不保留。实际重叠判据依据 featureCounts 的配对比对语义，不应声称它必然覆盖 insert 内的每个碱基。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# PREPARED INPUTS ONLY. No installation, alignment or peak calling here.
# consensus.saf has one unique region ID per region; coordinates validated.
# All six BAMs are independent samples, already filtered with the same policy.
atac_count_dir="$(mktemp -d "$PWD/atac-count-XXXXXX")" || exit 1
featureCounts -T 4 -F SAF -s 0 \
  -p --countReadPairs -B -C \
  -a /data/project/consensus.saf \
  -o "$atac_count_dir/counts.featureCounts.txt" \
  /data/project/bam/CONTROL_REP1.bam \
  /data/project/bam/CONTROL_REP2.bam \
  /data/project/bam/CONTROL_REP3.bam \
  /data/project/bam/TREATMENT_REP1.bam \
  /data/project/bam/TREATMENT_REP2.bam \
  /data/project/bam/TREATMENT_REP3.bam

这段只为本次计数创建全新目录，不删除任何数据。需检查退出状态、summary 中的成功与失败分配原因，并保留日志及参数。输出含注释列和 BAM 路径列名；导出下游矩阵时只保留计数列，按明确的一对一 BAM→sample ID 对照重命名，不能直接将整张 featureCounts 表当作数值矩阵。Subread 用户手册解释了双端计数、重叠分配和 SAF 输入的细节。

四、可直接运行的小练习：身份、边界和计数审计

以下只需要 Python 3 标准库，不读真实文件、不安装软件、不写磁盘。三个区域和十八个数都是人为构造，只用于验证结构。我们故意打乱元数据的顺序，程序按 sample ID 恢复对应关系。人工片段台账中每个 pair 只分配一个区域；重复写入同一 pair 会被拒绝。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

from collections import Counter, defaultdict

samples = ["CONTROL_REP1", "CONTROL_REP2", "CONTROL_REP3",
           "TREATMENT_REP1", "TREATMENT_REP2", "TREATMENT_REP3"]
conditions = ["control"] * 3 + ["treated"] * 3
batches = ["day1", "day2", "day1", "day1", "day2", "day2"]
metadata = [dict(sample_id=s, biological_id=f"demo_{i + 1}",
                 condition=c, batch=b, qc_pass=True)
            for i, (s, c, b) in enumerate(zip(samples, conditions, batches))]
metadata = list(reversed(metadata))
regions = [("region_01", "chr1", 100, 200),
           ("region_02", "chr1", 400, 500),
           ("region_03", "chr2", 100, 200)]
counts = [[10, 11, 9, 14, 12, 13],
          [20, 18, 21, 19, 22, 20],
          [5, 6, 4, 5, 7, 6]]

assert len(samples) == len(set(samples))
assert len(metadata) == len(samples)
assert len({r["sample_id"] for r in metadata}) == len(metadata)
assert {r["sample_id"] for r in metadata} == set(samples)
by_sample = {r["sample_id"]: r for r in metadata}
aligned = [by_sample[s] for s in samples]
assert [r["sample_id"] for r in aligned] == samples
assert all(r["qc_pass"] is True for r in aligned)
assert len({r["biological_id"] for r in aligned}) == len(aligned)
assert {r["condition"] for r in aligned} == {"control", "treated"}
assert all(n >= 3 for n in Counter(r["condition"] for r in aligned).values())
batch_conditions = defaultdict(set)
for row in aligned:
    batch_conditions[row["batch"]].add(row["condition"])
# This is a conservative teaching check, not a general design-rank test.
assert all(v == {"control", "treated"} for v in batch_conditions.values())

assert len({r[0] for r in regions}) == len(regions)
for region_id, chrom, start, end in regions:
    assert region_id and chrom in {"chr1", "chr2"}
    assert type(start) is int and type(end) is int and 0 <= start < end
for i, a in enumerate(regions):
    for b in regions[i + 1:]:
        assert not (a[1] == b[1] and a[2] < b[3] and b[2] < a[3])
assert len(counts) == len(regions)
assert all(len(row) == len(samples) for row in counts)
assert all(type(n) is int and n >= 0 for row in counts for n in row)

# Synthetic, already uniquely assigned pairs; not a BAM parser.
assignments = []
for (region_id, *_), row in zip(regions, counts):
    for sample, n in zip(samples, row):
        assignments.extend((sample, f"{region_id}_pair_{j}", region_id)
                           for j in range(n))
pair_keys = [(s, pair_id) for s, pair_id, region_id in assignments]
assert len(pair_keys) == len(set(pair_keys)), "Repeated pair assignment"
recounted = Counter((region_id, s) for s, pair_id, region_id in assignments)
assert all(recounted[(region_id, s)] == n
           for (region_id, *_), row in zip(regions, counts)
           for s, n in zip(samples, row))
print("PASS: 3 synthetic regions, 6 independent sample IDs, matched metadata")
print("Synthetic matrix column totals:",
      dict(zip(samples, [sum(row[j] for row in counts) for j in range(6)])))

检查输出应包括 PASS 和列总量 35、35、34、38、41、39。这些数是代码练习答案，不是深度、FRiP 或生物学效应。这里的区间断言也仅检查这三个教学区域，不能代替真实参考、blacklist 与 BAM 的审计。

试着复制一条 assignment，重复片段断言应失败；把一个计数改成 2.5，整数断言应失败；删除一个元数据样本，集合匹配应失败。最后将所有对照改为 day1、处理改为 day2，批次断言应拒绝继续。这些练习帮助区分“程序可以算”与“比较有意义”。

五、设计公式必须对应真实实验单位

本例的 day1 包含 CONTROL_REP1、CONTROL_REP3、TREATMENT_REP1，day2 包含其他三个样本，两组都出现在两个批次。采用 ~ batch + condition 的加性模型，估计考虑批次后的共同处理系数，不估计批次×处理交互。两组各三个不同来源，因此不加入一个每样本唯一的 subject 因子；那会耗尽可估计的信息。

若真实实验是同一供者各有处理与对照，则应携带供者身份并审查 ~ subject + condition 等配对设计，而不是照抄本例。若所有对照来自 day1、所有处理来自 day2，处理与批次完全混杂；删掉 batch 变量只会隐藏问题，不能创造可识别的处理效应。设计矩阵满秩是必要条件，不保证样本量充分、不存在残余混杂或因果解释成立。

下面的 R 模板仅接收已完成上述审计的真实、全量统一区域原始整数片段矩阵。不要把三行人工数据接入它来计算科学显著性。准备 raw_region_counts.tsv（第一列 region ID，其余列为 sample ID）与 metadata.tsv（第一列 sample ID，另有 biological_id、condition、batch、qc_pass；QC 合格写 TRUE）。已有相容 R/Bioconductor 与 DESeq2 时才可运行，不自动安装。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

# PREPARED REAL DATA ONLY: never feed the three-region teaching matrix here.
if (!requireNamespace("DESeq2", quietly = TRUE)) stop("Prepare DESeq2 first")
suppressPackageStartupMessages(library(DESeq2))
cts <- as.matrix(read.delim("raw_region_counts.tsv", row.names = 1,
                           check.names = FALSE))
meta <- read.delim("metadata.tsv", row.names = 1, check.names = FALSE)
stopifnot(is.numeric(cts), all(is.finite(cts)), all(cts >= 0),
          all(cts == floor(cts)), all(cts <= .Machine$integer.max),
          nrow(cts) > 0L, ncol(cts) > 0L,
          all(nzchar(rownames(cts))), all(nzchar(colnames(cts))),
          !anyDuplicated(rownames(cts)), !anyDuplicated(colnames(cts)),
          !anyDuplicated(rownames(meta)),
          setequal(colnames(cts), rownames(meta)),
          all(c("biological_id", "condition", "batch", "qc_pass") %in% names(meta)))
meta <- meta[colnames(cts), , drop = FALSE]
stopifnot(identical(colnames(cts), rownames(meta)),
          !anyNA(meta), all(meta$qc_pass == TRUE),
          all(nzchar(trimws(meta$biological_id))),
          all(nzchar(trimws(meta$batch))),
          !anyDuplicated(meta$biological_id),
          all(meta$condition %in% c("control", "treated")))
meta$condition <- factor(meta$condition, levels = c("control", "treated"))
meta$batch <- factor(meta$batch)
stopifnot(nlevels(meta$batch) >= 2L, all(table(meta$condition) >= 3L),
          all(table(meta$batch, meta$condition) > 0L))
design_formula <- ~ batch + condition
x <- model.matrix(design_formula, data = meta)
stopifnot(qr(x)$rank == ncol(x), nrow(x) > ncol(x))
storage.mode(cts) <- "integer"
dds <- DESeqDataSetFromMatrix(cts, colData = meta, design = design_formula)
keep <- rowSums(counts(dds) >= 10L) >= min(table(meta$condition))
stopifnot(any(keep))
dds <- dds[keep, ]
stopifnot(all(rowSums(counts(dds)) > 0L))

# Run this default normalization only after reviewing section 6's assumptions.
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"), alpha = 0.05)
tab <- data.frame(region_id = rownames(res), as.data.frame(res), row.names = NULL)
out <- tempfile("perry-atac-de-", tmpdir = getwd())
stopifnot(dir.create(out))
write.table(tab, file.path(out, "treated_vs_control_all_regions.tsv"),
            sep = "\t", row.names = FALSE, quote = FALSE)
saveRDS(list(dds = dds, res = res, metadata = meta, design = design_formula,
             contrast = c("condition", "treated", "control")),
        file.path(out, "atac_bundle.rds"))
capture.output(sessionInfo(), file = file.path(out, "sessionInfo.txt"))
message("Prepared-data outputs saved to: ", normalizePath(out))

本模板主动要求交叉批次例子和独立来源，不是所有设计的通用入口。预过滤的计数阈值与最小组样本数是预先声明的教学选择，应结合真实深度重新评估。保存的结果表对应进入模型的全部区域，而不是只保存显著行；原始计数、区域坐标和生成记录仍需另行归档。函数的计数入口、因子水平和显式 contrast 规则可与DESeq2 官方教程逐项核对。

六、归一化不是全局可及性测量仪

原始 counts 与测序深度有关，不能仅凭一列总数较大就宣称染色质整体更开放。常规的相对归一化依赖一组用于确定尺度的特征及其变化分布；大片区域同方向变化、处理引起细胞组成改变或建库效率改变时，尺度选择尤其需要审查。

如果研究问题是每个细胞的全局可及性变化，普通测序数据的相对组成不能自动提供绝对尺度。需要在实验阶段考虑可验证的外部控制、兼容协议的 spike-in、细胞数量与实验效率记录，或有充分依据的稳定参考区域；不能事后指定几个“看起来不变”的峰就声称问题已经解决。背景区域归一化也有自己的假设，不保证所有背景都稳定。DiffBind 提供多种归一化路线，选择工具本身不会消除这些判断。

应在结果解释前声明尺度依据，查看 QC、size factors、MA 图与其他诊断，必要时做有理由的敏感性分析。不能试遍归一化方式，只保留最符合预期的一套；更不能先对 counts 做 CPM，再把它送入原始计数入口。这一节也是上面的默认 DESeq2 模板的使用门槛，不满足时应暂停而不是硬跑。

七、看效应、精度和候选规则，不只看一颗星

本例固定 treated 相对 control 的方向：正 log2FC 表示在声明的模型与归一化尺度下信号更高。效应、标准误、p 值与 FDR 调整后的值回答不同问题；不显著不等于没有变化，低计数的大倍数也不自动值得优先解释。NA 可能与低信息、过滤或异常值有关，应追查原因，不能全部替换成零。

padj < 0.05 是声明模型下的多重检验规则，不是某个区域有“95% 概率真实”。普通检验后再按效应大小筛候选，不等于已经检验了“真实效应超过这个阈值”。可参考RNA-seq 第二章的设计与结果解释，其中统计概念能够借鉴，但 RNA 的计数单位与 ATAC 的归一化问题不能直接互换。

常见误区还有：峰出现／消失被当成定量差异；pooled BAM 被当成多个重复；图上不同自动纵轴被当成效应证据；整数矩阵被误当成正确片段矩阵。这些错误都可能得到漂亮结果，因而需要在结果之前检查，而不是靠图的美观程度判断。

八、从区域到机制，还有几道需要跨越的桥

motif 富集表示某类序列模式相对明确背景更集中，不证明对应转录因子在该条件下实际结合或激活。motif 家族相似、区域长度和 GC 背景都需要考虑；应记录候选与背景，而不是把一个富集名称当作功能结论。

最近基因注释只是位置上的关联，不等于已确认的调控靶点。远端区域可能作用于其他基因；ATAC 与 RNA 的方向一致可以帮助提出假设，却仍可能受到细胞组成、时间和共同上游因素影响。二者叠加并不会自动升级为因果机制，需要更具体的定位、时序、结合或干预证据。

最后的交付应能说清：比较哪几个独立实验单位、在哪套统一区域上按什么规则计数、QC 与归一化为何可信、设计和对比是什么，以及哪些结论只是待验证的解释。可重查的原始矩阵、区域 BED、元数据、参数、版本、完整表和排除记录，比只有一张火山图更有用。

上一章：从 FASTQ 到可审计的峰集合；课程起点：实验设计与质控口径。

官方参考与继续阅读

• nf-core/atacseq 2.1.2：consensus、计数与 QC 输出
• nf-core/atacseq 2.1.2：流程参数配置源码
• Subread：featureCounts 作者说明
• Subread：用户手册及配对、重叠、注释格式约定
• Bioconductor：DiffBind 官方教程，尤其归一化章节
• Bioconductor：DESeq2 官方计数输入、设计与 contrast 教程