software

同一种菌,Kraken 报 2.69%,Bracken 为什么报 5.53%?

作者:发布于 2026-11-14约 8 分钟阅读
Bracken 官方示例复算。大肠杆菌的 Kraken 分支占比为 2.69%,Bracken 物种级占比约 5.53%,按相同的全输入分母重算为 5.44%。

宏基因组分析结束,做物种组成图的人从 Kraken 报告里抄了大肠杆菌的百分比,整理补充表的人用的却是 Bracken 丰度列。两个数对不上。重新跑分类之前,可以先用项目作者提供的一份报告核账。同一个物种、同一批输入,沿着分类树往下看,再弄清两个百分比分别除以了什么,就能解释这道差额。

本文复算使用 Bracken 作者公开的 Kraken 示例报告、配套的 75 nt 分布文件和物种层级的估算程序,最低计数参数设为 10。示例是较早的 Kraken 报告,下面读列的方法同样适用于 Kraken 2 默认的六列报告。所有数字来自这次复算。读者可以用自己的文件替换文中的输入。作者发布的示例文件

Kraken 的同一行,为什么有两个计数

默认报告的前三列依次记录该分类分支占全部输入片段的百分比、落在该分支及其所有后代的片段数、直接分到当前节点的片段数。接下来的三列是分类层级、分类编号和名称。片段是本次分类的单位,成对模式下一对读段合计一个片段。先确认拿到的是样本的标准报告。逐条读段的分类输出、带额外 minimizer 列的报告以及 MPA 风格报告的列布局都不同。Kraken 2 样本报告说明

示例中的行 该分支及后代 直接落在本节点 第一列占比
未分类 425,595 425,595 1.60%
根节点 26,241,409 0 98.40%
埃希菌属,编号 561 717,629 0 2.69%
大肠杆菌,编号 562 717,629 449,674 2.69%

这里有两个很容易漏掉的事实。第一,物种行上的 717,629 已经包含更细层级的片段,直接停在种节点的只有 449,674。其余 267,955 个落到了这个物种以下。第二,属行也显示 717,629,是因为子树向上汇总。同一批片段会同时出现在属、种以及更高层的累计数里。把一张完整报告的第二列直接求和,等于沿分类树重复计数。要复算总输入,用根分支加上未分类即可,本例是 26,241,409 加 425,595,合计 26,667,004。物种分支的 717,629 除以这个总数,就得到 2.69%。

如果只想确认某个物种,按分类编号提取比按名字模糊搜索可靠。下面的命令读取制表符分隔的六列报告,并保留整行供核对。先看文件里是否确实有六列,再把编号换成目标物种。名称相似的菌株和上级属就不会混进来。

awk -F '\t' 'NF == 6 && $5 == 562 {print}' sample_test.report

Bracken 改变了分子,也改变了分母

近缘物种共享一部分序列。分类器可以把一条片段放到大家共有的祖先节点,单凭它暂时分不清来自哪个种。Bracken 使用与分类数据库及读长相匹配的分布文件,估计上层节点中的一部分片段该怎样分给下属物种。输出里的 added_reads 是模型重分配的估计量。它不会产生新的测序片段,也没有把每条模糊片段逐一证实到了某个物种。

这次重算,大肠杆菌一行的 kraken_assigned_reads 是 717,629,与种分支累计数一致。再加 added_reads 的 732,184,得到 new_est_reads 的 1,449,813。即便其中的属节点直接计数为零,更高层共同祖先的片段仍可能参与分配。因此,拿“属行第三列”当作该种所有新增计数的上限,也会对不上账。Bracken 的估算方法与参数

现在可以对照两个分母了。若把 1,449,813 除以全部输入片段数,结果约为 5.44%。Bracken 表里的 fraction_total_reads 则是 0.05525,换成百分比约为 5.53%。作者的实现用本次在目标层级纳入估算的总量归一化。本例物种层级估计总量约为 2,624 万,未分类片段不在这个分母里。按物种行导出的整数计数或舍入后的小数重加,末位可能略有出入。Bracken 丰度列的计算实现

Bracken 官方示例复算。大肠杆菌的 Kraken 分支占比为 2.69%,Bracken 物种级占比约 5.53%,按相同的全输入分母重算为 5.44%。
同一份公开数据,先核对分子的重分配,再比较两个分母。图中百分数四舍五入到小数点后两位。

