cobrapy-metabolic-modeling

COBRApy — 基于约束的代谢建模

概述

COBRApy 是一个 Python 包，用于对全基因组规模的代谢模型进行基于约束的重建与分析（COBRA）。它支持在 SBML 格式的代谢网络上进行通量平衡分析（FBA）、通量变异性分析（FVA）、基因与反应敲除筛选、通量采样、产物包络分析、补缺以及培养基优化。

适用场景

• 预测不同营养条件下的微生物生长速率（FBA）
• 通过单基因和双基因敲除筛选，识别必需基因或反应
• 确定通量范围和其他最优解（FVA）
• 对可行通量分布进行采样，以刻画代谢灵活性
• 设计最小生长培养基或优化碳源
• 计算代谢工程目标的产物包络
• 使用通用反应数据库为不完整的草稿模型补缺
• 对于动力学建模或基于动态 ODE 的模型，请改用 Tellurium
• 对于代谢图上的通路可视化，请改用 Escher

前置条件

• Python 包：cobra（包含 GLPK 求解器）、numpy、pandas
• 可选求解器：CPLEX 或 Gurobi（处理大型模型时速度更快，但需要许可证）
• 数据：SBML（.xml）、JSON 或 YAML 格式的代谢模型文件；可从 BiGG Models、AGORA 或 ModelSEED 获取

pip install cobra

预先访谈

编写任何分析代码之前，先与用户确认以下事项。

decisions:
  - id: D1
    param: modelSource
    kind: required
    source: user
    ask: "Which genome-scale reconstruction, and which version, should be used?"
    default: null

  - id: D2
    param: objectiveFunction
    kind: required
    source: user
    depends_on: [D1]
    ask: "What should the model be asked to maximize - biomass production, a specific product, or ATP yield?"
    default: "the reconstruction's own biomass reaction"

  - id: D3
    param: mediumConstraints
    kind: required
    source: user
    depends_on: [D1]
    ask: "Which nutrients are available, and at what uptake rates?"
    default: "the reconstruction's default medium"

  - id: D4
    param: analysisType
    kind: required
    source: user
    ask: "A single optimal flux distribution, the feasible range of each flux, or a sampled distribution over the solution space?"
    default: "single optimal distribution"

  - id: D5
    param: optimalityFraction
    kind: required
    source: user
    depends_on: [D4]
    ask: "What fraction of the optimum must solutions retain when reporting flux ranges?"
    default: 1.0
    skip_if: "not a flux-variability analysis"

  - id: D6
    param: looplessConstraint
    kind: optional
    source: user
    ask: "Should thermodynamically infeasible internal loops be excluded, at a large runtime cost?"
    default: "loops allowed"

  - id: D7
    param: sampleCount
    kind: optional_conditional
    source: user
    depends_on: [D4]
    ask: "How many flux samples should be drawn, and how far apart?"
    default: "not sampling"
    skip_if: "not a sampling analysis"

  - id: D8
    param: solverProcesses
    kind: never_ask
    source: data
    reason: "Parallel workers affect runtime only, not the fluxes"
    default: "min(4, available_cores)"

D2 和 D3 合起来才构成模型：通量平衡分析会根据给定的目标函数和培养基报告最优值；无论这两者是否正确，它都会同样给出完整且可行的通量分布。默认的生物量目标函数包含了原作者对生长的假设，但这些假设未必适用于当前研究条件。

快速入门

from cobra.io import load_model

model = load_model("textbook")  # E. coli core model
print(f"Model: {model.id} — {len(model.reactions)} rxns, {len(model.metabolites)} mets, {len(model.genes)} genes")

solution = model.optimize()
print(f"Growth rate: {solution.objective_value:.4f} /h")
print(f"Status: {solution.status}")
# Model: e_coli_core — 95 rxns, 72 mets, 137 genes
# Growth rate: 0.8739 /h
# Status: optimal

核心 API

1. 模型 I/O

加载内置模型，并读取或写入标准格式。

from cobra.io import load_model, read_sbml_model, write_sbml_model, load_json_model, save_json_model

# Bundled: "textbook" (95 rxns), "ecoli" (2583 rxns), "salmonella"
model = load_model("textbook")
# model = read_sbml_model("my_model.xml")   # from SBML file
# model = load_json_model("my_model.json")  # from JSON file

write_sbml_model(model, "output_model.xml")
save_json_model(model, "output_model.json")
print(f"Saved model: {model.id}")

2. 模型结构与组件

通过 DictList 容器访问反应、代谢物和基因。

from cobra.io import load_model
model = load_model("textbook")

# Inspect a reaction
rxn = model.reactions.get_by_id("PFK")
print(f"Reaction: {rxn.id} — {rxn.name}")
print(f"Equation: {rxn.reaction}")
print(f"Bounds: {rxn.bounds}, GPR: {rxn.gene_reaction_rule}")

