3. 输出文件解析：理解 profile、biom、sam 等输出格式及其字段含义
跑完 MetaPhlAn，你会拿到一堆文件。新手最容易犯的错——盯着屏幕看半天，不知道哪个文件是干嘛的。
我个人习惯，先搞清楚每种格式的定位，再去抠细节。这样内存占用小，效率高。
3.1 核心输出总览
MetaPhlAn 默认生成三个主要文件，外加一个可选的中间文件。我把它们的关系画了张图，你看完心里就有谱了。
3.2 Profile 文件——最常用的结果文件
这是你的主力文件。默认文件名是 sample_name_profile.txt。说白了，它就是个物种丰度表，但带了完整的分类学路径。
3.2.1 列字段详解
| 列名 | 含义 | 避坑提示 |
|---|---|---|
clade_name |
分类学谱系，从界到种的全路径 | 用 | 分隔，别漏看了层级 |
NCBI_tax_id |
NCBI 分类学 ID | 不同数据库版本可能不同 |
relative_abundance |
该分类的相对丰度（百分比） | 所有行加起来等于 100 |
estimated_number_of_reads |
该分类对应的估计 reads 数 | 是估计值，不是精确计数 |
additional_species |
如果当前节点是属，会列出下属的种 | 种级分辨率不够时参考 |
我的小技巧：拿到 profile 后，先看 relative_abundance 最高的那几行是不是 contaminants。我去年分析一个土壤项目，发现第一行是 Pseudomonas 占了70%，结果是因为水样被污染了。查了 lab log 才确认。
3.2.2 示例片段
#clade_name NCBI_tax_id relative_abundance estimated_number_of_reads additional_species
k__Bacteria|p__Proteobacteria|g__Pseudomonas 286 45.30 45200
k__Bacteria|p__Firmicutes|g__Clostridium 1506 22.10 22100
k__Bacteria|p__Bacteroidetes|g__Bacteroides 816 15.80 15800
k__Archaea|p__Euryarchaeota|g__Methanobrevibacter 2172 0.05 50
3.3 BIOM 文件——标准交换格式
BIOM 其实是 JSON 格式，结构比较复杂。不过我建议你别手动去编辑它，用 biom 命令行工具或者 Python 的 biom-format 包来读写。
3.3.1 核心结构
一个 BIOM 文件包含以下层级：
- matrix_type：矩阵类型，通常是 "dense" 或 "sparse"。
- shape：矩阵维度，[物种数量, 样本数量]。
- data：丰度数据，按列展开或按稀疏格式存储。
- rows：观测元数据（物种信息）。
- columns：样本元数据。
注意：MetaPhlAn 生成的 BIOM 文件默认只包含一种数据类型（OTU 表），不包含功能注释。如果需要功能信息，请使用 --add_viruses 或 --add_archaea 等扩展参数。
3.3.2 通过命令查看
biom summarize-table -i output.biom -o table_summary.txt你会看到类似这样的输出：
Num observations: 3421
Num features: 1
Num samples: 1
Total count: 100000.0
Table density: 100%嗯，这里要注意——密度 100% 说明每个物种都有非零丰度，但在宏基因组里不太可能。其实是因为 MetaPhlAn 把所有检测到的物种都平铺进去了。
3.4 SAM 文件——比对原始数据
这个文件大部分情况下你不需要碰。我只有在调试参数或者怀疑比对出问题时才会打开它。
3.4.1 SAM 格式的核心列
| 列号 | 字段 | 说明 |
|---|---|---|
| 1 | QNAME | 查询序列名称（read ID） |
| 2 | FLAG | 比对标志位（0 = 唯一比对，256 = 次优选） |
| 3 | RNAME | 参考序列名称（标记基因 ID） |
| 4 | POS | 比对位置（1-based） |
| 5 | MAPQ | 比对质量值（越高越可信） |
| 6 | CIGAR | 比对操作字符串（如 150M） |
| 10 | SEQ | 序列碱基 |
| 11 | QUAL | 质量值字符串 |
3.4.2 有用的标记（Tags）
SAM 文件的第 12 列开始是一些 key-value 对。MetaPhlAn 会写入几个重要的 tag：
- NM：编辑距离（错配数）。值越低越好。
- AS：比对得分。比对软件保留下来的。
- X1：MetaPhlAn 的标记基因分类信息。
我曾经犯过一个低级错误：用 SAMtools 统计唯一比对 reads 数时，没过滤掉 flag=256 的次优选。结果多算了一倍的 reads，丰度全乱套。记住——只统计 flag=0 或 flag=16 的。
3.5 如何选择输出格式
我的选择逻辑很简单，你参考一下：
- 日常分析：只用 profile 文件，文本格式好处理。
- 和同事/工具交换数据：用 BIOM 格式，标准化程度高。
- 调试性能或者查 bug：保留 SAM 文件，用
--sam_output output.sam指定路径。
核心原则：别在不需要的时候生成 SAM 文件。我见过有人每次跑都带 --sam_output，结果磁盘被撑爆了。一个 10M reads 的样本，SAM 文件能到 5GB 以上。
3.6 实战：用 Python 解析 Profile 文件
如果你喜欢用脚本处理，这段代码帮你快速上手：
import pandas as pd
df = pd.read_csv("profile.txt", sep="\t", skiprows=1)
df.columns = ["clade_name", "NCBI_tax_id",
"relative_abundance", "estimated_number_of_reads",
"additional_species"]
# 只看种水平
species = df[df["clade_name"].str.count("\|") == 6]
# 找 top10
top10 = species.nlargest(10, "relative_abundance")
print(top10[["clade_name", "relative_abundance"]])
运行后会显示丰度最高的10个物种，直接拿去画柱状图就完事了。
3.7 避坑清单
- Profile 文件第一行是注释，以
#开头，读取时记得跳过。 - 相对丰度加总不一定正好等于100，因为有未分类的 reads，浮点误差也可能。
- BIOM 文件用 gzip 压缩，别直接
less打开。可以用zcat output.biom | python -m json.tool看结构。 - SAM 文件不要随便用
grep，用samtools view -hS来读，否则 header 会乱。
明白了这些文件长什么样，你的 MetaPhlAn 之旅就成功了一半。剩下的一半，是让参数替你把活干好。