MetaPhlAn 原理：基于标记基因的分类方法
上一章聊了 Kraken2 怎么用 k-mer 打标签，这次咱们换个思路——MetaPhlAn。
说实话，我刚接触这工具时觉得挺神奇的。你想想看，同样是做宏基因组物种鉴定，Kraken2 把整个基因组都拿来比对，而 MetaPhlAn 只看几个“关键特征”——这就是标记基因的思路。
标记基因是什么？为什么选它？
标记基因，就是那些在物种内部高度保守、在物种之间又有足够差异的基因。比如 16S rRNA 基因。但 MetaPhlAn 用的不光是 16S，它从每个物种的全基因组里挑出一批单拷贝、物种特异、能唯一标识的基因。
我前些年帮一家医院做肠道菌群项目，当时老板要求同时跑 Kraken2 和 MetaPhlAn。结果发现，有些低丰度物种 Kraken2 给出好几个矛盾的界级分类，MetaPhlAn 却能稳定地识别到属。为什么会这样？因为标记基因的筛选原则就是避免交叉映射。
MetaPhlAn 的标记基因库筛选了大约 1.3 万个基因组，每个物种只保留几百个基因。数据库小，比对快。说白了，它只回答一个问题：哪些物种的标记基因片段出现在我的样本里？
核心流程：标定、比对、丰度估计
MetaPhlAn 的工作流可以分为三步，我习惯用这张图来团队讲解：
┌───────────────────────────────────────────────────┐
│ MetaPhlAn 工作流程 │
│ reads (FASTQ) │
│ │ │
│ ▼ │
│ 比对到标记基因数据库 (Bowtie2) │
│ │ │
│ ▼ │
│ 保留唯一比对到标记基因的 reads │
│ │ │
│ ▼ │
│ 基于标记基因计数，用贝叶斯/EM 推算物种相对丰度 │
│ │ │
│ ▼ │
│ 输出：物种丰度表 (TSV) │
└───────────────────────────────────────────────────┘这个流程分三步：
- 比对：用 Bowtie2 把 clean reads 往 MetaPhlAn 标记基因数据库里怼。这里有个关键——只保留唯一比对到某个标记基因的 reads。那些能同时比对上多个标记基因的 reads，会被直接扔掉。我曾经见过一个样本，30% 的 reads 因为多匹配被剔除，幸好数据库更新后比例降到了 5% 以下。
- 计数：统计每个物种的标记基因上覆盖了多少 reads。注意不是简单加和，因为有基因长度、拷贝数差异，需要归一化。
- 丰度估计：这是数学活。MetaPhlAn 用 EM 算法（Expectation-Maximization）来估算每种菌的相对丰度。原理很简单：根据所有标记基因的 reads 分布，通过迭代找到最可能的物种组成。我在项目里调试过这个步骤，发现如果样本某个物种丰度极低（比如 <0.1%），EM 算法可能收敛不到全局最优。这时候我会手动调高–-min_alignment_length 参数。
--add_viruses 和 --add_archaea 参数，覆盖更全。
标记基因库的进化
MetaPhlAn 数据库从 v2.0 到 v4.0 翻了几番。v2.0 只覆盖约 3500 个物种，v4.0 已经包含 1.3 万个物种、30 万个标记基因。每次更新我都会跑一遍 metaphlan --install，因为旧库会漏掉很多新测的菌株。
| 版本 | 标记基因数 | 覆盖物种 | 数据库大小 |
|---|---|---|---|
| v2.0 | ~2.5 万 | ~3500 | ~200 MB |
| v3.0 | ~7 万 | ~1 万 | ~600 MB |
| v4.0 | ~30 万 | ~1.3 万 | ~2.5 GB |
数据库大了，比对时间会变长。不过 MetaPhlAn 只比对标记基因，所以速度仍然比 Kraken2 的完整基因组索引快 5-10 倍。我实测过，一个 10M reads 的样本，用 8 核机器只跑 3 分钟。
实际跑一个看看
假设你已经安装了 MetaPhlAn 4（推荐用 conda：conda install -c biobakery metaphlan）。
命令很简洁：
metaphlan input.fastq --bowtie2db /path/to/metaphlan_databases --input_type fastq -o profile.txt输出文件 profile.txt 长这样：
#clade_name NCBI_tax_id relative_abundance coverage estimated_number_of_reads
k__Bacteria 2 99.85 29.12 1012345
k__Bacteria|p__Firmicutes 1239 45.12 32.54 456789
k__Bacteria|p__Firmicutes|c__Clostridia 186801 30.23 27.89 305678
k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales 816 28.56 26.11 289123
... ... ... ... ...每一行是一个分类等级的丰度信息，从界（k__）到种（s__）。你可以用 --profile_relative_abundance 控制输出格式。
--bowtie2out 把比对结果暂存到磁盘，或者升级到 16GB 内存。
MetaPhlAn vs Kraken2：选哪个？
两者不是一个赛道的。Kraken2 是“看见什么就喊什么”，MetaPhlAn 是“只认熟面孔”。我做了一张对比表，方便你快速决策：
| 维度 | MetaPhlAn | Kraken2 |
|---|---|---|
| 核心策略 | 基于标记基因的唯一比对 | 基于 k-mer 的精确匹配 |
| 数据库大小 | 较小（2.5 GB） | 较大（标准库 5-20 GB） |
| 速度 | 快（约 3 min/10M reads） | 极快（约 1 min/10M reads） |
| 对未知物种 | 无法识别（标记基因库必须包含） | 可给出较低层级分类 |
| 定量方式 | EM 算法产生相对丰度 | 直接计数后归一化 |
| 菌株级分辨率 | 支持（部分菌株） | 支持（需要全基因组索引） |
| 典型场景 | 已知菌群的精准量化 | 快速初筛、非标记物种检测 |
我个人的习惯是：先用 Kraken2 做泛筛查，再用 MetaPhlAn 做精准定量。 比如一个环境样本，我先跑 Kraken2 看看有哪些门类，把明显的人为污染或未知序列剔除。然后针对关注的菌群，用 MetaPhlAn 输出准确的丰度矩阵。这样既快又准。
SVG: MetaPhlAn 分类流程示意图
总结几个实用经验
- 数据库安装一定要用最新版： MetaPhlAn 数据库每月更新一次。我去年遇到过一个假阳性案例——因为库里的标记基因有冗余，导致把 E. coli 误判为 Shigella。更新后问题消失。
- 记得做质量控制： 跑 MetaPhlAn 前至少用 fastp 或 Trimmomatic 去除低质量 reads。如果样本包含大量宿主 DNA，我建议先用
kraken2过滤掉宿主序列，再把剩余 reads 喂给 MetaPhlAn。这样能节省 30-50% 的比对时间。 - 结果怎么看： 输出文件里
estimated_number_of_reads列是关键。如果某物种的 reads 数少于 100，统计可靠性下降。标记基因越多的物种，该值越可信。 - 当心病毒： MetaPhlAn 默认不包含病毒标记基因。需要加
--add_viruses，并下载专门的病毒数据库。我曾因此漏掉了样本里的噬菌体，后来补跑了一次。