4. 读取BIOM格式：数据结构初探
BIOM格式在微生物组领域几乎成了标准。不管你是做16S还是宏基因组，最后拿到的丰度表八成是.biom文件。
我个人第一次接触BIOM时有点懵。打开一看，全是JSON格式的嵌套结构。后来慢慢摸透了它的脾气，才发现其实挺清晰。
4.1 准备工作：安装biom-format
用Python读取BIOM文件，首选biom-format库。安装很简单：
pip install biom-format如果你在用conda环境（我强烈建议你这么做），可以这样装：
conda create -n biom_env python=3.9
conda activate biom_env
conda install -c conda-forge biom-format
小提示：
建议用Python 3.8以上版本。我曾经在Python 3.6上遇到过一个奇怪的编码问题，折腾了半天才发现是版本兼容性导致的。
4.2 BIOM v1.0 与 v2.1 的区别
BIOM有两个主版本。说白了就是：v1.0是纯JSON，v2.1是JSON加HDF5。
| 特性 | BIOM v1.0 | BIOM v2.1 |
|---|---|---|
| 存储格式 | 纯JSON文本 | JSON + HDF5 |
| 文件大小 | 较大（明文存储） | 较小（二进制压缩） |
| 读取速度 | 慢（全量加载） | 快（支持部分读取） |
| 适合场景 | 小数据集、调试 | 大规模数据、生产环境 |
你想想看，如果你的数据只有几十个样本，v1.0完全够用。但如果样本数上千，v2.1的HDF5格式能省下大量内存和时间。
4.3 读取BIOM文件：基础操作
读取BIOM文件其实就一行代码。但我还是建议你走完整流程，先看看文件头信息。
import biom
# 读取v1.0格式
table_v1 = biom.load_table('example_v1.biom')
# 读取v2.1格式（biom-format会自动识别）
table_v2 = biom.load_table('example_v2.biom')
# 查看基本信息
print(table_v1)
print(f"形状: {table_v1.shape}")
print(f"样本数: {table_v1.shape[0]}")
print(f"特征数（OTU/ASV）: {table_v1.shape[1]}")输出看起来像这样：
BIOM Table: 100 samples x 500 features
使用稠密矩阵存储
形状: (100, 500)
样本数: 100
特征数（OTU/ASV）: 500
核心概念：
BIOM对象的shape是(样本数, 特征数)。和丰度表转置后的概念一致。很多新手在这里搞反了行列，千万注意。
4.4 观察内部数据结构
BIOM对象内部存储了丰度矩阵、样本元数据和特征元数据。我们一层层拆开看看。
4.4.1 矩阵数据
# 获取丰度矩阵（稠密形式）
dense_data = table_v1.matrix_data.toarray()
print(dense_data[:5, :5]) # 前5行5列
# 获取稀疏矩阵（推荐大文件用）
sparse_data = table_v1.matrix_data
print(f"稀疏矩阵类型: {type(sparse_data)}")
注意：
对于大文件，千万别直接toarray()。我曾在处理一个1000样本×5000特征的文件时，直接toarray()导致内存爆了。建议用稀疏矩阵做计算。
4.4.2 样本ID和特征ID
# 样本ID
sample_ids = table_v1.ids(axis='sample')
print(f"前5个样本ID: {sample_ids[:5]}")
# 特征ID（OTU/ASV ID）
feature_ids = table_v1.ids(axis='observation')
print(f"前5个特征ID: {feature_ids[:5]}")4.4.3 元数据
# 样本元数据
sample_metadata = table_v1.metadata(axis='sample')
print(f"第一个样本的元数据: {sample_metadata[0]}")
# 特征元数据
feature_metadata = table_v1.metadata(axis='observation')
print(f"第一个特征的元数据: {feature_metadata[0]}")每个样本的元数据是个字典。比如：
{
'SampleID': 'S001',
'Group': 'Treatment',
'TimePoint': 'Day0',
'Description': '基线样本'
}4.5 实际项目中的坑
我曾经在一个项目中遇到这么个情况：读取BIOM文件后，发现丰度全为零。查了半天才发现，原来文件是用旧版Qiime生成的，存储格式有问题。
有几个常见坑你提前知道比较好：
- 索引不一致： 样本ID和特征ID的顺序，在不同BIOM文件中可能不同。做合并操作时务必对齐。
- 缺失值处理： BIOM规范允许缺失值（NaN）。你读取后最好检查一下。
- 数据类型： 丰度值可能是int、float甚至string。我见过某工具输出全是字符串的BIOM，坑得不行。
# 检查缺失值
import numpy as np
data = table_v1.matrix_data.toarray()
print(f"缺失值数量: {np.isnan(data).sum()}")
# 检查数据类型
print(f"数据类型: {data.dtype}")4.6 BIOM数据结构总览
下面这张图能帮你快速理解BIOM的结构：
4.7 实用小工具：快速查看BIOM内容
有时候你只是想知道BIOM文件里有什么，不想写代码。可以用命令行工具：
biom summarize-table -i example.biom -o summary.txt输出像这样：
Num samples: 100
Num observations: 500
Total count: 1250000
Table density: 0.15 (15% 非零值)
实用技巧：
用biom summarize-table查看密度(density)。如果密度很低（<5%），说明数据很稀疏，后续分析可能需要考虑零膨胀模型。
还有个命令能直接看元数据：
biom head -i example.biom -n 5这个命令会打印前5行样本的丰度和元数据，调试时特别好用。
4.8 总结一下
读取BIOM文件其实不难。关键是要理解它的三层结构：丰度矩阵、样本元数据、特征元数据。
我个人习惯是先看summary，再看元数据。等确认数据质量没问题后，才做后续分析。你想想看，如果源数据就有问题，后面一堆漂亮图表都是白搭。