BEDTools 实跑复盘，100 碱基的区间为什么会算出 110 碱基重叠

把一组峰与注释求交后，输出行数比输入峰数更多，重叠长度相加甚至超过峰本身的长度。这些现象未必意味着软件出错。它们首先提示我们，正在统计的可能是区间之间的配对关系，而不是被注释覆盖的独立碱基，或者有注释的峰的数量。

本文实际安装并运行 BEDTools 2.31.1，用自行构造的四个查询区间和六个注释区间，依次核对配对、覆盖、比例门槛、链方向、相邻区间合并、转录本分块和排序检查。每一个数字都来自本次输出，读者可以用纸笔复算。这里没有真实测序数据，也没有通过小例子的耗时推测全基因组性能。

一、先写清楚，一行到底代表什么

查询文件中的四个区间分别叫作甲、乙、丙、丁，对应脚本里的 P1 到 P4。它们位于同一条虚构染色体，起止坐标依次是 100 到 200、200 到 300、500 到 600、700 到 800，长度都为一百。甲、丙、丁标记正链，乙标记负链。

注释文件中，前三段为 90 到 140、130 到 180、180 到 200，都与甲相交。随后两段为 200 到 250、250 到 300，都与乙相交。最后一段为 550 到 560，只与丙相交。丁没有对应注释。脚本保留了每个注释的名称和链方向，不需要依靠行号猜测身份。

这里使用 BED 的零起点、右端不包含约定。100 到 200 的长度是终点减起点，等于一百。200 到 250 与它仅边界相接，没有共同覆盖的碱基。本次单独运行边界对照，软件返回的重叠条数确实为零。将两个区间都看作双端包含，会凭空创造一个交点。

实际项目还需要检查参考组装版本、染色体命名和坐标来源。把染色体名称补上前缀，不能把旧版本坐标变成新版本坐标。对由浏览器、注释文件或手工表格导出的区域，最好留下原始格式与转换步骤，再检查一两个已知位置。格式正确与生物学位置正确是两项检查。

二、六行输出，不代表有六个峰得到注释

本次带双方完整记录的求交产生六行，含义是六对查询与注释匹配。甲出现三次，乙出现两次，丙出现一次。输出中重复出现甲，是因为它与三条注释分别建立关系，这三行承担了不同的信息。

如果问题是“哪些查询至少命中过一条注释”，使用只输出命中查询的模式，得到甲、乙、丙三个区间。如果问题是“每个查询各有几条命中”，计数模式给出三、二、一、零，仍然保留四个查询。没有命中的丁也很重要，否则计算注释率时，分母就已经悄悄变成了筛选后的集合。

bedtools intersect -a A.bed -b B.bed -wa -wb
bedtools intersect -a A.bed -b B.bed -u
bedtools intersect -a A.bed -b B.bed -c
bedtools intersect -a A.bed -b B.bed -v

这几条命令使用同样的输入，只改变输出要回答的问题。第一条适合保存注释关联表，第二条适合筛出有命中的查询，第三条适合把命中数量作为查询的一列属性，第四条适合追踪未命中区域。将第一条输出直接按行计算峰数，会把注释结构的复杂度误当成峰数量增加。

还有一个实际跑过的细节。将同一条查询在输入中写两次，再使用只保留命中查询的参数，输出仍然有两条相同记录。这个参数避免的是一条输入因多个匹配而重复输出，不负责把输入文件本来就存在的重复记录合并。清理重复坐标、重复名称和保留不同转录本，应在明确业务含义之后分别处理。

三、110 碱基来自重复累计，覆盖宽度仍是 100

甲与前三条注释的重叠长度，实际输出分别为四十、五十、二十，相加是一百一十。前两条注释在 130 到 140 的十个碱基范围内重合，因此这十个位置被计入了两次。两两求交的长度之和适合描述配对累计量，但不能直接当作甲的非重复覆盖长度。

运行覆盖统计后，甲的非零覆盖碱基数为一百，覆盖比例为百分之百。软件同时报告三条注释命中、一百碱基被覆盖、一百碱基查询长度以及覆盖比例。四个输出量各有含义，不能只因为都来自同一行，就将命中数当深度、将长度当比例。

bedtools intersect -a A.bed -b B.bed -wo
bedtools coverage -a A.bed -b B.bed

前一条便于追溯每一对重叠的来源，后一条回答查询本身有多少位置被覆盖。假如你的研究问题是“启动子被某类状态覆盖的比例”，应优先明确是否需要区间并集意义上的覆盖。如果问题是“每个启动子对应多少条注释记录”，则可以计算匹配数，但要说明同一基因的多个转录本是否算多个记录。

这里也可以形成一个简单的质控规则。对于非重复覆盖比例，结果不应超过一。若出现大于一的数值，先检查是否把每对匹配的重叠长度累加，再除以查询长度。不要急着截断到一，因为截断只会隐藏重复累计，让错误统计口径继续进入下一步。

四、覆盖八成的门槛，到底施加在谁身上

最容易混淆的是单次匹配比例与总覆盖比例。甲的注释并集覆盖百分之百，但任何一条注释对甲的覆盖都没有达到八成。本次将最小比例设为八成后，甲的合格匹配数变成零。乙也一样，两条各覆盖一半的注释合起来覆盖全部，却没有任何单条达到八成。

bedtools intersect -a A.bed -b B.bed -f 0.8 -c
bedtools intersect -a short.bed -b long.bed -f 0.9 -u
bedtools intersect -a short.bed -b long.bed -f 0.9 -r -u

第一条命令在本例返回四个零。它并没有否定前面的百分之百覆盖，而是对每一对候选匹配单独应用门槛。如果需求是注释并集覆盖查询至少八成，应先计算覆盖比例，再按该比例筛选。若业务允许先合并注释区间，也需要确认合并是否会丢掉基因、类别或链方向等需要区分的信息。

