总体流程
HUMAnN (HMP Unified Metabolic Analysis Network) 是 bioBakery 工作流中用于分析宏基因组功能的核心工具[[12]]。通路丰度的计算是一个多步骤的递归过程：
计算步骤：
- 基因家族丰度 → 2. 反应丰度 → 3. 通路丰度
详细计算原理
第1步：基因家族丰度（Gene Family Abundance）
从原始测序 reads 开始：
- 使用 BLASTX 将 reads 比对到参考数据库（如 UniRef）
- 根据比对质量、覆盖度、序列长度进行加权
- 生成 RPK（Reads Per Kilobase）值
公式：
基因丰度 = Σ(比对权重) / 基因长度(kb)其中每个 read 的总权重为 1.0，根据比对质量分配到多个基因匹配上[[9]]。
第2步：反应丰度（Reaction Abundance）
每个生化反应由一个或多个基因催化：
反应丰度 = Σ(催化该反应的所有基因丰度)第3步：通路丰度（Pathway Abundance）
这是最关键的一步。通路包含多个反应，反应之间有不同的关系：
核心原则： 通路丰度由”最弱环节”（weakest link）决定[[1]]
计算方法：
- 串联反应（必须全部存在）：使用调和平均数（harmonic mean）
- 并联反应（可选路径）：使用最大值（max）
- 可选反应：只有当其丰度大于必需反应的调和平均数时才计入[[1]]
最终通路丰度 = 通路中丰度最低的关键反应
具体示例
示例场景：糖酵解通路（Glycolysis）
假设糖酵解通路包含 5 个关键反应（R1-R5）：
葡萄糖 → R1 → G6P → R2 → F6P → R3 → ... → 丙酮酸基因-反应关系：
- R1: 由基因 GK1 和 GK2 催化（冗余）
- R2: 由基因 PGI 催化
- R3: 由基因 PFK 催化
- R4: 由基因 ALDO 催化
- R5: 由基因 GAPDH 催化
测序后得到的基因丰度（RPK单位）：
GK1: 8.0
GK2: 4.0
PGI: 10.0
PFK: 6.0
ALDO: 7.0
GAPDH: 5.0计算步骤：
① 计算反应丰度：
R1 = GK1 + GK2 = 8.0 + 4.0 = 12.0 (冗余基因相加)
R2 = PGI = 10.0
R3 = PFK = 6.0
R4 = ALDO = 7.0
R5 = GAPDH = 5.0② 计算通路丰度： 由于糖酵解是串联反应（所有步骤必须完成），使用”最弱环节”原则：
通路丰度 = min(R1, R2, R3, R4, R5)
= min(12.0, 10.0, 6.0, 7.0, 5.0)
= 5.0 RPK解释： 该样本中糖酵解通路的丰度为 5.0 RPK，意味着”最弱环节”（R5/GAPDH）的覆盖度为 5.0。这表示通路中至少有 5.0 个”完整拷贝”的活性[[1]]。
归一化处理
为什么需要归一化？
原始 RPK 值受测序深度影响，不能直接跨样本比较[[1]]。
示例：
- 样本 A：总 reads = 1000万，通路丰度 = 5.0 RPK
- 样本 B：总 reads = 2000万，通路丰度 = 5.0 RPK
虽然都是 5.0 RPK，但样本 A 的相对丰度更高！
归一化方法：
CPM（Counts Per Million）或 RPKM：
归一化丰度 = (原始RPK / 总RPK) × 1,000,000HUMAnN 输出的 pathabundance_relab.tsv 文件已经是归一化后的相对丰度[[11]]。
输出文件说明
HUMAnN 生成两个关键文件：
| 文件 | 含义 | 取值范围 |
|---|---|---|
pathcoverage.tsv |
通路覆盖度（存在/缺失概率） | 0-1 |
pathabundance.tsv |
通路丰度（相对拷贝数） | 0 及以上 |
覆盖度 vs 丰度：
- 覆盖度 = 0.8：表示 80% 置信度该通路存在
- 丰度 = 5.0：表示该通路约有 5 个完整拷贝（在归一化前）[[9]]
分层分析（Stratification）
HUMAnN 还能计算每个物种对通路的贡献：
总通路丰度 = 物种A贡献 + 物种B贡献 + 未分类贡献例如：
糖酵解通路总丰度: 5.0 RPK
├─ Escherichia coli: 3.0 RPK (60%)
├─ Bacteroides fragilis: 1.5 RPK (30%)
└─ Unclassified: 0.5 RPK (10%)这帮助您理解哪些微生物在驱动特定功能[[21]]。
关键要点总结
✅ 通路丰度 = 最弱关键反应的丰度（木桶原理）
✅ 单位：RPK（归一化前）或 CPM/RPKM（归一化后）
✅ 跨样本比较：必须使用归一化值（pathabundance_relab.tsv）[[15]]
✅ 解释：丰度 5.0 表示通路中最弱环节有 5.0 RPK 的覆盖度
✅ 准确性：考虑了基因冗余、反应关系、可选步骤等生物学复杂性
参考文献
- bioBakery Forum – Pathway abundance calculation [[1]]
- HUMAnN SOP – HMP Data Coordination Center [[9]]
- HUMAnN3 Documentation – Huttenhower Lab [[12]]
- bioBakery 3 publication (Nature Methods, 2021) [[13]]
HUMAnN Pathway Abundance 计算方法
基本原理
HUMAnN (HMP Unified Metabolic Analysis Network) 是 biobakery 工具集的一部分，用于从宏基因组数据中计算功能通路的丰度。
计算步骤
1. 基因家族定量 (Gene Family Quantification)
- 首先将测序 reads 比对到 UniRef 基因家族数据库
- 计算每个基因家族的丰度（reads per kilobase, RPK）
2. 通路映射 (Pathway Mapping)
- 将基因家族映射到 MetaCyc 通路
- 一个通路通常包含多个反应步骤，每个步骤可能由多个基因家族催化
3. 通路丰度计算 (Pathway Abundance Calculation)
关键公式：
通路丰度 = min(该通路中所有反应步骤的丰度)其中：
- 每个反应步骤的丰度 = 该步骤中所有基因家族丰度的总和
- 通路的最终丰度 = 所有反应步骤丰度的最小值（瓶颈原则）
具体示例
假设情况
假设有一个简单的代谢通路 “Glycolysis”（糖酵解），包含 3 个反应步骤：
反应步骤 1：葡萄糖 → 葡萄糖-6-磷酸
- 由基因家族 UniRef90_A 和 UniRef90_B 催化
- UniRef90_A 丰度 = 100 RPK
- UniRef90_B 丰度 = 50 RPK
- 步骤 1 丰度 = 100 + 50 = 150 RPK
反应步骤 2：葡萄糖-6-磷酸 → 果糖-6-磷酸
- 由基因家族 UniRef90_C 催化
- UniRef90_C 丰度 = 80 RPK
- 步骤 2 丰度 = 80 RPK
反应步骤 3：果糖-6-磷酸 → 果糖-1,6-二磷酸
- 由基因家族 UniRef90_D 和 UniRef90_E 催化
- UniRef90_D 丰度 = 200 RPK
- UniRef90_E 丰度 = 120 RPK
- 步骤 3 丰度 = 200 + 120 = 320 RPK
通路丰度计算
Glycolysis 通路丰度 = min(步骤1, 步骤2, 步骤3)
= min(150, 80, 320)
= 80 RPK为什么用最小值？
- 这遵循”木桶原理”（瓶颈效应）
- 通路的整体通量受限于最慢/最少的反应步骤
- 步骤 2 只有 80 RPK，即使其他步骤丰度很高，整个通路的最大通量也只能是 80
归一化 (Normalization)
HUMAnN 还会进行归一化处理：
1. RPK (Reads Per Kilobase)
RPK = (比对到基因的 reads 数) / (基因长度 kb)2. RPKM/CPM (可选)
RPKM = RPK / (总百万 reads)3. 相对丰度 (Relative Abundance)
通路相对丰度 = (通路 RPKM) / (所有通路 RPKM 总和) × 1,000,000单位：CPM (Counts Per Million)
您的数据示例
从您的输出看：
otu_table() OTU Table: [ 1097 taxa and 2 samples ]这 1097 个 “taxa” 实际上是 1097 个 MetaCyc 通路，每个通路的丰度值是经过以下处理的：
- 原始丰度：基于上述最小值原则计算
- 归一化：转换为相对丰度（sum = 1 或 100%）
- 输出文件：
pathabundance_relab.tsv中的值就是相对丰度
关键要点
| 特点 | 说明 |
|---|---|
| 计算方法 | 取通路中所有反应步骤的最小丰度 |
| 单位 | 通常是相对丰度（0-1 或 0-100%） |
| 生物学意义 | 反映通路的潜在代谢能力 |
| 优势 | 考虑了通路的完整性，不是简单加和 |
| 局限性 | 无法区分活跃/非活跃通路（需要转录组验证） |
注意事项
⚠️ 重要提醒：
- 通路丰度反映的是基因潜力（gene potential），不是实际代谢活性
- 一个通路存在 ≠ 该通路正在被使用
- 需要结合转录组（RNA-seq）或代谢组数据才能确定实际活性
- 对于您的 n=1 样本，只能做描述性比较，无法统计推断