cobrapy

COBRApy - 基于约束的重建与分析

概述

COBRApy 是一个用于代谢模型的基于约束的重建与分析（COBRA）的 Python 库，对系统生物学研究至关重要。使用全基因组规模代谢模型，执行细胞代谢的计算模拟，开展代谢工程分析，并预测表型行为。

**版本说明：**示例适用于 PyPI 上的 cobra 0.31.1（导入 cobra）。文档：cobrapy.readthedocs.io。代码库：opencobra/cobrapy。

何时使用此技能

在以下情况下使用此技能：

• 加载、构建或导出全基因组规模代谢模型（SBML、JSON、YAML）
• 对 COBRA 模型运行 FBA、pFBA、FVA 或通量采样
• 执行基因或反应敲除筛选以及产物包络分析
• 设计或优化生长培养基和交换约束
• 为不可行模型补充缺失反应，或验证模型一致性

安装

uv pip install "cobra==0.31.1"

MATLAB 模型 I/O（可选）：

uv pip install "cobra[array]==0.31.1"

COBRApy 使用 optlang 作为求解器接口。GLPK 会通过 swiglpk 自动安装。对于大型 MILP/QP，cobra 0.29+ 增加了混合求解器（HIGHS/OSQP）；现在将 model.solver = "osqp" 设置为通过 hybrid 路由，未来版本中在普通 LP 上可能会报错——若可用，请优先使用 model.solver = "hybrid"。

核心功能

COBRApy 提供全面的工具，分为几个关键领域：

1. 模型管理

从存储库或文件加载现有模型：

from cobra.io import load_model

# Bundled locally (no network): textbook, iJO1366, salmonella
model = load_model("textbook")      # alias for e_coli_core (95 reactions)
model = load_model("e_coli_core")   # same core E. coli model
model = load_model("iJO1366")       # genome-scale E. coli (bundled)
model = load_model("salmonella")    # Salmonella iYS1720 (bundled)

# Remote (BiGG / BioModels; requires network, cached after first fetch)
model = load_model("iML1515")       # E. coli genome-scale on BiGG

# Load from files
from cobra.io import read_sbml_model, load_json_model, load_yaml_model
model = read_sbml_model("path/to/model.xml")
model = load_json_model("path/to/model.json")
model = load_yaml_model("path/to/model.yml")

以各种格式保存模型：

from cobra.io import write_sbml_model, save_json_model, save_yaml_model
write_sbml_model(model, "output.xml")  # Preferred format
save_json_model(model, "output.json")  # For Escher compatibility
save_yaml_model(model, "output.yml")   # Human-readable

2. 模型结构和组件

访问并检查模型组件：

# Access components
model.reactions      # DictList of all reactions
model.metabolites    # DictList of all metabolites
model.genes          # DictList of all genes

# Get specific items by ID or index
reaction = model.reactions.get_by_id("PFK")
metabolite = model.metabolites[0]

# Inspect properties
print(reaction.reaction)        # Stoichiometric equation
print(reaction.bounds)          # Flux constraints
print(reaction.gene_reaction_rule)  # GPR logic
print(metabolite.formula)       # Chemical formula
print(metabolite.compartment)   # Cellular location

3. 通量平衡分析（FBA）

执行标准 FBA 模拟：

# Basic optimization
solution = model.optimize()
print(f"Objective value: {solution.objective_value}")
print(f"Status: {solution.status}")

# Access fluxes
print(solution.fluxes["PFK"])
print(solution.fluxes.head())

# Fast optimization (objective value only)
objective_value = model.slim_optimize()

# Change objective
model.objective = "ATPM"
solution = model.optimize()

简约 FBA（最小化总通量）：

from cobra.flux_analysis import pfba
solution = pfba(model)

几何 FBA（寻找中心解）：

from cobra.flux_analysis import geometric_fba
solution = geometric_fba(model)

4. 通量变异分析（FVA）

确定所有反应的通量范围：

from cobra.flux_analysis import flux_variability_analysis

# Standard FVA
fva_result = flux_variability_analysis(model)

# FVA at 90% optimality
fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)

# Loopless FVA (eliminates thermodynamically infeasible loops)
fva_result = flux_variability_analysis(model, loopless=True)

# FVA for specific reactions
fva_result = flux_variability_analysis(
    model,
    reaction_list=["PFK", "FBA", "PGI"]
)

5. 基因和反应删除研究

执行敲除分析：

from cobra.flux_analysis import (
    single_gene_deletion,
    single_reaction_deletion,
    double_gene_deletion,
    double_reaction_deletion
)

# Single deletions
gene_results = single_gene_deletion(model)
reaction_results = single_reaction_deletion(model)

