跳转至

宏基因组代谢通路重建

一句话概述

使用KEGG module completeness方法评估宏基因组中微生物代谢通路的完整性，通过HUMAnN3/MinPath/DRAM等工具重建社区级别的代谢能力图谱，揭示微生物群落的功能潜力。

核心知识点总览

知识点关键内容重要程度
KEGG Module概念功能单元/模块完整性定义⭐⭐⭐⭐⭐
Module completeness通路步骤的覆盖度评估⭐⭐⭐⭐⭐
HUMAnN3社区级代谢通路重建⭐⭐⭐⭐⭐
DRAM基因组/MAG代谢注释⭐⭐⭐⭐
MinPath最小通路集推断⭐⭐⭐⭐
KofamScanKEGG KO注释⭐⭐⭐
通路丰度归一化CPM/RPK/stratified abundance⭐⭐⭐
物种贡献解析哪些物种贡献了特定通路⭐⭐⭐

各步骤详解

第一步：KEGG Module系统与通路完整性概念

白话解释： KEGG Module是比KEGG Pathway更精细的功能单元——每个Module是完成一个特定代谢功能所需的最少基因集合。比如"M00001 Glycolysis"需要10个步骤的10个酶，如果一个微生物只有其中8个酶，则module completeness = 80%。这比简单统计KO数量更能反映真实的代谢能力。

技术细节： KEGG Module定义语法： - 每个module由若干"步骤"组成 - 每个步骤可以有多个替代酶（OR关系，用逗号分隔） - 步骤之间是AND关系（用空格分隔） - 模块完整性 = 完成的步骤数 / 总步骤数

# KEGG Module completeness 计算示例
# Module M00001 (Glycolysis, Embden-Meyerhof pathway)
# 定义：K00844,K12407,K00845 K01810 K00850,K16370 K01623,K01624 ...
# 表示：步骤1需要K00844或K12407或K00845，步骤2需要K01810，...

def calculate_module_completeness(module_definition, present_kos):
    """计算模块完整性"""
    steps = module_definition.split()
    completed = 0
    for step in steps:
        # 每个步骤中的alternatives用逗号分隔
        alternatives = step.split(',')
        if any(ko in present_kos for ko in alternatives):
            completed += 1
    return completed / len(steps)

# 示例
module_def = "K00844,K12407,K00845 K01810 K00850,K16370 K01623,K01624 K01803 K00134,K00150 K00927 K01834 K00873 K00163,K00161,K00162"
present_kos = {"K00845", "K01810", "K00850", "K01623", "K01803", "K00134", "K00927", "K01834", "K00873", "K00163"}
completeness = calculate_module_completeness(module_def, present_kos)
print(f"Glycolysis module completeness: {completeness:.0%}")  # 100%

第二步：KO注释——从序列到KEGG功能

白话解释： 第一步是给宏基因组中的基因分配KEGG Orthology（KO）编号。KO是KEGG对功能基因的分类系统。可以通过KofamScan（基于hmmer profile的最准确方法）或eggNOG-mapper进行注释。

技术细节：

# === 方法1：KofamScan（官方推荐，最准确）===
# 下载KOfam profiles和ko_list
wget ftp://ftp.genome.jp/pub/db/kofam/profiles.tar.gz
wget ftp://ftp.genome.jp/pub/db/kofam/ko_list.gz

# 运行KofamScan
exec_annotation \
    -f detail-tsv \
    --cpu 16 \
    -o kofamscan_results.txt \
    -p profiles/ \
    -k ko_list \
    predicted_proteins.faa

# === 方法2：HUMAnN3（自动化社区级功能注释）===
# HUMAnN3直接从reads进行通路重建
humann --input sample_R1.fq.gz \
    --output humann3_output/ \
    --threads 16 \
    --nucleotide-database chocophlan/ \
    --protein-database uniref/ \
    --metaphlan-options "--bowtie2db mpa_db"