后两条用独立的长短区间检验分母。短区间为 100 到 200，长区间为 0 到 1000。短区间完全位于长区间内，因此只要求短区间被覆盖九成时能够保留。加入双向比例要求后没有输出，因为同样的一百碱基只占长区间的十分之一。

选择单向还是双向条件，取决于要判断“落在较大区域内”还是“两个区域范围接近”。窄峰与大型调控区相交时，双向严格门槛可能丢掉原本符合问题定义的记录。比较两个峰调用结果的边界一致性时，单向完全包含又可能显得过于宽松。门槛应由比较目的解释，而不是只凭数值看起来严格。

五、链方向改变的不只是行数

本例甲的三条匹配中，前两条是正链，第三条是负链。限制同链后，甲保留两条匹配，非重复覆盖长度从一百降到八十。乙保留一条负链匹配，覆盖从一百降到五十。丙只有反链注释，同链条件下覆盖变成零。

这种变化为检查实验设计提供了一个可见例子。如果查询表达的是具有方向性的转录事件，忽略链可能把反向重叠算作支持证据。如果查询只是无方向的开放染色质区间，机械套用链条件又可能排除本应保留的位置关系。不能因为输入文件恰好有第六列，就默认这一列具备适用于当前分析的方向含义。

还需要区分基因注释的链、比对记录的链和文库协议对应的转录方向。它们并不在所有测序方案中都一一相同。本次使用的是我们自己写入的方向标签，已经足够验证区间工具如何执行条件，但不能验证真实文库的链特异性。真实数据应先借助已知基因或专门质控核实方向约定。

六、相交、合并和转录本分块，不能互相替代

甲与乙没有共同碱基，但是默认合并操作把它们合成 100 到 300 的连续区间。这并不矛盾。相交要求共享位置，默认合并则允许边界相接。我们通过附加名称与方向汇总保留了来源，输出明确列出甲和乙，以及正负两个方向。

限制同链合并后，甲与乙保持分开。这个对照尤其适合检查“先合并全部注释，再算覆盖”的做法是否符合需求。合并后的一个大区间只表示空间上连续，不能据此断言其内部属于同一基因、同一峰事件或同一种调控状态。需要保留来源时，应该明确选择汇总字段，并检查汇总结果是否仍能回答原问题。

bedtools merge -i A.bed -c 4,6 -o distinct
bedtools merge -i A.bed -s -c 4,6 -o distinct
bedtools intersect -a transcript.bed -b intron.bed -u
bedtools intersect -a transcript.bed -b intron.bed -split -u

转录本例子进一步说明了“区间的内部结构”。我们写入一个从 100 到 400 的两块转录本，两个外显子分别位于 100 到 150、350 到 400，另放一个位于 250 到 260 的查询。按整个跨度判断时，两者相交。按分块判断时，没有相交，因为查询完全位于两块之间。

因此，想统计外显子覆盖，就必须让输入结构和命令共同表达外显子。仅把转录本总起止坐标当作普通三列区间，已经丢掉了内部断开的信息。反过来，如果研究的是整个基因座范围，分块后只检查外显子也可能改变问题。软件不会替我们选择这两个生物学对象。

七、排序是算法前提，不是装饰性参数

本次对已经有序的输入启用排序模式，并显式提供染色体长度文件，结果与普通求交逐字一致。随后故意把乙放在甲前面，破坏起点顺序，排序模式返回非零退出码并指出输入有序性问题。复现脚本将这次失败列为预期的诊断实验，不把它混进成功结果。

这说明批量处理时不能只看输出文件是否存在。任务失败时可能留下空文件，也可能留下部分内容。应同时保存命令、退出码与错误信息，让后续步骤在失败时停下。对于按染色体拆分的数据，还要确保不同文件遵循同一套染色体顺序，而不只是每条染色体内部按起点排列。

排序模式适合需要在有序区间上进行扫描的场景，但本次只有几个区间，没有测量峰值内存，也没有报告加速倍数。实际规模扩大后，应分别记录排序和求交的开销。如果为了省掉一次排序而反复读取混乱的大文件，或者每个小任务重复排序同一份注释，整个流程的代价可能主要来自外围准备工作。

八、把输出检查变成可复用习惯

一次区间操作完成后，可以先抽查三个位置。一个明确命中的查询，用来核对坐标和输出字段。一个边界相接但不相交的查询，用来核对端点约定。一个没有命中的查询，用来确认流程是否保留了必要的零值与分母。对于存在多重注释的项目，再加一个一对多匹配，通常能及早暴露错误计数方式。

BEDTools 的优势是能够把很多区间问题拆成清晰的小操作，并保留可以追查的文本结果。限制也来自这里，格式与命令都正确，不保证组合后的统计问题正确。尤其在峰注释、可变剪接和多轨道覆盖分析中，必须先定义记录身份、长度分母及内部结构，才能确定操作顺序。

下载本组三篇复现脚本与实测记录 后，在有 Docker 的 Linux 环境进入本工具目录，执行下面的命令。脚本使用固定镜像摘要，重建所有小文件，并记录成功操作及故意制造的排序失败。读取结果时，应核对具体案例的预期退出码，不应把所有非零码统一解释为环境安装失败。

python3 run.py output

参数定义可对照 intersect 官方文档、coverage 官方文档 与 merge 官方文档。本文的坐标、统计表与对照结果均由随文脚本构造和实际执行，适合学习与检查语义，不代表真实数据的分析结论。

当区间定义与统计口径确定后，若要将同样的流程用于大量样本或注释轨道，可携带文件规模、命令及复现记录，通过云生信生物分析服务器服务评估运行环境与批量任务安排。