# Double deletions (uses multiprocessing)
double_gene_results = double_gene_deletion(
    model,
    processes=4  # Number of CPU cores
)

# Manual knockout using context manager
with model:
    model.genes.get_by_id("b0008").knock_out()
    solution = model.optimize()
    print(f"Growth after knockout: {solution.objective_value}")
# Model automatically reverts after context exit

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

管理生长培养基：

# View current medium
print(model.medium)

# Modify medium (must reassign entire dict)
medium = model.medium
medium["EX_glc__D_e"] = 10.0  # Set glucose uptake
medium["EX_o2_e"] = 0.0       # Anaerobic conditions
model.medium = medium

# Calculate minimal media
from cobra.medium import minimal_medium

# Minimize total import flux
min_medium = minimal_medium(model, minimize_components=False)

# Minimize number of components (uses MILP, slower)
min_medium = minimal_medium(
    model,
    minimize_components=True,
    open_exchanges=True
)

7. 通量采样

对可行通量空间进行采样：

from cobra.sampling import sample

# Sample using OptGP (default, supports parallel processing)
samples = sample(model, n=1000, method="optgp", processes=4)

# Sample using ACHR
samples = sample(model, n=1000, method="achr")

# Validate samples
from cobra.sampling import OptGPSampler
sampler = OptGPSampler(model, processes=4)
sampler.sample(1000)
validation = sampler.validate(sampler.samples)
print(validation.value_counts())  # Should be all 'v' for valid

8. 产物包络

计算表型相平面：

from cobra.flux_analysis import production_envelope

# Standard production envelope
envelope = production_envelope(
    model,
    reactions=["EX_glc__D_e", "EX_o2_e"],
    objective="EX_ac_e"  # Acetate production
)

# With carbon yield
envelope = production_envelope(
    model,
    reactions=["EX_glc__D_e", "EX_o2_e"],
    carbon_sources="EX_glc__D_e"
)

# Visualize (use matplotlib or pandas plotting)
import matplotlib.pyplot as plt
envelope.plot(x="EX_glc__D_e", y="EX_o2_e", kind="scatter")
plt.show()

9. 缺失反应补充

添加反应以使模型可行：

from cobra.flux_analysis import gapfill

# Provide a universal reaction database (SBML/JSON); not bundled in cobra 0.31+
from cobra.io import read_sbml_model
universal = read_sbml_model("path/to/universal_reactions.xml")

# Perform gapfilling
with model:
    # Remove reactions to create gaps for demonstration
    model.remove_reactions([model.reactions.PGI])

    # Find reactions needed
    solution = gapfill(model, universal)
    print(f"Reactions to add: {solution}")

10. 模型构建

从头构建模型：

from cobra import Model, Reaction, Metabolite

# Create model
model = Model("my_model")

# Create metabolites
atp_c = Metabolite("atp_c", formula="C10H12N5O13P3",
                   name="ATP", compartment="c")
adp_c = Metabolite("adp_c", formula="C10H12N5O10P2",
                   name="ADP", compartment="c")
pi_c = Metabolite("pi_c", formula="HO4P",
                  name="Phosphate", compartment="c")

# Create reaction
reaction = Reaction("ATPASE")
reaction.name = "ATP hydrolysis"
reaction.subsystem = "Energy"
reaction.lower_bound = 0.0
reaction.upper_bound = 1000.0

# Add metabolites with stoichiometry
reaction.add_metabolites({
    atp_c: -1.0,
    adp_c: 1.0,
    pi_c: 1.0
})

# Add gene-reaction rule
reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"

# Add to model
model.add_reactions([reaction])

# Add boundary reactions
model.add_boundary(atp_c, type="exchange")
model.add_boundary(adp_c, type="demand")

# Set objective
model.objective = "ATPASE"

常见工作流

工作流 1：加载模型并预测生长

from cobra.io import load_model

# Load model (textbook = fast tutorial; iJO1366 / iML1515 for genome-scale)
model = load_model("textbook")

# Run FBA
solution = model.optimize()
print(f"Growth rate: {solution.objective_value:.3f} /h")

# Show active pathways
print(solution.fluxes[solution.fluxes.abs() > 1e-6])

工作流 2：基因敲除筛选

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

# Load model
model = load_model("textbook")
baseline = model.slim_optimize()

# Perform single gene deletions
results = single_gene_deletion(model)

# Find essential genes (growth < threshold)
essential_genes = results[results["growth"] < 0.01]
print(f"Found {len(essential_genes)} essential genes")

# Find genes with minimal impact
neutral_genes = results[results["growth"] > 0.9 * baseline]

