alterlab-cobrapy

COBRApy - 制約ベースの再構築と解析

概要

COBRApyは、システム生物学の研究に不可欠な、代謝モデルの制約ベース再構築と解析（COBRA）を行うためのPythonライブラリです。ゲノムスケール代謝モデルを扱い、細胞代謝の計算シミュレーションを実行し、代謝工学解析を行い、表現型の挙動を予測できます。

このスキルを使用する場面

ユーザーが次のことを求めている場合に、このスキルを使用します。

• ゲノムスケール代謝モデル（SBML/JSON/YAML）を読み込み、FBA/pFBAで増殖をシミュレーションする。
• ノックアウトの表現型、必須遺伝子、または最小培地を予測する。
• FVA、フラックスサンプリング、または生産エンベロープを使ってフラックス空間を調べる。
• 制約ベースモデルを構築、ギャップフィリング、またはキュレーションする。

対象外

シナリオ	代わりに使用
シミュレーションではなく、経路や反応のアノテーションを調べる	alterlab-kegg、alterlab-reactome
反応の酵素動力学パラメーター（kcat、Km）	alterlab-brenda
フラックスではなく、代謝ネットワークのグラフ解析（中心性、コミュニティ）	alterlab-networkx
メタボロミクス測定値とその同定	alterlab-metabolomics-wb、alterlab-hmdb
フラックスバランス問題ではない多目的最適化	alterlab-pymoo

インストールと要件

uv pip install cobra              # 0.32.x as of 2026-09
uv pip install 'cobra[array]'     # scipy extra, for array-based model handling
uv pip install 'cobra[chrr]'      # hopsy, enables the fast CHRR sampler

要件： optlangとソルバーに依存します。GLPKはswiglpk経由でデフォルトで同梱されます。一方、CPLEXとGurobiはオプションで、ゲノムスケールモデルでは大幅に高速です。pandasは必須依存関係です（< 3に固定）。matplotlib/seabornはプロットにのみ必要です。フラックスサンプリングと並列デリーションではmultiprocessingを使用するため、Windows/macOSのspawnではスクリプトをif __name__ == "__main__":で保護してください。

同梱モデル名。 load_modelが認識する同梱エイリアスは3つのみです。"textbook"（E. coliのコアモデル、e_coli_core）、"iJO1366"（E. coliの完全なゲノムスケールモデル）、および"salmonella"です。それ以外の識別子（例："ecoli"、"e_coli_core"）は同梱エイリアスではなく、BiGG/BioModelsへのリモート検索が行われます。これにはネットワークアクセスが必要で、失敗する場合があります。独自のモデルにはread_sbml_model("path.xml")を使用してください。

主な機能

COBRApyは、いくつかの主要な領域に整理された包括的なツールを提供します。

1. モデル管理

リポジトリまたはファイルから既存のモデルを読み込みます。

from cobra.io import load_model

# Load bundled test models
model = load_model("textbook")  # E. coli core model
model = load_model("iJO1366")   # Full E. coli genome-scale model (BiGG)
model = load_model("salmonella")

# 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

# method defaults to "auto" in current cobrapy; the three options are:
#   "chrr"  - coordinate hit-and-run with rounding, via the C++ hopsy library
#             (install with cobra[chrr]); usually the fastest and best-mixing
#   "optgp" - parallel, wants large n (>1000) to be efficient
#   "achr"  - artificial centering hit-and-run, better for small n
samples = sample(model, n=1000)                              # auto
samples = sample(model, n=1000, method="optgp", processes=4)  # explicit OptGP
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. ギャップフィリング

モデルを実行可能にするために反応を追加します。

import cobra
from cobra.flux_analysis import gapfill

# Build a universal model of candidate reactions to draw from.
# In practice, load a curated reaction database (e.g. the BiGG
# universal reactions JSON via load_json_model) or assemble reactions
# into an empty cobra.Model — there is no load_model("universal").
universal = cobra.Model("universal_reactions")
# ...populate `universal` with candidate cobra.Reaction objects...

with model:
    # Create a gap for demonstration
    model.remove_reactions([model.reactions.PGI])

    # Find the minimal set of reactions from `universal` to restore feasibility
    solution = gapfill(model, universal, demand_reactions=False)
    for reaction in solution[0]:
        print(reaction.id)

注： gapfillは反応リストのリストを返します（iterations > 1の場合、各反復につき1つ）。そのため、solutionを直接出力するのではなく、solution[0]を反復処理してください。

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
model = load_model("iJO1366")

# 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("iJO1366")

# 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 * solution.objective_value]

ワークフロー3：培地の最適化

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

# Load model
model = load_model("iJO1366")

# 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("iJO1366")

# 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オブジェクト

モデルでは、反応、代謝物、遺伝子にDictListオブジェクトを使用します。これはリストと辞書の両方のように動作します。

# Access by index
first_reaction = model.reactions[0]

# Access by ID
pfk = model.reactions.get_by_id("PFK")

# Query methods
atp_reactions = model.reactions.query("atp")

フラックスの制約

反応の境界値によって、実行可能なフラックス範囲が定義されます。

• 不可逆：lower_bound = 0, upper_bound > 0
• 可逆：lower_bound < 0, upper_bound > 0
• 不整合を避けるため、.bounds を使って両方の境界値を同時に設定します

遺伝子・反応ルール（GPR）

遺伝子と反応を結び付けるブール論理：

# AND logic (both required)
reaction.gene_reaction_rule = "gene1 and gene2"

# OR logic (either sufficient)
reaction.gene_reaction_rule = "gene1 or gene2"

# Complex logic
reaction.gene_reaction_rule = "(gene1 and gene2) or (gene3 and gene4)"

交換反応

代謝物のインポート／エクスポートを表す特殊な反応：

• 慣例として、EX_ を接頭辞に付けた名前を使用します
• 正のフラックス＝分泌、負のフラックス＝取り込み
• model.medium 辞書を通じて管理します

ベストプラクティス

1. コンテキストマネージャーを使用することで、一時的な変更による状態管理の問題を回避します
2. 解析の前に model.slim_optimize() を使ってモデルを検証し、実行可能性を確認します
3. 最適化後に解のステータスを確認します。optimal は正常に解が得られたことを示します
4. 熱力学的な実行可能性が重要な場合は、ループレス FVA を使用します
5. FVA で最適でない領域を探索するには、fraction_of_optimum を適切に設定します
6. 計算コストの高い操作（サンプリング、二重欠失）を並列化します
7. モデルの交換と長期保存には、SBML 形式を優先します
8. パフォーマンス向上のため、目的関数値のみが必要な場合は slim_optimize() を使用します
9. 数値的な安定性を確保するため、フラックスサンプルを検証します

トラブルシューティング

実行不可能な解：培地の制約、反応の境界値、モデルの整合性を確認してください
最適化が遅い：model.solver を通じて、別のソルバー（GLPK、CPLEX、Gurobi）を試してください
非有界の解：交換反応に適切な上限が設定されていることを確認してください
インポートエラー：正しいファイル形式と有効な SBML 識別子を使用していることを確認してください

参考資料

詳細なワークフローと API の使用パターンについては、以下を参照してください。

• references/workflows.md - ステップごとの包括的なワークフロー例
• references/api_quick_reference.md - 一般的な関数シグネチャと使用パターン

公式ドキュメント：https://cobrapy.readthedocs.io/en/latest/

AlterLab Academic Skills スイートの一部です。