732. 宏基因组代谢建模GEM（Genome-scale Metabolic Models）¶

数学表示：
maximize: c^T * v        # 目标：最大化生物量
subject to: S * v = 0    # 约束：代谢物浓度不变（稳态）
            lb ≤ v ≤ ub  # 约束：反应速率有上下限

# 克隆metaGEM仓库
git clone https://github.com/franciscozorrilla/metaGEM.git  # 下载metaGEM源代码

# 进入目录
cd metaGEM  # 切换到metaGEM目录

# 使用conda创建环境（metaGEM自带环境配置文件）
conda env create -f envs/metagem.yaml  # 从yaml文件创建conda环境

# 激活环境
conda activate metagem  # 激活metaGEM专用环境

# 验证安装
snakemake --version  # 检查Snakemake是否正确安装

原始reads → 质控(fastp) → 组装(MEGAHIT) → 分箱(MetaBAT2/MaxBin2) 
→ 优化(MetaWRAP) → 质检(CheckM) → 分类(GTDB-Tk) 
→ 建模(CarveMe) → 质评(MEMOTE) → 互作(SMETANA)

# ===== 第一步：用CarveMe从MAG构建GEM =====
# CarveMe是自动化GEM重建工具，输入基因组蛋白序列，输出SBML格式模型

# 安装CarveMe
pip install carveme  # 安装CarveMe（需要CPLEX或GLPK求解器）

# 预测蛋白序列（如果没有的话）
prodigal -i mag_bin1.fa \  # 输入：MAG的fasta文件
         -a mag_bin1.faa \  # 输出：蛋白序列（氨基酸fasta）
         -o mag_bin1.gff \  # 输出：基因预测结果GFF格式
         -p meta  # 使用宏基因组模式（更适合不完整基因组）

# 用CarveMe构建GEM
carve mag_bin1.faa \  # 输入：蛋白序列文件
     -o mag_bin1_gem.xml \  # 输出：SBML格式的代谢模型
     --solver glpk  # 指定求解器（免费的用glpk，商业的用cplex更快）

# 用MEMOTE测试GEM质量
pip install memote  # 安装MEMOTE

# 运行质量测试
memote run mag_bin1_gem.xml \  # 输入：GEM模型文件
       --filename report.html  # 输出：HTML格式的质量报告

# 报告包含：
# - 模型一致性（反应平衡、代谢物完整性）
# - 注释覆盖率（基因、反应、代谢物是否有数据库ID）
# - 生物量反应检查
# - 总体得分（满分100）

# ===== 用COBRApy进行通量平衡分析 =====
import cobra  # 导入COBRA工具包

# 加载GEM模型
model = cobra.io.read_sbml_model("mag_bin1_gem.xml")  # 从SBML文件读取模型

# 查看模型基本信息
print(f"反应数: {len(model.reactions)}")  # 打印模型中代谢反应的总数
print(f"代谢物数: {len(model.metabolites)}")  # 打印代谢物总数
print(f"基因数: {len(model.genes)}")  # 打印关联基因数

# 设置培养基条件（模拟肠道环境）
# 通过设置exchange reaction的下界来定义可用营养物
for rxn in model.exchanges:  # 遍历所有交换反应（与外界交换物质的反应）
    rxn.lower_bound = 0  # 先关闭所有营养物输入

# 开放常见肠道营养物
model.reactions.get_by_id("EX_glc__D_e").lower_bound = -10  # 葡萄糖摄入上限
model.reactions.get_by_id("EX_o2_e").lower_bound = -2  # 氧气（厌氧环境，限制氧气）
model.reactions.get_by_id("EX_nh4_e").lower_bound = -10  # 铵盐（氮源）

# 运行FBA
solution = model.optimize()  # 求解最优代谢方案
print(f"预测生长速率: {solution.objective_value:.4f} /h")  # 打印最大生物量产生速率

