cobrapy-metabolic-modeling

COBRApy — 制約ベースの代謝モデリング

概要

COBRApy は、ゲノムスケールの代謝モデルを対象とする制約ベースの再構築および解析（COBRA）用の Python パッケージです。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. モデルの入出力

同梱モデルの読み込みと、標準形式の読み書きを行います。

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 object	フラックスエンベロープとともに炭素収率を計算します
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: COnstraints-Based Reconstruction and Analysis for Python」 — BMC Systems Biology