2. 第2章

好，这一章咱们要把 Kraken2 的“脑子”拆开看看。说白了，它到底是怎么从一堆乱七八糟的测序片段里认出来哪个细菌、哪个病毒的？我自己开始用这个工具时也是一头雾水，直到把它的四个核心机制摸透——k-mer 匹配、LCA 算法、分类树结构和置信度评分，才算真正能调好参数。

嗯，这一章不求快，咱们一个一个啃。

k-mer 匹配：打碎序列找指纹

Kraken2 的第一步，就是把每条 reads（读长）切成固定长度的小片段。这个片段就是 k-mer。比如你设 k=31，一条 150 bp 的 reads 会被切成 120 个 31-mers（150-31+1）。

我常跟团队说：k-mer 就是序列的“指纹”。每个物种都有独特的 k-mer 组合。Kraken2 建立了一个巨大的 k-mer 到分类单元的映射数据库，就像字典一样。

匹配过程很简单：

对每个 k-mer，在数据库中精确查找
找到就记录该 k-mer 所属的物种（或更高分类层级）
找不到就记为未分类

这里有个关键参数——k-mer 长度。我做过测试，k 越小，匹配越容易但特异性差；k 越大，特异性强但敏感度下降。下表是典型选择：

k-mer 长度	优点	缺点	常见场景
20–25	敏感，对短 reads 友好	多匹配到共有序列，LCA 偏向高分类层级	短读长（< 75 bp）、低丰度物种
31（默认）	平衡灵敏度和特异性	对短 reads 不友好（150 bp 只给 120 个 k-mer）	常规宏基因组（100–150 bp）
35–50	高特异性，区分亲缘物种	敏感度下降，数据库更大	高深度、长读长（> 250 bp）

曾经有个项目，我用默认 k=31 跑口腔样本，结果链球菌和口腔链球菌被混在一起。后来把 k 提到 35，才把它们分开——代价是内存需求涨了 30%。这就是权衡。

LCA 算法：解决“k-mer 打架”问题

一个 k-mer 可能同时匹配多个近缘物种。比如某个 31-mer 在数据库里既匹配大肠杆菌又匹配志贺氏菌，怎么办？Kraken2 不是选一个，而是用 最低公共祖先（LCA）算法搞定。

LCA 的规则很简单：

取出该 k-mer 对应的所有分类单元（如属、种、株）
沿分类树向上回溯，找到这些单元最近的那个共同父节点
把这个父节点作为该 k-mer 的“投票”结果

举个例子：k-mer 12345 同时匹配大肠杆菌（Escherichia coli）和志贺氏菌（Shigella dysenteriae），它们的 LCA 是 肠杆菌科（Enterobacteriaceae） 这个层级。所以这个 k-mer 的投票指向肠杆菌科。

我的一点经验： LCA 是 Kraken2 比旧版本准确的关键。但它也有副作用——当数据库里近缘物种太多时，投票结果容易“向上漂移”，很多 reads 都被归到科甚至纲。这时候需要结合置信度评分来校正。

分类树结构：Kraken2 的“家谱”

Kraken2 用 NCBI 的分类树，但做了简化。它只保留有参考序列的节点，不保留空节点。这是一个树状结构，每个节点代表一个分类单元（域、门、纲、目、科、属、种）。

匹配过程汇总后，对每条 reads，Kraken2 统计所有 k-mer 投票给哪个分类节点，然后把 read 分配给得票最多的那个节点——前提是得票数超过阈值。

我整理了一个分类树片段示例（从细菌到种）：

root (1)
  Bacteria (2)
    Proteobacteria (1224)
      Gammaproteobacteria (1236)
        Enterobacterales (91347)
          Enterobacteriaceae (543)
            Escherichia (561)
              Escherichia coli (562)

注意：树越深，分类粒度越细。Kraken2 默认输出从“种”到“界”所有层级的占比，但你也可以用 --report-minimizer-data 控制报告里只显示哪些层级。

置信度评分：给结果上个“保险”

这是 Kraken2 最容易被忽略但极其重要的参数。我给它个比喻：信心指数。

每条 read 被分配到某个分类节点后，Kraken2 会计算一个置信度分数：

置信度 = 1 - (第二多票节点得票数 / 最多票节点得票数)

如果第二多的票数和最多的很接近，置信度就低。满足该 read 置信度大于等于你设定的 --confidence 阈值（默认 0），才会被纳入最终丰度计算。

注意踩坑： 我曾经把置信度设为 0.7，结果 half 的 reads 被丢弃，原本差异显著的样本变得没有区别。反过来，设太低又会被近缘物种干扰。建议先跑一遍统计置信度分布，再选一个让 超过 80% 的 reads 保留的阈值。

下表给你一个参考（基于我测试的一个土壤宏基因组）：

置信度阈值	保留 reads 百分比	分类精度	应用场景
0.0（默认）	100%	低，多 false positive	快速筛查
0.1	~90%	中等	一般分析
0.3	~75%	高	关键物种鉴定
0.5	~55%	极高	只信任最佳匹配

嗯，总结一下这四者的关系：先通过 k-mer 匹配得到一堆候选节点，再用 LCA 合并近缘冲突，接着用分类树定位节点层级，最后用置信度过滤掉不确定的结果。整个过程就像个流水线——

下面这张图帮你理清流程：

最后再说一句：这四个原理你一定要动手验证。我自己的习惯是拿一小段已知物种的 reads 跑一下 Kraken2，打开 --classified-out 选项看每条 read 的归属，再对照 LCA 和置信度影响，才能真正理解它为什么这样工作。