# 输出文件：
# sample_genefamilies.tsv     - 基因家族丰度
# sample_pathabundance.tsv    - 通路丰度(MetaCyc)
# sample_pathcoverage.tsv     - 通路覆盖度

# === 方法3：从组装基因进行注释 ===
# 先用Prodigal预测基因
prodigal -i assembly.fa -a proteins.faa -p meta -f gff

# 再用eggNOG-mapper注释
emapper.py -i proteins.faa --output eggnog_results \
    -m diamond --cpu 16 \
    --dmnd_db eggnog_db/eggnog_proteins.dmnd

第三步：Module completeness评估

白话解释： 有了KO注释后，对每个KEGG Module计算其在样本中的完整性。完整性100%意味着社区中存在完成该代谢通路的所有必需酶（不一定在同一个物种中）。这反映了群落的整体代谢能力。

技术细节：

# === KEGG Module Completeness 批量计算 ===
import pandas as pd
from collections import defaultdict

# 读取KO注释结果
ko_annotations = pd.read_csv("kofamscan_results.txt", sep="\t")
# 过滤高置信度注释
ko_confident = ko_annotations[ko_annotations['score'] > ko_annotations['threshold']]
present_kos = set(ko_confident['KO'].unique())

# 读取KEGG Module定义
# 从KEGG API或本地文件获取
module_definitions = {}  # {module_id: definition_string}
with open("module_definitions.txt") as f:
    for line in f:
        parts = line.strip().split("\t")
        module_definitions[parts[0]] = parts[1]

# 批量计算completeness
results = {}
for module_id, definition in module_definitions.items():
    completeness = calculate_module_completeness(definition, present_kos)
    results[module_id] = completeness

# 排序输出
results_df = pd.DataFrame.from_dict(results, orient='index', columns=['completeness'])
results_df = results_df.sort_values('completeness', ascending=False)
complete_modules = results_df[results_df['completeness'] >= 0.75]
print(f"Modules with >=75% completeness: {len(complete_modules)}")

# === 使用MicrobiomeProfiler或KEGGModuleParser ===
library(KEGGREST)

# 从KEGG API获取module信息
module_info <- keggGet("M00001")

# 或使用专用R包
# BiocManager::install("MicrobiomeProfiler")
library(MicrobiomeProfiler)

# 计算多样本的module completeness
module_completeness <- calculate_module_score(
  ko_abundance_matrix,  # 行=KO, 列=样本
  module_db = kegg_modules
)

第四步：HUMAnN3通路重建与分层分析

白话解释： HUMAnN3不仅计算通路是否存在，还估计每条通路的丰度（RPK，reads per kilobase），并能告诉你每条通路中各物种的贡献比例。这种"分层"（stratified）分析是理解群落功能分工的关键。

技术细节：

# === HUMAnN3 分层分析 ===

# 合并多样本结果
humann_join_tables -i humann3_output/ \
    -o merged_pathabundance.tsv \
    --file_name pathabundance

# 归一化
humann_renorm_table -i merged_pathabundance.tsv \
    -o merged_pathabundance_cpm.tsv \
    -u cpm

# 分离分层和非分层结果
humann_split_stratified_table -i merged_pathabundance_cpm.tsv \
    -o stratified/
# 输出：stratified/merged_pathabundance_cpm_stratified.tsv
#        stratified/merged_pathabundance_cpm_unstratified.tsv

# 重组到KEGG KO
humann_regroup_table -i merged_genefamilies.tsv \
    -o merged_kegg_orthologs.tsv \
    -g uniref90_ko

# 注意：humann_regroup_table 只能映射到KO层级，不能直接输出KEGG Module
# 要获取Module级别丰度，需将KO丰度表结合KEGG Module定义自行计算
# 或使用第三方工具如KEGGDecoder/KEGG-Decoder进行Module汇总

# 分层结果解析——哪些物种贡献了特定通路
import pandas as pd

