拿到分箱结果后,很多人首先按完整度百分比排序,再决定保留哪些基因组。这个步骤有用,但百分比对应的对象必须讲清楚。一个基因组箱拿回了多少已知参考碱基,与标记基因模型在其中找到了多少预期信号,使用的是两种不同证据。我们把两条公开细菌参考完整切成片段,交给 MaxBin 2 实际分箱,再逐条核对回收结果。一次运行完整找回了大肠杆菌参考的全部四百多万碱基,软件报告的完整度仍为百分之九十八点一。这个看似矛盾的结果,适合用来理解分箱质量表究竟能支持什么结论。
从原始片段到质量百分比,中间经过了哪些判断
MaxBin 2 利用片段的序列组成和覆盖丰度进行分箱,标记基因又参与候选来源及质量评估。输入需要有统一的片段身份,丰度可以由原始读段回贴估计,也可以按要求提供现成的丰度文件。多个样本能提供同一来源随样本共同变化的线索。最终一个箱对应哪些片段,是由这些信息共同影响的。因此在解释结果时,要同时保留分组成员、丰度输入和使用的标记集合,不能只有汇总百分比。MaxBin 2 方法论文 · 作者发布页与说明
标记基因之所以能估计完整度,是因为某些保守基因在对应类群的完整基因组中具有可预测的存在情况。工具能识别到多少,就提供了关于缺失程度的线索。实际识别还经过编码序列预测、蛋白匹配和阈值筛选。片段边界落在一个基因中间,可能影响其完整识别,序列本身与模型差异较大,也可能使匹配不够充分。即使参考序列的每个碱基都在输入里,工具所看到的标记信号仍然可能随着片段化方式发生变化。

