第3章 读段比对与物种分类归纳总结
这一章我们聊聊宏基因组分析里最核心的一步：把测序读段比对到参考数据库，然后归纳出每个物种的丰度。我个人刚开始做宏基因组时，觉得这一步就是跑个命令而已，结果各种坑踩下来才发现，这里面的门道比想象中多得多。
3.1 读段比对——别急着跑命令
说白了，读段比对就是为每一条测序读段找到一个「标尺」。这个标尺就是参考基因组或基因数据库。你想想看，宏基因组样本里混杂着数百种微生物的DNA片段，就像一锅大杂烩。比对的目的就是把每一块碎片归类到它原来的物种身上。
我建议大家在选择比对工具时，先问自己三个问题：
- 读段长度是多长？短读（<150bp）和长读（>500bp）用的工具不同
- 参考数据库有多大？几GB还是几百GB？这直接影响运行时间和内存
- 要的是物种级分辨率还是功能级？不同问题选择不同策略
3.2 常用比对工具对比
下面这张表是我自己整理的，是这么多年项目里反复试错后的经验总结。
| 工具 | 比对策略 | 适用场景 | 优缺点 |
|---|---|---|---|
| Bowtie2 | FM-index 双索引比对 | 短读段到参考基因组的快速比对 | 速度极快，但对高度相似区域区分力弱 |
| BWA-MEM | Burrows-Wheeler 变换 | 较长读段或同源性搜索 | 灵敏度和特异性均衡，适合全长基因组 |
| Kraken2 | k-mer 精确匹配 + LCA 算法 | 宏基因组学快速物种分类 | 准确度高，但需要较大内存（~50GB+） |
| Bracken | 基于 Kraken2 结果的贝叶斯校正 | 物种到属级的丰度估计 | 比 Kraken2 原始计数更接近真实丰度 |
Bowtie2是我个人比较爱用的一款短读段比对器。它专为宏基因组而生，比对速度快，内存占用也合理。我曾经用它在一天之内处理了200个环境样本的比对任务，那台服务器也只有64GB内存，跑下来非常顺利。
3.3 物种分类策略
分类策略本质上分两大派：
- 有参考基因组派：先把读段比对到已知物种的基因组，然后统计每个基因组上的覆盖深度和覆盖比例。这种方法对已知物种效果好，但遇到新物种就抓瞎。
- 无需参考基因组派：用k-mer方法直接对读段进行分类，Kraken2就是代表。它把每条读段拆成短k-mer，然后去参考数据库里找精确匹配，再用最长公共祖先算法来确定分类级别。
我个人习惯优先使用Bowtie2进行初步比对，过滤掉宿主和低质量匹配后，再对剩余读段用Kraken2做精细分类。这两步组合起来，能覆盖宏基因组样本中大部分已知和未知微生物。
--very-sensitive-local参数虽然慢，但能提高多拷贝基因的检出率。我在处理重复序列较多的样本时，用这个参数总能有意外发现。
3.4 归纳总结实战
比对完成后的归纳总结，其实关键不在于跑哪个工具，而在于跑完之后你怎么整理。我见过很多同行跑完命令就把结果丢在那里不管了，这样容易漏掉重要信息。
下面是一个典型的Bowtie2比对命令：
bowtie2 -x /ref_db/hg38 -1 sample_R1.fastq -2 sample_R2.fastq \
-S sample_aligned.sam \
--very-sensitive \
-p 16 \
2> sample_bowtie2.log这是Kraken2分类调用示例：
kraken2 --db /kraken_db/standard \
--threads 16 \
--output sample_kraken_output.txt \
--report sample_kraken_report.txt \
--use-names \
sample_R1.fastq sample_R2.fastq接着用Bracken做丰度校正：
bracken -d /kraken_db/standard \
-i sample_kraken_report.txt \
-o sample_bracken_species.txt \
-r 150 \
-l S-r 150必须和你的实际读段长度一致。我上次帮同事排查问题，发现他跑出来结果异常，就是因为读段是100bp，但他写成了150bp。Bracken内部会根据读段长度调整k-mer覆盖度估计，参数不对结果白搭。
结果归纳时，我建议做三件事：
- 生成物种丰度表（行是物种，列是样本）
- 计算α多样性（Shannon指数、Chao1等）
- 做β多样性分析（主坐标分析PCoA）
举个例子，从Bracken的物种丰度表里，你可以直接看每个物种的相对丰度。比如某环境样本中，假单胞菌属占35%，芽孢杆菌属占28%，这都是可以直接读出来的信息。
3.5 实战建议
你可能会问，这么多工具，我到底选哪个？别急，我给一个快速决策指南：
- 样本数量多（>50个），优先选Kraken2+Bracken，速度快，内存优化好
- 需要高精度分析，比如临床诊断，用Bowtie2+BWA二阶段比对，再做LCA分类
- 长读段样本（PacBio/Nanopore），用minimap2+Bracken的修改版
我曾经在处理一组100个临床样本的时候，一开始用了Bowtie2全文库比对，结果跑了一周都没跑完。后来换了Kraken2的快速模式，四小时搞定。选择工具真要考虑数据规模。
最后，别忘了比对后的日志文件。那里藏着比对率、多匹配率、未比对率等关键指标。如果比对率低于10%，说明样本质量或者数据库选择有问题，需要回头检查。
这一章就告一段落。希望这些实战经验能帮你少走弯路。比对分类做得好，后续分析就有了扎实的基础。