宏基因组分类基础：分类学层级与 k‑mer 原理
做宏基因组分析，你最先要搞明白的就是：我们到底怎么把一堆乱七八糟的 DNA 片段，贴上“这是大肠杆菌”、“那是拟杆菌”的标签？说白了，核心就两件事——分类学层级和 k‑mer 原理。Kraken2 能又快又准地干活，靠的就是这两套东西组合在一起。
我个人刚开始学的时候，总觉得分类层级太简单，不就是界门纲目科属种嘛。但真跑到项目里，你才会发现，每一步都可能踩坑。比如有一次我拿 Kraken2 跑一个肠道样本，结果看到大量“未分类”的序列，后来才发现是 k‑mer 长度设得不对，导致很多 k‑mer 没匹配上数据库。嗯，从那以后我就特别重视这两个基础概念了。
1. 分类学层级——从域到株
分类学层级，就是生物学家给所有物种排的家谱。从上到下一般是：域 → 界 → 门 → 纲 → 目 → 科 → 属 → 种 → 株。Kraken2 的数据库把每个物种的基因组都切成 k‑mer，然后给每个 k‑mer 打上它所属的分类层级标签。你想想看，一条 reads 里有几十个 k‑mer，每个 k‑mer 都有各自的分类标签，Kraken2 再通过路径一致性算法，把这个 reads 归到最可能的分类节点上。
关键点：Kraken2 的输出既可以是“种”级别，也可以是“属”级别，甚至“纲”级别。具体到哪个层级，取决于你数据库的覆盖度和 k‑mer 的唯一性。我一般建议先看“属”级结果，因为很多近缘种的 k‑mer 太像了，容易混淆。
举个例子：Escherichia coli 和 Shigella flexneri 在 16S rRNA 上可能差别很小，但 Kraken2 如果用了完整的基因组 k‑mer 索引，就能在 k‑mer 层面区分开。我在一个项目里遇到过，某次跑临床样本，大量序列被标成 E. coli，后来发现其实是数据库里 Shigella 的 k‑mer 太少，导致归类错误。避坑指南就是——数据库的完整性和代表性比啥都重要。
| 层级 | 例子 | Kraken2 中的意义 |
|---|---|---|
| 域 (Domain) | 细菌、古菌、真核 | 最顶层，通常不用做最终分类 |
| 门 (Phylum) | Proteobacteria | 样本中高丰度门可快速识别 |
| 纲 (Class) | Gammaproteobacteria | 区分大类时有用 |
| 目 (Order) | Enterobacterales | 中间层级，辅助判断 |
| 科 (Family) | Enterobacteriaceae | 常与属一起看 |
| 属 (Genus) | Escherichia | 我最推荐的报告单位 |
| 种 (Species) | E. coli | 准确但依赖数据库 |
| 株 (Strain) | O157:H7 | 需要专门数据库才能分辨 |
2. k‑mer 原理——核心中的核心
什么是 k‑mer？就是一个长度为 k 的短 DNA 子串。比如序列 ATGCGTA，如果 k=3，那它的 3‑mer 就是 ATG、TGC、GCG、CGT、GTA 这五个。Kraken2 干的事情很简单：把参考基因组（细菌、病毒、真菌等）全部切成 k‑mer，建一个巨大的哈希表。然后你拿一条测序 reads 过来，也切成同样的 k‑mer，去哈希表里查，看它落在哪个分类节点上。
为什么要用 k‑mer 而不是直接比对？说白了，速度。全局比对（比如 BLAST）太慢了，不可能处理几千万条 reads。而 k‑mer 查哈希表是 O(1) 的复杂度，加上 Kraken2 用了最小万花筒（Minimizer）技术压缩索引，内存和速度都控制得很好。我做过一个 100G 的肠道样本，Kraken2 只用 8 分钟就跑完了，而 BLAST 可能要跑几天。
小技巧：k 值的选择很关键。Kraken2 默认 k=31（对于细菌），因为 31‑mer 在细菌基因组中通常有足够的唯一性。如果你分析的是病毒或者很短的 reads（比如 50bp），我建议把 k 调小到 21 或 25，否则很多 k‑mer 在数据库里根本没出现过。
下面是一个简单的 Python 脚本，演示如何从一条序列生成 k‑mer：
def get_kmers(seq, k):
kmers = []
for i in range(len(seq) - k + 1):
kmers.append(seq[i:i+k])
return kmers
seq = "ATGCGTA"
k = 3
kmers = get_kmers(seq, k)
print(kmers)
# 输出: ['ATG', 'TGC', 'GCG', 'CGT', 'GTA']不过，实际 Kraken2 里不会保存所有 k‑mer，而是用 minimizer（最小化 k‑mer）来节约内存。每个 k‑mer 会选出一个最短的、字典序最小的子串作为代表。我曾踩过一个坑：用默认的 minimizer 长度 15，结果在高度重复的区域，不同物种的 minimizer 完全相同，导致分类错误。后来我改用更长的 minimizer（比如 18），精度立马提升了。
3. Kraken2 中的 k‑mer 索引流程
Kraken2 建数据库的完整流程可以概括为：
- 把每个参考基因组的染色体序列切成长度为 k 的 k‑mer（k 默认 31）
- 对每个 k‑mer，找到它的 minimizer（默认长度 15）
- 以 minimizer 为键，存储 (k‑mer, 分类ID) 对到哈希表
- 同一个 minimizer 对应的多个 k‑mer 会被合并，降低索引大小
你运行 kraken2-build 时，它会自动做上面这些事。但如果你自己定制数据库，比如只包含某个门的细菌，那需要注意 minimizer 参数要跟运行 kraken2 时的一致，否则索引没法用。
下面是一个我曾经踩过的坑，写成 warning 分享给大家：
曾经我犯过的错：用 kraken2-build --download-library bacteria --kmer-len 31 --minimizer-len 15 建库，但在跑分类时忘了指定一样的参数，结果 Kraken2 报了一堆 "incompatible index" 错误。后来我养成了习惯——建库和运行用同一份配置文件，或者在脚本里用变量统一参数。
另外，k‑mer 的数量跟基因组大小关系很大。一个细菌基因组大约 5 Mbp，k=31 时产生的 k‑mer 数大约是 5M 个（去除重复后约 1‑2M 个唯一 k‑mer）。人类基因组有 3 Gbp，k‑mer 数量是几十亿级别，所以 Kraken2 通常只用于微生物，因为人类基因组太大了，索引会爆内存。
上面这张图就是 Kraken2 从读段到分类节点的完整流程。你仔细看，从切 k‑mer 到计算 minimizer，再到查哈希表，每一步都是影响速度和精度的关键点。
总结一下核心要点
- 分类层级：掌握从域到株的七级体系，理解 Kraken2 输出中不同层级的意义。我在项目中坚持用“属”级作为默认报告，因为种级假阳性偏多。
- k‑mer 原理：k 和 minimizer 长度是调试参数的灵魂。k 越大，唯一性越好但覆盖度越低；k 越小，覆盖度越高但噪声越多。我通常对 150bp 双端 reads 用 k=31，对 50bp 单端用 k=21。
- 索引一致性：建库和运行必须用同样的 k 和 minimizer 长度，否则直接报错。这个坑我掉进去过两次，现在用 Makefile 统一管理参数，再也没出过问题。
嗯，到这里分类学层级和 k‑mer 原理就讲完了。下节课我们会直接上手安装 Kraken2 并构建第一个数据库。但在此之前，我强烈建议你把上面的流程图在脑子里过一遍，想清楚每个 k‑mer 是怎么变成分类节点的——这对你后面调试参数、评估结果质量会非常有帮助。