实跑经验 · 两条公开参考、两种切片长度和两套标记
这次运行 MaxBin 2.2.7,使用大肠杆菌 K-12 MG1655 和枯草芽孢杆菌 168 的公开参考,登录号分别为 NC_000913.3 和 NC_000964.3。读入的序列长度分别为四百六十四万一千六百五十二和四百二十一万五千六百零六碱基。脚本保存了规范化序列的校验值,方便复跑时确认取到的是同一版碱基序列。两条参考都完整进入了切片步骤,之间没有人为挖去任何区段。大肠杆菌参考记录 · 枯草芽孢杆菌参考记录
第一种输入从每条参考开头依次切出两万碱基的片段,最后不足两万的部分保留,总计四百四十四条。第二种每五千碱基切一次,总计一千七百七十三条。两种输入总长完全相同。为两种来源设计两列覆盖,大肠杆菌分别为三十五和八,枯草芽孢杆菌分别为七和三十八。相同来源的各片段使用相同的设计数值。这个实验从已知序列和可控丰度开始,没有从测序原始读段估计覆盖,也不涉及真实群落组装的误差。运行脚本、参考校验值与三轮结果
三轮分别使用两万切片加一百零七套标记、五千切片加同样标记,以及两万切片加四十套标记。最短参与片段均为一千碱基,线程数固定为二。实际通过 MaxBin 的正常命令入口进行基因预测、标记搜索和分箱。依赖程序包含 FragGeneScan、HMMER 和 Bowtie2,本例已经提供丰度,读段回贴步骤无需重复执行。安装后先检查入口帮助和这些依赖能否正常启动,再提交完整输入,能够较早发现路径或运行库缺失。
python3 108-maxbin2/run.py "$(command -v run_MaxBin.pl)" \
EC.fasta BS.fasta work108
# 其中一轮的等价命令,另选新结果前缀
run_MaxBin.pl -contig work108/chunks20k_107/assembly.fa \
-abund_list work108/chunks20k_107/abund_list.txt \
-out work108/recheck -min_contig_length 1000 -markerset 107 -thread 2
| 输入与标记设置 | 大肠杆菌回收碱基 | 大肠杆菌完整度 | 枯草芽孢杆菌回收碱基 | 枯草芽孢杆菌完整度 |
|---|---|---|---|---|
| 两万切片,107 标记 | 4,641,652 | 98.1% | 4,215,606 | 100.0% |
| 五千切片,107 标记 | 4,641,652 | 97.2% | 4,215,000 | 96.3% |
| 两万切片,40 标记 | 4,641,652 | 97.5% | 4,215,606 | 100.0% |
三轮都输出两个箱,成员按已知来源核对后,各箱均只含一个参考来源的片段。完整度直接取自软件汇总文件,参考碱基量由实际写出的成员求和。因为切片互不重叠,且事先知道每条片段从哪里来,这里求和就能对应回收参考长度。如果使用有重叠的片段或含多处可定位序列的真实组装,就应改用参考坐标并集并报告多重归属,避免把重复覆盖算成额外恢复的基因组。
分数下降了,究竟少了什么
在两万切片条件下,大肠杆菌的两百三十三条输入全部进入同一个箱,回收量与完整参考长度相等,报告完整度却是百分之九十八点一。改为五千切片后,大肠杆菌的九百二十九条输入仍全部回收,完整度下降到百分之九十七点二。这组对照说明,对于这条已知参考,标记指标会受到输入片段结构的影响。它不能被直接换算成“丢了百分之二点八的碱基”。若据此推算缺失的基因组总长,会和我们逐条核对的输入成员相冲突。
要进一步解释每一个百分点,应该打开标记列表,找出两轮哪些标记发生了变化,再把这些标记映射到原始参考和切片边界。若一个标记附近恰好出现断点,可以查看预测蛋白是否被分成两部分、各部分的匹配覆盖和得分如何。若它远离边界,则继续检查预测与模型识别。本次报告给出了碱基回收和标记分数的对照,还没有逐个定位所有变化标记的具体失败原因,因此不能把全部分差都归结为同一种边界效应。这样的追查顺序也适用于真实项目中两个软件质量指标不一致的情况。
枯草芽孢杆菌在五千切片条件下,确实少了六百零六碱基。这个数字可以直接从原始长度除以五千的余数得到。最后一段小于当前一千碱基门槛,进入了过短片段文件。它没有进入后续主要分箱,也没有出现在未分类结果文件里。该轮未分类文件为空,仍然不能说明所有原始输入都获得了分箱归属。把过短片段、未分配片段和所有写出箱的成员分别汇总,才能从输出端还原完整的输入账目。
这里还可做一次量级检查。六百零六个碱基占四百多万参考很小的一部分,而报告完整度从百分之百变为百分之九十六点三。两者数量级不同,就值得继续检查标记识别层。缺失一个小片段可能碰巧影响一个标记,切片也可能影响其他位置的标记预测。实际解释需要查具体标记的坐标,不能仅凭总长差就把两项变化一一对应。这个简单的量级核对能及时发现讨论中把两个不同分母混在一起的问题。
换一套标记,能不能直接比较哪个结果更好
第三轮只改成四十套标记,两个箱的成员和回收长度与第一轮完全一致,大肠杆菌的完整度却变为百分之九十七点五。这说明不同标记集合给同一个序列集合提供了不同的观察尺度。若一项筛选规则恰好卡在这两个数值之间,同一个箱就可能在两张表里得到不同的保留结论。项目中比较质量时应固定评估软件、模型或标记集合及版本,同时保存成员文件的校验值,让读者能区分分数变了还是箱内序列变了。
在一套标记很少的评价里,一个标记有无被识别对百分比影响会更大,这也是阅读分数时应保留分母的原因。百分比精确到一位小数,不等于对真实完整性的估计也具有同样精度。对临界箱,可以查看具体缺失标记、分类归属以及是否存在明显异常片段,再给出保留、继续修整或暂作候选的依据。若研究结论依赖一个关键代谢区域,则应直接核查该区域的序列及读段支持,不能只拿整体完整度替代功能位点的证据。
完整度还需要和混入情况一起看。两个不同来源的片段混在一起,可能增加某些标记的可检测数,也可能产生重复标记。只追求把完整度提高到某个数字,容易保留错误混合的箱。本例通过已知来源可以直接验证没有交叉混入。真实样本缺少这种标签时,应结合重复标记、各片段分类一致性、组成异常与覆盖模式来识别可疑成员,再检查移除后关键区域是否仍有足够支持。
丰度文件同样需要独立核验。先把每个文件首列与最终 FASTA 的片段标识做集合比较,再确认各列都是同一种深度定义。一个样本用平均覆盖,另一个样本用按总读数归一后的比例,会让本来可比较的变化方向混入单位差异。含有零值的片段应核对实际比对结果,尤其要确认数据缺失是否被误写成零。若多个样本的丰度文件依次写在列表里,应保存这一顺序与样本名称的对应关系,后续查看各箱丰度时才不会张冠李戴。
本例的两列覆盖是为了让来源容易区分而设计的数值,不能从中讨论培养条件、疾病状态或时间变化。替换为真实样本后,还需要检查样本之间是否独立、读段处理是否一致,以及某个批次的测序差异会不会恰好和生物学分组重合。先确保分箱输入有一致含义,再用得到的基因组箱讨论后续功能或群落变化,能减少把技术差异解释成生物差异的机会。
如何交付一份经得起复核的分箱质量表
每个箱至少保留一张成员清单,列出片段标识、长度、各样本丰度、主要分类依据和是否触及重点功能区。另留一张评估表,明确标记集合、软件版本、完整度和其他质量指标。对于有可靠参考的阳性对照,增加参考坐标覆盖及混入来源统计。对于真实未知对象,则保留所用判断的证据范围。这样即使之后更换评估器,原有成员集合仍然可以独立复查,比较结果也不至于只剩几列无法解释的数字。
流程结束时还应检查输入守恒。原始记录可以分别进入过短集合、未分配集合与各个箱,核对标识交集能够发现重复归属,核对集合并集能够发现被遗漏的记录。如果某个工具输出的未分类文件只覆盖通过过滤的输入,就单独读取过滤文件补足这部分。我们这次发现的六百零六碱基就属于这种情况。把这个小检查纳入脚本,比在下游找不到目标基因时再回头翻所有中间文件省事。
最后把结论写到实际证据所在的位置。本次已知参考的完整碱基回收,并未强制得到百分之百的标记完整度,切片方式和标记集合都影响了分数,短片段过滤还造成了一项可以定位的序列损失。真实分析中沿着同样顺序检查,便能逐步区分输入过滤、分箱归属和标记识别各自的问题。报告中每一项改进也有明确的复核对象,后续补数据或重新评估时更容易判断是否真正解决了原来的缺口。
若准备把这种逐箱核对扩展到大批样本,可用参考对照的实际输出量与云生信生物分析服务器团队讨论批量运行安排。