工作流 3：培养基优化

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

# Load model
model = load_model("textbook")

# Calculate minimal medium for 50% of max growth
target_growth = model.slim_optimize() * 0.5
min_medium = minimal_medium(
    model,
    target_growth,
    minimize_components=True
)

print(f"Minimal medium components: {len(min_medium)}")
print(min_medium)

工作流 4：通量不确定性分析

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

# Load model
model = load_model("textbook")

# First check flux ranges at optimality
fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# For reactions with large ranges, sample to understand distribution
samples = sample(model, n=1000)

# Analyze specific reaction
reaction_id = "PFK"
import matplotlib.pyplot as plt
samples[reaction_id].hist(bins=50)
plt.xlabel(f"Flux through {reaction_id}")
plt.ylabel("Frequency")
plt.show()

工作流 5：用于临时更改的上下文管理器

使用上下文管理器进行临时修改：

# Model remains unchanged outside context
with model:
    # Temporarily change objective
    model.objective = "ATPM"

    # Temporarily modify bounds
    model.reactions.EX_glc__D_e.lower_bound = -5.0

    # Temporarily knock out genes
    model.genes.b0008.knock_out()

    # Optimize with changes
    solution = model.optimize()
    print(f"Modified growth: {solution.objective_value}")

# All changes automatically reverted
solution = model.optimize()
print(f"Original growth: {solution.objective_value}")

关键概念

DictList 访问模式、通量边界约定、基因-反应规则（GPR）以及 EX_ 交换反应符号约定，详见 references/api_quick_reference.md 中的“Key Concepts”。

最佳实践

1. 使用上下文管理器进行临时修改，避免状态管理问题
2. 验证模型：分析前使用 model.slim_optimize() 确保模型可行
3. 检查求解状态：优化后检查状态；optimal 表示成功求解
4. 在热力学可行性很重要时使用无环 FVA
5. 在 FVA 中适当设置 fraction_of_optimum，以探索次优空间
6. 并行化计算开销较大的操作（采样、双重删除）——在全基因组规模模型上，从较小的 n 和 processes=1 开始
7. 模型交换和长期存储优先使用 SBML 格式
8. 仅需目标值时使用 slim_optimize()，以提升性能
9. 验证通量样本，确保数值稳定性
10. 确认输出路径，再从工作流示例中写入 CSV/PNG 文件

故障排查

不可行解：检查培养基约束、反应边界和模型一致性 优化速度慢：通过 model.solver 尝试不同的求解器（GLPK、CPLEX、Gurobi） 无界解：确认交换反应设置了适当的上界 导入错误：确保文件格式正确，且 SBML 标识符有效

参考资料

详细的工作流和 API 模式请参阅：

• references/workflows.md - 全面的分步工作流示例
• references/api_quick_reference.md - 常用函数签名和模式

官方文档：https://cobrapy.readthedocs.io/en/latest/

Agent 操作流程

1. 检查环境。 确认工具版本、参考基因组／注释版本以及输入格式（FASTQ、BAM、VCF、h5ad）。
2. 明确输入。 根据数据或用户确认格式、标识符和参数。任何会改变结果的值都应询问，不要猜测。
3. 先运行小规模版本。 先在小型子集上运行流水线（一个样本、一条染色体、几千个细胞）。
4. 按照上述说明和参考资料执行完整任务。
5. 验证结果。 在解读结果前，检查 QC 指标、样本身份、基因组版本一致性和批次效应。
6. 报告。 说明运行了什么（版本、命令、参数）、检查了什么，以及仍有哪些不确定之处。

如果发生这种情况	请执行以下操作
输入之间的基因组版本或标识符不匹配	停止并先进行协调（liftover、ID 映射），再继续。
这些说明中的函数、标志或端点在已安装版本中不存在	查阅已安装版本自己的文档（help()、--help、官方文档），进行调整并告知用户。切勿编造 API。
必需的输入、标识符或参数不明确	询问用户，或在运行前明确说明所采用的假设。

完整性规则

• 切勿捏造结果、参数、标识符、引用或统计数据。如果无法运行或验证某项内容，请直说。
• 不要解读统计数据所不支持的生物学意义；报告多重检验校正情况。
• 此处与版本相关的细节可能已经过时：请对照已安装版本的官方文档进行确认。
• 对于会产生费用、占用共享 GPU 或云配额、涉及个人或患者数据，或无法撤销的操作，请先征求同意。

相关技能

• 13c-metabolic-flux：使用经过验证的原子映射和 mfapy 同位素……，根据稳态碳 13 同位素示踪测量结果估算细胞内代谢通量。