# Inspect a metabolite
met = model.metabolites.get_by_id("atp_c")
print(f"Metabolite: {met.id}, Formula: {met.formula}, Compartment: {met.compartment}")

# Query and list exchange reactions
atp_rxns = model.reactions.query("atp", attribute="name")
print(f"ATP-related reactions: {len(atp_rxns)}, Exchange reactions: {len(model.exchanges)}")

3. 通量平衡分析（FBA）

通过最大化目标函数来预测最优通量分布。

from cobra.io import load_model
from cobra.flux_analysis import pfba

model = load_model("textbook")

# Standard FBA
solution = model.optimize()
print(f"Growth: {solution.objective_value:.4f} /h, Active fluxes: {(solution.fluxes.abs() > 1e-6).sum()}")

# Parsimonious FBA — same growth, minimal total flux
pfba_sol = pfba(model)
print(f"pFBA total flux: {pfba_sol.fluxes.abs().sum():.1f} vs standard: {solution.fluxes.abs().sum():.1f}")

# Change objective; slim_optimize for speed
from cobra.io import load_model
model = load_model("textbook")

with model:
    model.objective = "ATPM"
    print(f"Max ATPM flux: {model.optimize().objective_value:.2f}")

print(f"Growth (slim): {model.slim_optimize():.4f}")  # no flux vector, faster

4. 通量变异性分析（FVA）

确定最优或接近最优时的可行通量范围。

from cobra.io import load_model
from cobra.flux_analysis import flux_variability_analysis

model = load_model("textbook")

fva = flux_variability_analysis(model, fraction_of_optimum=1.0)
fva_90 = flux_variability_analysis(model, fraction_of_optimum=0.9)
fva["range"] = fva["maximum"] - fva["minimum"]
fva_90["range"] = fva_90["maximum"] - fva_90["minimum"]
print(f"Mean range at 100%: {fva['range'].mean():.2f}, at 90%: {fva_90['range'].mean():.2f}")

# Loopless FVA on specific reactions
from cobra.io import load_model
from cobra.flux_analysis import flux_variability_analysis

model = load_model("textbook")
fva_ll = flux_variability_analysis(
    model, loopless=True, reaction_list=["PFK", "PGI", "FBA", "TPI", "GAPD"],
)
print(fva_ll)

5. 基因与反应敲除

通过敲除模拟筛选必需基因或反应。

from cobra.io import load_model
from cobra.flux_analysis import single_gene_deletion, double_gene_deletion

model = load_model("textbook")
wt_growth = model.slim_optimize()

# Single gene deletions
gene_results = single_gene_deletion(model)
gene_results["growth_fraction"] = gene_results["growth"] / wt_growth
essential = gene_results[gene_results["growth_fraction"] < 0.01]
print(f"Essential genes: {len(essential)} / {len(model.genes)}")

# Double deletions (synthetic lethality) — use multiprocessing
double_results = double_gene_deletion(model, processes=4)
print(f"Double deletion results: {double_results.shape}")

6. 生长培养基与最小培养基

调整营养物质的可用性并计算最小培养基。

from cobra.io import load_model
from cobra.medium import minimal_medium

model = load_model("textbook")

# View current medium
for rxn_id, flux in sorted(model.medium.items()):
    print(f"  {rxn_id}: {flux}")

# Anaerobic switch via context manager
with model:
    medium = model.medium
    medium["EX_o2_e"] = 0.0
    model.medium = medium
    print(f"Anaerobic growth: {model.slim_optimize():.4f} /h")

# Minimal medium
min_med = minimal_medium(model, minimize_components=True, open_exchanges=True)
print(f"Minimal medium: {len(min_med)} components")

7. 通量采样

对可行通量分布进行采样，用于变异性分析。

from cobra.io import load_model
from cobra.sampling import sample

model = load_model("textbook")
samples = sample(model, n=500, method="optgp")
print(f"Samples shape: {samples.shape}")  # (500, n_reactions)
print(f"PFK flux: mean={samples['PFK'].mean():.2f}, std={samples['PFK'].std():.2f}")

8. 产物包络与补缺

计算表型相平面并填补模型缺口。

from cobra.io import load_model
from cobra.flux_analysis import production_envelope

model = load_model("textbook")
envelope = production_envelope(
    model,
    reactions=model.reactions.get_by_id("EX_ac_e"),
    carbon_sources=model.reactions.get_by_id("EX_glc__D_e"),
)
print(f"Envelope: {len(envelope)} points")
print(envelope[["flux_minimum", "flux_maximum", "carbon_yield_minimum", "carbon_yield_maximum"]].head())

