📚 课程目录第3章 Kraken2原理简介

Kraken2 原理简介：基于 k-mer 精确匹配的分类

好，咱们来聊聊 Kraken2 的原理。说实话，我第一次接触它的时候，心里想的是：“这不就是字符串匹配嘛，有啥高级的？” 但后来我发现，它把简单粗暴做到了极致，反而在某些场景下比那些花里胡哨的模型还管用。

Kraken2 的核心挺简单的。说白了就是基于 k-mer 的精确匹配。你给一条读长进来，它切成固定长度的小片段，也就是 k-mer，然后去数据库里一个个比对。全部匹配上了，就知道这条读长属于哪个物种。

核心思想：精确匹配 + LCA 打分

你想想看，一条 150 bp 的读长，切成 31 bp 的 k-mer（Kraken2 默认是 k=31），大概能切出 120 个片段。每个片段在数据库里查一遍，看它出现在哪些基因组里。最后再拿这些结果做一个 LCA 投票。

我个人习惯把这个过程叫作「指纹比对」。每个基因组都有自己独特的 k-mer 指纹。Kraken2 干的活儿，就是把未知读长的指纹跟已知的指纹库做精确匹配。匹配的越多，就越确定。

精确匹配：不考虑错配，完全相同的 k-mer 才算 match
LCA 投票：同一个 k-mer 出现在多个物种里，就往上一级祖先投票
最终标记：统计所有 k-mer 的投票结果，得分最高的那个物种就是分类标签

关键区别：Kraken2 不做序列比对，也不做概率模型。它完全依赖 k-mer 的精确出现与否。这意味着它非常快，但也意味着它对数据库的完整性要求极高。数据库里没有的物种，它绝对认不出来。

k-mer 索引：其实是个哈希表

为什么 Kraken2 这么快？因为它把数据库里所有参考基因组的 k-mer 都做成了哈希表。每次查询就是 O(1) 的时间复杂度。我记得有次在 96 核服务器上跑全宏基因组，Kraken2 花了不到 5 分钟就处理完 5000 万条读长。同时跑的 MetaPhlAn 还在那儿吭哧吭哧地比对标记基因呢。

具体的索引结构是这样的：

哈希表键（Key）: 每一个 k-mer（31 bp 的序列）
哈希表值（Value）: 该 k-mer 所属的 LCA 分类号（NCBI taxonomy ID）

索引构建时的处理流程：
1. 把参考基因组切成所有可能的 k-mer
2. 每个 k-mer 找到它所在的物种分类
3. 如果某个 k-mer 出现在多个物种里，取这些物种的 LCA
4. 把 (k-mer, taxonomy_id) 对写入哈希表

经验之谈：我刚开始用 Kraken2 时，一直以为它的数据库越大越好。后来发现，数据库太大反而会引入大量冗余 k-mer，导致分类结果变模糊。比如人肠道样本，如果数据库里加了大量环境菌的 k-mer，很多读长就被分配到更高级的分类了，因为 k-mer 在多个物种里都出现了。

数据库构建：从基因组到 k-mer 映射

Kraken2 的数据库是怎么来的？不是你想的那样直接拿 FASTA 文件就行了。你需要先下载参考基因组，然后用 kraken2-build 工具来建库。

咱们直接看命令吧：

# 下载标准数据库（细菌、古菌、病毒、人类）
kraken2-build --standard --db $DBNAME

# 或者你只想要细菌的数据库
kraken2-build --download-library bacteria --db $DBNAME

# 手动添加某个参考基因组
kraken2-build --add-to-library custom_genome.fna --db $DBNAME

# 最后一步，构建哈希索引
kraken2-build --build --db $DBNAME

整个建库过程，说白了就是把参考基因组里所有长度为 k 的子串都提取出来。举个例子，一个只有三个序列的简单库：

参考序列：ATCGATCGATCG

k=4 时提取的 k-mer 列表：

ATCG
TCGA
CGAT
GATC
ATCG
TCGA
CGAT

注意这里 ATCG 出现了两次。Kraken2 在构建索引时会做去重，所以哈希表里只会保留一份 (ATCG, taxonomy_id) 条目。

分类过程：读长进来，k-mer 投票

嗯，这里要注意。分类并不是跑完所有的 k-mer 就直接出结果。Kraken2 里面有一套置信度机制。简单来说就是：

读长上所有 k-mer 先做 LCA 投票。每个 k-mer 得到一个分类号，所有分类号按出现次数排序。
计算覆盖度。匹配上的 k-mer 数量占读长总 k-mer 数的比例。如果一条 150 bp 读长有 120 个 k-mer，但只有 30 个匹配上了，那覆盖度就是 25%。
决定丢弃线。Kraken2 默认有一个最低覆盖度阈值（--minimum-hit-groups），低于这个值的读长会被标记为 unclassified。

Kraken2 分类结果的典型格式
字段	说明	示例
C/U	分类状态 (Classified/Unclassified)	C
读到 ID	原始 FASTQ 序列 ID	@E00526:123:ABCDE
分类号	NCBI taxonomy ID	562 (Escherichia coli)
读长长度	bp 数	150
LCA 映射	每个 k-mer 对应的分类号列表	562:18 560:2 561:5

我记得在项目里遇到过一个问题：Kraken2 分类结果里有很多「Unclassified」，当时以为是数据库太旧了。后来发现是读长质量太差，大量 k-mer 因为测序错误而匹配不上。你想想看，k=31 的时候，只要有一个碱基错了，这个 k-mer 就算白搭。

局限性与优化技巧

Kraken2 虽然快，但并不是万能的。我总结了几条避坑指南：

k 值的选择：默认 31 bp 是平衡了灵敏度和特异度的结果。k 越大，特异性越高，但有些短的读长根本切不出几个完整的 k-mer。k 越小（比如 19），灵敏度高但容易误分类。
数据库大小：标准数据库大概 8-10 GB。如果你只关注呼吸道病原体，完全没必要建全库，用 --quick 模式配合自定义库更省内存。
内存占用：哈希表会整体载入内存。我曾经在 64 GB 的服务器上跑标准库，结果内存直接拉爆。后来改成 --memory-mapping 才勉强能跑。
分类粒度：Kraken2 在种水平上表现不错，但到了菌株水平就基本没戏了。因为菌株之间共享的 k-mer 太多了，LCA 会直接拉回到种水平。

避坑指南：我曾经在评估 Kraken2 和 MetaPhlAn 的一致性时，发现两个工具在低丰度物种上分歧很大。后来定位到问题是：Kraken2 对低丰度物种的 k-mer 匹配数量太少，置信度不足，直接标记成 Unclassified 了。这时候如果你用默认参数输出分类报告，那些低丰度物种就彻底丢了。后来我调整了 --minimum-hit-groups 参数，才把结果找回来。

这张图概括了 Kraken2 的完整分类链路。你看，从读入序列到输出分类标签，中间没有冗余的比对步骤。这也是它能在几分钟内搞定全宏基因组样本的原因。