stratified = pd.read_csv("stratified/merged_pathabundance_cpm_stratified.tsv", sep="\t")

# 格式：PATHWAY|SPECIES  sample1  sample2 ...
# 如：PWY-6305: putrescine biosynthesis IV|g__Bacteroides.s__Bacteroides_uniformis

# 解析物种贡献
target_pathway = "LACTOSECAT-PWY"  # 乳糖降解通路
pathway_data = stratified[stratified.iloc[:, 0].str.startswith(target_pathway)]

# 分离物种名
pathway_data['species'] = pathway_data.iloc[:, 0].str.split('|').str[1]
species_contribution = pathway_data.groupby('species').mean()
print(species_contribution.sort_values(by=species_contribution.columns[0], ascending=False).head())

第五步：DRAM代谢注释（MAG水平）

白话解释： 如果已经有宏基因组装配的bins/MAGs（Metagenome-Assembled Genomes），DRAM可以对每个MAG做全面的代谢注释，包括碳/氮/硫代谢、能量代谢、有机物降解等，并生成漂亮的代谢功能热图。

技术细节：

# === DRAM (Distilled and Refined Annotation of Metabolism) ===

# 1. DRAM注释
DRAM.py annotate \
    -i 'mags/*.fa' \
    -o dram_annotation/ \
    --threads 32 \
    --min_contig_size 2500

# 2. DRAM蒸馏（生成代谢功能汇总）
DRAM.py distill \
    -i dram_annotation/annotations.tsv \
    -o dram_distill/ \
    --trna_path dram_annotation/trnas.tsv \
    --rrna_path dram_annotation/rrnas.tsv

# 输出：
# dram_distill/metabolism_summary.xlsx  - 代谢功能汇总表
# dram_distill/genome_stats.tsv         - 基因组统计
# dram_distill/product.html             - 交互式热图

# 3. DRAM生成KEGG module completeness
# annotations.tsv中包含KEGG KO信息，可用于计算每个MAG的module completeness

第六步：差异代谢通路分析

白话解释： 比较不同条件（如健康vs疾病）的代谢通路丰度差异——哪些通路在疾病中增强了（如致病相关代谢），哪些减弱了（如有益的短链脂肪酸产生）。

技术细节：

# === 差异通路分析 ===
library(Maaslin2)

# 使用MaAsLin2进行差异分析
fit <- Maaslin2(
  input_data = pathway_abundance,  # 通路×样本矩阵
  input_metadata = sample_metadata,
  output = "maaslin2_pathways",
  fixed_effects = c("group", "age", "sex"),
  normalization = "TSS",
  transform = "LOG",
  analysis_method = "LM"
)

# 或使用LEfSe/ANCOM-BC
# 或Wilcoxon检验 + FDR校正
library(coin)
pathway_results <- data.frame()
for (pw in rownames(pathway_abundance)) {
  pval <- wilcox.test(
    as.numeric(pathway_abundance[pw, metadata$group == "Disease"]),
    as.numeric(pathway_abundance[pw, metadata$group == "Control"])
  )$p.value
  fc <- mean(as.numeric(pathway_abundance[pw, metadata$group == "Disease"])) /
        mean(as.numeric(pathway_abundance[pw, metadata$group == "Control"]))
  pathway_results <- rbind(pathway_results,
                            data.frame(pathway = pw, p_value = pval, fold_change = fc))
}
pathway_results$padj <- p.adjust(pathway_results$p_value, method = "BH")

实战命令速查

# HUMAnN3完整流程
humann --input sample.fq.gz --output output/ --threads 16
humann_join_tables -i output/ -o merged.tsv --file_name pathabundance
humann_renorm_table -i merged.tsv -o merged_cpm.tsv -u cpm

# KofamScan
exec_annotation -f detail-tsv --cpu 16 -o results.txt -p profiles/ -k ko_list proteins.faa

# DRAM
DRAM.py annotate -i 'mags/*.fa' -o dram_out/ --threads 32
DRAM.py distill -i dram_out/annotations.tsv -o distill/