# 查看关键反应通量
print(f"葡萄糖消耗: {solution.fluxes['EX_glc__D_e']:.4f}")  # 葡萄糖的消耗量
print(f"乳酸产生: {solution.fluxes.get('EX_lac__D_e', 0):.4f}")  # 乳酸的分泌量

# ===== 安装SMETANA =====
pip install smetana  # 安装SMETANA群落互作分析工具

# 运行群落代谢互作分析
smetana models/*.xml \  # 输入：多个物种的GEM模型文件
       -o community_interactions.tsv \  # 输出：互作预测结果表格
       --flavor bigg \  # 使用BiGG数据库的命名规范
       --solver glpk \  # 指定求解器
       -v  # 显示详细运行信息

# SMETANA输出关键指标：
# MIP (Metabolic Interaction Potential): 群落互补性，值越大说明互依赖越强
# MRO (Metabolic Resource Overlap): 资源竞争程度，值越大说明竞争越强
# SCS (Species Coupling Score): 物种间耦合度

# ===== 用Python分析SMETANA结果 =====
import pandas as pd  # 导入pandas用于数据处理
import matplotlib.pyplot as plt  # 导入绑图库
import seaborn as sns  # 导入seaborn增强绑图

# 读取SMETANA结果
df = pd.read_csv("community_interactions.tsv", sep="\t")  # 读取互作结果表

# 筛选高置信度互作
significant = df[df["score"] > 0.5]  # 只保留得分>0.5的显著互作

# 绘制互作热图
pivot = significant.pivot_table(
    index="receiver",  # 行：接收物质的物种
    columns="donor",   # 列：提供物质的物种
    values="score",    # 值：互作强度
    aggfunc="sum"      # 聚合方式：求和
)

plt.figure(figsize=(12, 10))  # 设置画布大小
sns.heatmap(pivot, cmap="YlOrRd", annot=True, fmt=".2f")  # 绘制热图
plt.title("Species Metabolic Interaction Network")  # 设置标题
plt.tight_layout()  # 自动调整布局
plt.savefig("metabolic_interactions_heatmap.png", dpi=300)  # 保存高分辨率图片
plt.show()  # 显示图片

# ===== IMIC流程示意 =====
# IMIC的核心改进：用转录组数据约束代谢模型

# 1. 传统方法：只用基因组 → 知道"能做什么"但不知道"实际在做什么"
# 2. IMIC方法：基因组+转录组 → 同时知道"能做什么"和"正在做什么"

# 伪代码示意IMIC的工作原理
"""
输入：
  - MAG构建的GEM模型
  - 宏转录组（metatranscriptome）数据

步骤：
  1. 将转录组reads比对到MAG的基因上
  2. 计算每个基因的表达量（TPM）
  3. 用表达量作为约束，限制对应反应的通量范围
     - 高表达基因 → 对应反应的通量上限提高
     - 低表达基因 → 对应反应的通量上限降低
  4. 在新约束下重新运行FBA

优势：
  - 更准确预测特定条件下的代谢状态
  - 可以比较不同环境/疾病状态下的代谢差异
"""

# ===== GEM分析速查 =====

# 1. MAG → 蛋白
prodigal -i mag.fa -a mag.faa -p meta

# 2. 蛋白 → GEM
carve mag.faa -o mag_gem.xml --solver glpk

# 3. GEM质检
memote run mag_gem.xml --filename report.html

# 4. FBA模拟
python -c "
import cobra
m = cobra.io.read_sbml_model('mag_gem.xml')
print(m.optimize().objective_value)
"

# 5. 群落互作
smetana models/*.xml -o interactions.tsv --solver glpk

# ===== 关键Python包 =====
pip install cobra carveme memote smetana

# ===== 关键数据库 =====
# BiGG Models: http://bigg.ucsd.edu/
# ModelSEED: https://modelseed.org/
# VMH: https://www.vmh.life/ (人类代谢+微生物模型)