物种柱子翻了一倍,先做一个分母实验

假设另有两份教学样本,各投入十万个片段。目标种在两次估算中的片段数都为一万。甲有八万个片段进入物种层级的估算,乙只有四万个,此时目标种在物种组成图上的占比分别是百分之十二点五和百分之二十五。如果两张图只展示经过归一化的彩色柱子,乙看起来恰好翻倍。用全部输入再算,两者却同为百分之十。分子从头到尾没有增加,差异发生在能进入这张组成图的片段集合。

碰到自己的表有类似变化,先按样本列出总输入、未分类片段、参与种级估算的总量和目标种的估计数。接着挑未分类率高的一份检查,读段质量、参考库是否覆盖该环境、上游是否做过不同的宿主过滤。若只重新运行了物种丰度估算,应保留第一次的分类报告,这样可以把分类时发生的变化与后续重新分配的变化分开追查。

实跑经验|先用五条序列看懂最低共同祖先,再复算公开报告

在 Ubuntu 24.04 中实际运行 Kraken 2.1.3,用两条人工参考序列建了一个小库。两条参考各有一段各自独有的序列,末尾另有三百五十个碱基完全相同。从甲的独有区取两条读段,从乙的独有区取一条,从共有区取一条,再加一条独立生成的随机读段,每条七十五个碱基。逐条输出依次落到甲种、甲种、乙种、共同的属、未分类。这个五条输入的实验直接说明,属这一支一共收了四条,其中只有一条停留在属上。它们在同一行显示为两个不同的计数。

随后在相同环境实际编译了 Bracken,从项目固定提交中取官方数据库的七十五碱基分布文件和配套报告,使用物种层级、阈值十调用完整命令。程序脚本标示为 3.0.2,大肠杆菌一行给出原始计数 717,629,重新分配增加 732,184,合计 1,449,813,输出比例 0.05525,与上文算例一致。第一次调用遇到一个能复现的安装问题,该提交的外层脚本要调用名为 python 的解释器,而这份 Ubuntu 只有 python3。把分析环境内的 python 命令指向同一版本的 Python 3 后即可完成,无须改动数据库。实际排错时,先看最先出错的解释器或文件名,再去怀疑生物学分母。本次使用的 Bracken 命令入口及解释器调用

换成自己的样本,先核对这三个条件

第一,看分类和估算是否用同一套数据库。目录名相似不足以确认,最好保留建库日期、物种范围和实际调用路径。第二,看分布文件的读长是否和处理后的数据相符。下面的 75 来自作者配套示例。如果你的文库读长不同,要换成与之匹配的设置和库文件。第三,检查比较的是同一层级。用种级估计画一组、用属级分支占比画另一组,两边列名即使都写“相对丰度”,数值的含义也不一致。

如果已准备好相匹配的已建库目录,可用官方命令行入口。先把示例路径换成自己的数据库目录,这里演示种级估算。参数 -t 10 影响低计数分类单元是否进入估算,它不是物种存在的生物学判定线。做多样本比较时,把这项设置与数据库、读长、分析层级一并写入方法记录。

DB=/path/to/matched-kraken-db
bracken -d "$DB" -i sample.kreport -o sample.bracken \
  -r 75 -l S -t 10

最后留两列可以省去很多解释。物种组成图用 Bracken 的物种级相对丰度,并在图注中写明归一化口径。样本质控表另留 Kraken 的未分类比例。尤其在两个样本未分类比例相差很大时,明显变高的物种百分比可能部分源自分母变化。若研究问题是某种菌的绝对数量有没有增加,还需要总菌量、外源标准或独立定量实验等额外证据。把目前能回答的问题写准确,比再给柱状图加一位小数更重要。

把这一份样本的数字核清,再固定数据库和参数去跑整批队列,后面的表才可以直接比较。如果团队需要长期保留这套数据库、反复处理新增样本,可了解云生信的生物信息分析服务器。

把知识用于当前任务

下一步可以直接开始分析

按研究目标查找在线工具、课程、计算资源或完整分析服务。

按目标找方案仍有疑问,联系客服 ↗