📚 课程目录第3章 Kraken2工作原理

Kraken2 工作原理：k-mer 匹配与 LCA 算法

Kraken2 的核心思路其实很直接——把测序 reads 切成短小的 k-mer，然后在预先建好的数据库里查这些 k-mer 属于哪个物种。说得再直白点：就像用一堆小拼图去比对一张巨大的物种地图，每块小拼图都指向地图上的一个节点，最后把节点投票最多的那条分支作为最终分类。

从 read 到 k-mer：第一步切分

我习惯把 k-mer 长度设在 31 bp，这是 Kraken2 的默认值，也是很多项目里平衡速度与精度的折中点。你想想看：k-mer 越短，匹配到的数据库条目就越多，但特异性下降；k-mer 越长，特异性高，但能匹配上的 k-mer 数量骤减，灵敏度可能掉得很厉害。

k-mer 长度	匹配数量	特异性	典型场景
21 bp	多	低	高错误率数据（如 Nanopore）
31 bp	中	中	Illumina 标准推荐
51 bp	少	高	高覆盖度、需要种级分辨率

记得有一次我在做肠道宏基因组分析，默认 31 bp 跑下来，大批 reads 被打上了“Bacteria（细菌域）”的标签。我后来一查，原来是样本中很多菌株太新，数据库里只有属以上水平的 k-mer 覆盖。把 k-mer 长度改到 21 后，部分 reads 成功落到了科级，反而比域级更有生物学意义。嗯，这里要注意——参数不是越大越好，得看你的参考数据库覆盖深度。

k-mer 匹配：查表，快就一个字

Kraken2 把数据库建成一个紧凑的哈希表，key 是 k-mer 的 minimizer（一种简化指纹），value 是这个 minimizer 对应的 LCA 分类 ID。匹配时只需要对每条 read 滑动窗口提取每个 k-mer，计算 minimizer，再去哈希表里取分类 ID。整个过程 O(1) 查表，速度远快于 BLAST 式的比对。

我个人很喜欢这个设计——它不需要做序列比对，只有哈希命中或未命中。代价是，如果数据库里没有该 k-mer，这条信息就丢了。这也是为什么 Kraken2 对数据库完整性要求高，建库时如果漏掉了某些物种的 k-mer，结果里就会出现大量 unclassified。所以我的建议是：建库前一定要检查物种代表性，别只塞公共数据。

# 典型的 Kraken2 建库命令（针对标准五数据库）
kraken2-build --standard --threads 8 --db $DB_DIR

# 对测序 reads 进行分类（输出为 Kraken 格式）
kraken2 --db $DB_DIR --threads 8 --output sample.kraken \
        --report sample.report sample_R1.fastq

输出报告里每一行都是一个 read 对应的分类 ID 和匹配到的 k-mer 数量。如果匹配到的 k-mer 太少（默认阈值1），会让 Kraken2 直接丢掉这条 read。

LCA 算法：把碎片拼成完整图谱

每条 read 会对应几十到上百个 k-mer 匹配结果，每个结果都指向一个分类节点。Kraken2 的 LCA 算法不是简单取多数，而是用一种“路径累积分数”的方法：

对每个 k-mer 匹配到的分类节点，把它到根节点路径上的所有节点都加 1 分；
统计所有路径得分，选出得分最高的节点；
如果最高分节点不止一个，取它们的最近共同祖先（LCA）。

为什么要这么做？因为一个 read 里如果大部分 k-mer 都指向大肠杆菌（E. coli），那整个 read 很可能就是大肠杆菌的片段；但如果有少量 k-mer 指向了沙门氏菌，它们也会给大肠杆菌的祖先节点（比如肠杆菌科）贡献分数，使得肠杆菌科的总分有时会超过大肠杆菌本身。Kraken2 会选择得分最高的那个节点，这就是最终分类。

核心一句：Kraken2 的 LCA 不是“所有匹配节点的最近祖先”，而是“得分最高的节点或它的最近共同祖先”。这个细节很多人容易搞混。

我曾经踩过一个坑：在分析一个混合样本时，发现 read 明明有大肠杆菌和志贺氏菌的 k-mer 同时命中，结果却被归到了“Shigella flexneri”。后来查源码才明白，志贺氏菌的 k-mer 数量更多、路径得分更高，而大肠杆菌的 k-mer 反而成了少数派。这类近缘物种的区分，Kraken2 会偏向 k-mer 数更多的那个。如果两个分类的 k-mer 数量几乎相等，LCA 才会回退到它们的共同祖先。

可视化流程：一图胜千言

下面这张图概括了从 read 到最终分类的全链条。我把每个阶段的输入输出都标出来了，你照着看就能把原理串起来。

参数调优与避坑

我的个人习惯：先拿 1 万条 reads 跑一次默认参数，观察 unclassified 比例和分类深度。如果 unclassified > 20%，降低 minimizer 长度（–minimizer-len）或 k-mer 长度；如果整体分类深度都停留在科级，说明数据库对该样本覆盖不足，考虑增加定制库。

我曾经犯过的错误：在一次病毒宏基因组分析中，我用默认的 31 bp k-mer 跑完，结果 90% 以上的 reads 都成了 unclassified。后来发现病毒基因组 GC 含量极低，31 bp k-mer 在那段区域几乎找不到完全匹配。调低到 25 bp 后，unclassified 降到了 30%。这类低复杂度区域多的样本（病毒、高 AT 生物）一定要降低 k-mer 长度。

还有一个容易被忽略的点：Kraken2 的 LCA 算法对重复区域敏感。如果某个保守基因在多个物种中完全相同，那么对应 k-mer 的 LCA 就会上浮到这些物种的共同祖先。这其实是合理的——同一个 k-mer 无法区分大肠杆菌和志贺氏菌，LCA 给出肠杆菌科是最保守的正确选择。但你要是希望得到种级分辨率，就得指望那些种间有差异的 k-mer 数量能压过保守区。

小结

Kraken2 的思路其实很工程化：用哈希表实现最快匹配，用路径累积投票代替复杂比对。k-mer 长度决定灵敏度和特异性的天平，LCA 算法则负责把碎片化的信号合成一个可信的分类结论。理解了这两块，你就掌握了 Kraken2 的核心。