# Gapfilling: restore growth after reaction removal
from cobra.io import load_model
from cobra.flux_analysis.gapfilling import gapfill

model = load_model("textbook")
universal = load_model("textbook")  # In practice, use a universal reaction DB
with model:
    model.remove_reactions([model.reactions.get_by_id("PFK")])
    print(f"Growth after removing PFK: {model.slim_optimize():.4f}")
    for rxn in gapfill(model, universal)[0]:
        print(f"  Gapfill suggests: {rxn.id}")

核心概念

DictList 对象

反应、代谢物和基因存储在 DictList 中；它支持排序、索引访问，也可通过 ID 访问。

rxn = model.reactions[0]                    # by index
rxn = model.reactions.get_by_id("PFK")      # by ID
matches = model.reactions.query("phospho")  # keyword search

交换反应

带有 EX_ 前缀的反应表示系统边界。正通量 = 分泌；负通量 = 摄取。通过 model.medium 字典管理。

基因-反应规则（GPR）

连接基因与反应的布尔表达式：(b0726 and b0727) or b1234。基因敲除会根据 GPR 逻辑传播。

上下文管理器

with model: 会创建状态快照，并在退出时还原所有更改（边界、目标函数、培养基、敲除）。支持嵌套。

常见工作流

工作流 1：基因敲除筛选

目标：识别必需基因、降低生长的基因和中性基因。

from cobra.io import load_model
from cobra.flux_analysis import single_gene_deletion

model = load_model("textbook")
wt_growth = model.slim_optimize()

results = single_gene_deletion(model)
results["growth_fraction"] = results["growth"] / wt_growth

essential = results[results["growth_fraction"] < 0.01]
reduced = results[(results["growth_fraction"] >= 0.01) & (results["growth_fraction"] < 0.9)]
neutral = results[results["growth_fraction"] >= 0.9]

print(f"Essential: {len(essential)}, Reduced: {len(reduced)}, Neutral: {len(neutral)}")
for idx in essential.index:
    print(f"  Essential gene: {list(idx)[0]}")

工作流 2：培养基优化

目标：寻找不同生长目标下的最小培养基；比较有氧与无氧条件。

from cobra.io import load_model
from cobra.medium import minimal_medium
import pandas as pd

model = load_model("textbook")
results = []
for frac in [0.1, 0.5, 0.8, 1.0]:
    with model:
        target = model.slim_optimize() * frac
        model.reactions.get_by_id("Biomass_Ecoli_core").lower_bound = target
        try:
            mm = minimal_medium(model, minimize_components=True, open_exchanges=True)
            results.append({"growth_frac": frac, "n_components": len(mm)})
        except Exception:
            results.append({"growth_frac": frac, "n_components": None})
print(pd.DataFrame(results).to_string(index=False))

for label, o2 in [("Aerobic", 1000.0), ("Anaerobic", 0.0)]:
    with model:
        medium = model.medium
        medium["EX_o2_e"] = o2
        model.medium = medium
        print(f"{label} growth: {model.slim_optimize():.4f} /h")

工作流 3：生产菌株设计

目标：设计能够最大限度生产目标代谢物的菌株。结合模块 3、5 和 8。

1. 加载模型并计算目标代谢物的野生型生产包络
2. 筛选单基因敲除，找出高产菌株（筛选目标分泌量增加的情况）
3. 组合排名靠前的敲除，并通过生产条件下的 FBA 进行验证
4. 如果敲除导致不可行，则进行补缺
5. 计算最终生产包络，并与野生型进行比较

关键参数

参数	模块/函数	默认值	范围 / 选项	影响
fraction_of_optimum	flux_variability_analysis	1.0	0.0-1.0	需要维持的最大目标值比例；值越低，通量范围越宽
loopless	flux_variability_analysis	False	True、False	消除热力学上不可行的循环；速度较慢
method	sample	"optgp"	"optgp"、"achr"	采样算法；optgp 支持并行处理
n	sample	必填	100-10000	要抽取的通量样本数量
processes	sample、double_gene_deletion	1	1-N_cores	并行工作进程数
minimize_components	minimal_medium	False	True、False	True = 营养物质数量最少（MILP）；False = 总通量最小
open_exchanges	minimal_medium	False	True、False	将所有交换反应都作为候选营养物质
carbon_sources	production_envelope	None	Reaction 对象	在计算通量包络的同时计算碳产率
thinning	sample	100	1-1000	保留样本之间的步数；值越高，相关性越低

最佳实践

1. 对临时更改使用上下文管理器：退出 with model: 后会还原所有修改。

   with model:
       model.reactions.PFK.knock_out()
       print(model.slim_optimize())  # modified
   # model is restored here
   