面试常问点

Q1: KEGG Module completeness与通路丰度有什么区别？

A: Module completeness是二元指标——评估通路是否"可能运作"（所有必需酶是否存在），不考虑丰度。通路丰度（如HUMAnN3输出的RPK）反映该通路的活跃程度（参与基因的reads数）。一个完整度100%但丰度低的通路可能只有少量微生物在维持；而丰度高说明多个微生物在大量表达相关酶。

Q2: HUMAnN3如何处理"社区级"通路重建？

A: HUMAnN3假设通路可以由不同物种的酶"拼凑"完成（community-level pathway）。先用MetaPhlAn做物种profiling，再将reads比对到物种特异性pangenome，然后用MinPath算法推断最小通路集。这允许一个通路的不同步骤由不同物种贡献——反映了微生物社区的代谢分工。

Q3: 为什么需要MinPath？直接统计KO不行吗？

A: 直接统计所有KO并映射到通路会严重高估活跃通路数（一个KO可属于多条通路，但可能只参与其中一条）。MinPath使用整数线性规划找到能解释所有观测KO的最小通路集，减少假阳性通路预测。这是一种保守但更准确的方法。

Q4: 宏基因组功能注释的局限性？

A: (1) DNA水平功能不等于活跃表达——宏转录组才能反映真实代谢活性；(2) 注释数据库不完善——约30-50%的宏基因组基因功能未知；(3) 社区级通路假设——不同物种的酶未必能在空间上协作；(4) 短reads的错误注释——短片段匹配可能不准确。

Q5: 如何从通路分析推断微生物群落的代谢互作？

A: (1) 分层分析确定各物种对通路的贡献；(2) 计算通路"互补性"——如果物种A有通路前半段酶，物种B有后半段，暗示代谢合作（syntrophy）；(3) 构建代谢网络图——连接有代谢物exchange的物种；(4) 结合宏转录组确认合作关系的活跃性。

易错点

1. 忽略通路的物种来源

社区级通路重建可能将不同物种的酶拼成一条"虚假"完整通路。应检查分层结果确认通路步骤是否来自能物理接触的物种。

2. 低丰度物种的通路假阳性

覆盖度很低的物种可能偶然比对到少量KO，被算入module completeness。应设置最小基因覆盖/表达阈值。

3. KEGG数据库更新问题

KEGG Module定义会随版本更新。确保分析使用一致的数据库版本，或标注所用版本便于复现。

4. 混淆基因存在与基因表达

宏基因组检测的是基因存在（DNA水平），不代表基因被表达或通路在运作。代谢活性分析需要宏转录组或宏蛋白组数据。

5. 不同注释工具结果差异大

KofamScan（hmm-based）和eggNOG-mapper（diamond-based）对同一基因可能给出不同KO。建议以KofamScan为准或取多工具共识。

补充知识

相关工具与数据库

  • KEGG Mapper：在线通路可视化
  • MetaCyc：HUMAnN3使用的通路数据库（比KEGG更细）
  • SEED/RAST：微生物基因组功能注释
  • antiSMASH：次级代谢产物生物合成基因簇预测
  • gutSMASH：肠道微生物特化代谢基因簇

引用推荐

  • HUMAnN3: Beghini et al., eLife, 2021
  • DRAM: Shaffer et al., Nucleic Acids Research, 2020
  • MinPath: Ye & Doak, PLoS Computational Biology, 2009
  • KofamScan: Aramaki et al., Bioinformatics, 2020

📚 相关文章推荐

宏基因组代谢通路重建：MinPath/Pathway Tools/KEGG Mapper
38% 相关
KEGG通路分析详解
35% 相关
181_微生物组功能通路差异
31% 相关
微生物组代谢预测（PICRUSt2 / Tax4Fun / MIMOSA）
29% 相关
宏基因组功能注释（HUMAnN3 / eggNOG-mapper）
28% 相关