同一份片段文件,文件行数没有变,MACS3 日志中的片段总量却可以不同,最终峰列表也会改变。原因不一定是程序重复读取了输入,而可能是指定格式后是否使用次数列。本次用一份五列人工片段文件,实际运行四组峰调用,把文件行数、次数贡献、条码筛选和输出坐标逐项核对。
这篇文章适合需要整理片段文件、按细胞集合调用峰或复查峰表的读者。实验使用完全虚构的坐标和条码,没有真实细胞、比对或测序质量信息。它展示格式与处理规则如何影响结果,不评价真实样本的检出率,也不把人工峰当作生物学发现。

先确认软件把一行理解成什么
本次安装并运行 MACS3 3.1.0,明确指定输入格式。五列文件依次包含染色体、起点、终点、条码和次数。使用 FRAG 格式时,程序会读取条码与次数,而简化的 BEDPE 模式主要按片段区间读取。对于这种有额外信息的输入,文件后缀相同并不保证解释方式相同。
还要注意,MACS3 所说的简化 BEDPE,与另一些区间工具使用的双区间格式并不是同一套字段定义。不能看到名字相同就不检查列。官方格式及选项要求见 MACS3 峰调用说明。本实验通过运行日志验证实际读取数量,而不是仅根据文档推算预期值。
这种检查应该发生在分析开始时。先抽看文件头、确认列数,再将日志中的总量与自己按定义计算的数量对照。若两者不同,应先解决计数口径问题,不要立即调整显著性阈值。否则可能在错误理解的输入上,逐步调出一份看起来合理的峰列表。
人工文件中安排了三个富集区
我们只使用一条人工染色体,并将有效基因组大小设为十万。所有片段长度均为一百五十碱基。这个数值是实验条件,不对应任何真实物种。真实分析需要采用适合参考版本和实验问题的有效基因组大小,不能为了得到更多峰照搬这里的小分母。
甲条码在一万附近有二十五个不同区间,每行次数为四。乙条码在三万附近也有二十五个不同区间,每行次数为四。另在一万五千附近为甲条码放置五个区间,每行次数为二十。此外,每个条码各有一百行分散背景,次数均为一。
| 输入部分 | 文件行数 | 每行次数 | 次数贡献合计 |
|---|---|---|---|
| 甲条码第一区域 | 25 | 4 | 100 |
| 甲条码第二区域 | 5 | 20 | 100 |
| 乙条码区域 | 25 | 4 | 100 |
| 两个条码的背景 | 200 | 1 | 200 |
| 合计 | 255 | 不同 | 500 |
这样,二百五十五行对应五百的次数贡献,甲条码单独占一百三十行、三百的次数贡献。所有数字都可以从生成脚本直接核对。两个甲区域虽然总次数相同,但一个由二十五行贡献,另一个只由五行贡献,适合观察限制次数后的差别。
次数列在本实验中只是人为设置的输入权重,不能自动解释为独立分子数。真实片段文件还需要结合其生成流程理解这一列,包括是否聚合了重复观测。把次数相加得到的量,和去重后的独立证据数量,不能仅因为都被称为片段就混为一谈。
四组调用保持了哪些条件
四组使用同一文件、相同有效基因组大小和零点零五的校正显著性门槛,均保存信号轨道和窄峰结果。第一组按 FRAG 读取全部条码,第二组只选甲条码,第三组把每行次数上限设为一,第四组按简化 BEDPE 读取并保留所有记录。
这里没有输入对照样本,没有运行宽峰分析或 HMMATAC,也没有进行插入位点移动。程序在片段模式下使用实际区间,日志中的片段长度为一百五十,并跳过了常规峰模型构建。不能把本次运行描述成已经验证了从真实双端比对开始的全部处理。
macs3 callpeak -t fragments.tsv -f FRAG \
-g 100000 -q 0.05 -n frag_all \
--outdir frag_all --threads 2 -B
macs3 callpeak -t fragments.tsv -f FRAG \
-g 100000 -q 0.05 -n frag_cap1 \
--outdir frag_cap1 --threads 2 -B --max-count 1
本次使用固定镜像,并记录软件版本、完整帮助信息与每组命令。两线程是本实验的运行预算,不能据此推算真实项目吞吐。日志还显示,这种单染色体小输入的峰扫描实际采用串行路径。因此,不能仅因为命令里写了线程参数,就声称已经测得多线程加速效果。
峰数量怎样随计数口径改变
按 FRAG 使用全部次数时,三个设计区域都进入峰表,峰顶处的堆叠高度均为一百。把每行次数限制为一后,日志中的总量变为二百五十五,只保留一万和三万附近的两个区域,两处峰顶堆叠高度都变为二十五。
| 运行条件 | 使用的次数贡献 | 峰数量 | 主要变化 |
|---|---|---|---|
| 全部 FRAG | 500 | 3 | 三个区域保留 |
| 仅甲条码 | 300 | 2 | 乙条码区域消失 |
| 每行次数上限为一 | 255 | 2 | 五行组成的区域消失 |
| 简化 BEDPE,保留所有行 | 255 | 2 | 与上限为一的结果一致 |
在本次输入上,简化 BEDPE 与次数上限为一的峰区间和数值一致,只是结果名称不同。这是因为两者都让每行贡献一次,且实验区间没有额外的重复行。不能推广成这两种处理在所有真实文件中都等价,尤其当输入包含其他重复结构或格式差异时。
五行组成的区域在限制次数后没有进入最终峰表。我们没有将原因单独归为某一个未经隔离验证的门槛,因为峰结果同时受信号、背景、显著区间长度和合并规则影响。可以确定的是,修改次数处理改变了输入贡献,继而改变了完整调用输出。若想判断具体由哪条规则淘汰,还需进一步检查相应中间轨道和条件。
这个对照最直接的经验是,不要在看到少一个峰后,马上认定软件漏检。先比较程序真正使用了多少记录贡献,以及改变的是输入、背景还是过滤规则。只统计文件行数,无法解释本实验中五百和二百五十五之间的差额。
条码筛选要看筛选之后的日志
只选择甲条码时,日志先报告读入的总贡献为五百,随后明确写出提取了三百。若只搜索第一处总片段数量,就会误认为筛选没有生效。完整日志还显示,此时用于扫描的区间位置数量是一百三十,与甲条码的文件行数对应。
最终表格中只剩甲条码的两个区域,乙条码对应的区域消失。这验证了本次条码列表确实参与了筛选。但筛选成功不等于细胞身份正确,本实验的条码名称都是人工指定,没有进行细胞识别、双细胞检测或分群。
实际使用细胞分群结果时,应核对条码是否与片段文件完全一致,包括可能附加的样本前缀或后缀。多个样本之间若复用了相同条码字符串,还需要确认上游如何区分样本来源。只看条码列表长度,不能保证选到的是预期的那一组细胞。
macs3 callpeak -t fragments.tsv -f FRAG \
--barcodes barcodes_A.txt \
-g 100000 -q 0.05 -n frag_CellA \
--outdir frag_CellA --threads 2 -B
本次日志还提示 FRAG 模式会保留次数贡献。若要测试该模式中的次数上限,应使用对应选项,而不是假设其他格式的重复处理参数具有同样效果。是否限制次数应由文件定义和分析目标决定,不能因为本例中峰更少,就把某个设置统一称作更严格或更正确。
同一处峰的高度相同,校正值仍可能不同
甲条码第一个区域在全部输入和条码子集中都保留,峰顶堆叠高度同为一百,其负对数原始 P 值也相同。但负对数校正值从约一百二十二点三五三变为一百二十点五二七。候选位置未变,不代表所有统计列必须逐位相同。
这提醒我们,校正结果与当前分析范围和统计背景有关。不能把不同细胞集合分别调用峰后得到的校正值,直接当作同一标尺上的细胞活性强弱。即使某一局部信号没有变化,结果提取的整体情境也可能改变。
本实验没有进一步运行差异可及性分析,因此不能根据两个峰表的有无或校正值大小,得出甲乙两类细胞之间存在统计学差异的结论。需要比较条件时,应明确样本或生物学重复、共同检测区域及计数模型,不能把独立调用的两个门槛简单相减。
同理,极小的校正值只描述当前模型下的信号与背景关系,不代表这一区域已经验证了某种调控功能。人工数据中的高度富集是由我们设计出来的,恰好提醒读者不要把数值极端与生物学意义强画等号。
输出坐标需要用同一套规则解释
本次同一条峰在制表符结果中的起点是一万零四,在窄峰文件中的起点是一万零三,而终点均为一万零二百一十六。前者采用从一开始的闭区间表示,后者采用从零开始的左闭右开表示。两种写法对应同一段区域,不是程序把起点移动了一碱基。
| 核对内容 | 本次输出 | 如何核对 |
|---|---|---|
| 制表符峰区间 | 10004 到 10216 | 长度 213 |
| 窄峰区间 | 10003 到 10216 | 终点减起点为 213 |
| 窄峰内峰顶偏移 | 108 | 加到零起点得到 10111 |
| 峰顶文件区间 | 10111 到 10112 | 表示一个碱基位置 |
坐标转换最适合用这样的小例子验证。把两个起点直接比较会看到差一,但同时检查区间长度和峰顶位置,就能确认是否只是坐标体系不同。若后续要与基因注释求交,应先统一参考版本和坐标定义,再解释重叠数量。
文件名后缀也容易误导。本次带表格后缀的结果实际是含注释行的制表符文本,不能因为后缀看起来像电子表格,就按二进制工作簿读取。读取时应处理开头说明行,并保留其中的版本和参数信息,不宜只留下数据区后丢失分析条件。
窄峰文件中的相关显著性列是负十进制对数,不能直接当作原始概率。峰顶偏移也不是基因组绝对坐标。自行合并结果之前,先为每一列写清单位与定义,能够避免后续坐标注释和排序中出现看似合理、实际错位的结果。
怎样复现与检查自己的流程
完整脚本和实际输出见 本组实跑资料。驱动程序生成片段与条码文件,运行四组峰调用,保存日志、峰表、峰顶以及关键轨道区间。进入 MACS3 对应目录后,在有 Docker 的 Linux 环境中执行下面的命令。
python3 run.py run161
复现时先核对三个数字,全部文件二百五十五行、次数合计五百、甲条码次数合计三百。再检查峰数量分别为三、二、二、二,最后核对坐标及峰顶偏移。这个顺序能把格式问题与峰调用问题分开,避免一开始只盯着最后一个峰的显著性。
若替换成自己的数据,应先从少量可人工核对的记录验证解析方式,再扩大运行范围。尤其是格式转换、条码合并和次数聚合之后,要重新检查总量。一次转换即使没有丢行,也可能改变每行的含义,因而不能仅以行数一致作为转换正确的证据。
软件当前功能与安装要求可查 MACS3 项目文档。本次固定版本的程序已实际运行,但没有将在线文档中提到的其他子命令和性能改进当作本篇的实测结果。
优点、局限与选用
MACS3 的实用优势是可以明确处理不同片段格式,并保存峰调用所依据的轨道与日志。对于需要按指定条码集合生成候选区域的流程,这种可追溯输出有助于核对输入与结果的关系。本次四组对照之所以能够解释,是因为没有只保存最终峰列表。
它的局限是不会自动判断输入中的次数是否代表独立观测,也不会替代样本质量、参考选择和实验设计。调用出的富集区并不自动等于真实结合位点或因果调控区域。若输入背景和计数口径不合适,再小的显著性数值也不能消除这些问题。
最值得带走的习惯,是先说明一行是什么、次数如何使用、条码筛选发生在哪里,再讨论峰是否可信。把这条证据链保留下来,后续出现峰数量变化时,就能从实际处理步骤查原因,而不是反复尝试阈值直到得到预期的名单。
当已经核对的片段与峰调用流程需要扩展到更多样本或细胞集合时,可按数据规模了解 云生信生物信息学分析服务器服务。