2. 分析前使用 slim_optimize() 验证：在执行耗时操作（FVA、采样）前快速检查可行性。

3. 优化后检查 solution.status：在解读通量前，务必确认状态为 "optimal"。

4. 需要考虑热力学可行性时使用无循环 FVA：标准 FVA 可能包含不可行的内部循环，从而扩大通量范围。

5. 并行执行耗时操作：采样和双基因敲除支持 processes 参数。

6. 模型交换时优先使用 SBML：这是 COBRA 工具都支持的社区标准。

7. 在循环中使用 slim_optimize()：跳过完整通量向量的构建，在筛选时速度显著更快。

8. 验证通量样本：使用 sampler.validate(samples) 检查化学计量和边界约束。

常用操作方案

操作方案：参数扫描（葡萄糖摄取速率）

from cobra.io import load_model

model = load_model("textbook")
print("Glucose_uptake | Growth_rate")
for glc_uptake in [1, 2, 5, 10, 15, 20]:
    with model:
        model.reactions.get_by_id("EX_glc__D_e").lower_bound = -glc_uptake
        growth = model.slim_optimize()
        print(f"  {glc_uptake:>13} | {growth:.4f}")

操作方案：批次条件分析

from cobra.io import load_model
import pandas as pd

model = load_model("textbook")
conditions = [
    {"name": "Rich aerobic", "EX_o2_e": 1000, "EX_glc__D_e": 10},
    {"name": "Anaerobic", "EX_o2_e": 0, "EX_glc__D_e": 10},
    {"name": "Low glucose", "EX_o2_e": 1000, "EX_glc__D_e": 1},
]
results = []
for c in conditions:
    with model:
        medium = model.medium
        medium["EX_o2_e"], medium["EX_glc__D_e"] = c["EX_o2_e"], c["EX_glc__D_e"]
        model.medium = medium
        results.append({"condition": c["name"], "growth": round(model.slim_optimize(), 4)})
print(pd.DataFrame(results).to_string(index=False))

操作方案：模型验证清单

from cobra.io import load_model
from cobra.flux_analysis import find_blocked_reactions

model = load_model("textbook")

# Feasibility, mass balance, dead-ends, blocked reactions
print(f"Growth feasible: {model.slim_optimize() > 0}")
print(f"Missing formula: {sum(1 for m in model.metabolites if m.formula is None)}")
print(f"Dead-end metabolites: {sum(1 for m in model.metabolites if len(m.reactions) == 1)}")
print(f"Blocked reactions: {len(find_blocked_reactions(model))} / {len(model.reactions)}")

故障排除

问题	原因	解决方案
solution.status == "infeasible"	约束条件无法同时满足	检查培养基是否含有所需营养物质；验证反应边界；使用 model.medium 恢复默认值
solution.status == "unbounded"	通量没有上限	为交换反应设置有限的上界
优化非常缓慢	模型较大且使用默认的 GLPK 求解器	安装 CPLEX 或 Gurobi：model.solver = "cplex"
设置边界时出现 ValueError	暂时出现 lower_bound > upper_bound	将边界设置为元组：rxn.bounds = (new_lb, new_ub)
基因敲除返回 NaN	敲除导致模型不可行	对必需基因而言这是预期情况；将其归类为必需基因
读取 SBML 时出现 IOError	SBML 无效或缺少命名空间	在 sbml.org 上验证；尝试 cobra.io.sbml.validate_sbml_model(path)
通量样本验证失败	数值求解器容差	增大 thinning 参数；尝试 method="achr"

捆绑资源

1 个参考文件：

• references/api_workflows.md — 汇总 API 速查和高级工作流。涵盖：详细的函数签名、求解器配置（GLPK/CPLEX/Gurobi）、高级分析（find_blocked_reactions、find_essential_genes/find_essential_reactions）、模型操作（添加反应/代谢物/基因）、通量样本验证，以及 5 个工作流示例（敲除与可视化、培养基设计、通量空间探索、生产菌株设计、模型验证）。基础 FBA/FVA/删除/采样内容已移至文中（核心 API 模块 3-7）。未涵盖：几何 FBA 内部机制、MIP 间隙配置——请参阅 COBRApy 文档。

相关技能

• escher — 交互式代谢图可视化；使用 COBRApy 中的 JSON 模型
• bioservices — 从 BiGG、BioModels 或 KEGG 获取模型，以导入 COBRApy
• networkx — 代谢网络拓扑分析；将化学计量矩阵导出为图

参考资料

• COBRApy 文档 — 官方 API 参考和教程
• BiGG Models — 精选的基因组规模代谢模型库
• COBRApy GitHub — 源代码和问题跟踪器
• Ebrahim 等人（2013）《COBRApy：基于约束的 Python 重建与分析》— BMC Systems Biology