cobrapy

Par mkurman · zorai

Modélisation métabolique par contraintes (COBRA). FBA, FVA, knockouts de gènes, échantillonnage de flux, modèles SBML, pour l'analyse en biologie des systèmes et en ingénierie métabolique.

npx skills add https://github.com/mkurman/zorai --skill cobrapy

COBRApy - Constraint-Based Reconstruction and Analysis

Overview

COBRApy est une bibliothèque Python pour la reconstruction et l'analyse basée sur les contraintes (COBRA) de modèles métaboliques, essentielle pour la recherche en biologie des systèmes. Travaillez avec des modèles métaboliques à l'échelle du génome, effectuez des simulations informatiques du métabolisme cellulaire, menez des analyses d'ingénierie métabolique et prédisez les comportements phénotypiques.

Core Capabilities

COBRApy fournit des outils complets organisés en plusieurs domaines clés :

1. Model Management

Chargez les modèles existants à partir de référentiels ou de fichiers :

from cobra.io import load_model

# Load bundled test models
model = load_model("textbook")  # E. coli core model
model = load_model("ecoli")     # Full E. coli model
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")

Enregistrez les modèles dans différents formats :

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. Model Structure and Components

Accédez et inspectez les composants du modèle :

# 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. Flux Balance Analysis (FBA)

Effectuez une simulation FBA standard :

# 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 parcimonieuse (minimiser le flux total) :

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

FBA géométrique (trouver la solution centrale) :

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

4. Flux Variability Analysis (FVA)

Déterminez les plages de flux pour toutes les réactions :

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. Gene and Reaction Deletion Studies

Effectuez des analyses de knockout :

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. Growth Media and Minimal Media

Gérez le milieu de croissance :

# 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. Flux Sampling

Échantillonnez l'espace de flux réalisable :

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. Production Envelopes

Calculez les plans de phase phénotypique :

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. Gapfilling

Ajoutez des réactions pour rendre les modèles réalisables :

from cobra.flux_analysis import gapfill

# Prepare universal model with candidate reactions
universal = load_model("universal")

# 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. Model Building

Construisez des modèles à partir de zéro :

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"

Common Workflows

Workflow 1: Load Model and Predict Growth

from cobra.io import load_model

# Load model
model = load_model("ecoli")

# 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])

Workflow 2: Gene Knockout Screen

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

# Load model
model = load_model("ecoli")

# 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]

Workflow 3: Media Optimization

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

# Load model
model = load_model("ecoli")

# 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)

Workflow 4: Flux Uncertainty Analysis

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

# Load model
model = load_model("ecoli")

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

Workflow 5: Context Manager for Temporary Changes

Utilisez des context managers pour effectuer des modifications temporaires :

# 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}")

Key Concepts

DictList Objects

Les modèles utilisent des objets DictList pour les réactions, métabolites et gènes - se comportant à la fois comme des listes et des dictionnaires :

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

Flux Constraints

Les bornes de réaction définissent les plages de flux réalisables :

  • Irréversible : lower_bound = 0, upper_bound > 0
  • Réversible : lower_bound < 0, upper_bound > 0
  • Définissez les deux bornes simultanément avec .bounds pour éviter les incohérences

Gene-Reaction Rules (GPR)

Logique booléenne reliant les gènes aux réactions :

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

Exchange Reactions

Réactions spéciales représentant l'import/export de métabolites :

  • Nommées avec le préfixe EX_ par convention
  • Flux positif = sécrétion, flux négatif = absorption
  • Gérées via le dictionnaire model.medium

Best Practices

  1. Utilisez les context managers pour les modifications temporaires afin d'éviter les problèmes de gestion d'état
  2. Validez les modèles avant l'analyse en utilisant model.slim_optimize() pour assurer la réalisabilité
  3. Vérifiez le statut de la solution après l'optimisation - optimal indique une résolution réussie
  4. Utilisez FVA loopless quand la réalisabilité thermodynamique importe
  5. Définissez fraction_of_optimum correctement dans FVA pour explorer l'espace sous-optimal
  6. Parallélisez les opérations coûteuses en calcul (échantillonnage, doubles deletions)
  7. Préférez le format SBML pour l'échange de modèles et le stockage à long terme
  8. Utilisez slim_optimize() quand seule la valeur objective est nécessaire pour la performance
  9. Validez les échantillons de flux pour assurer la stabilité numérique

Troubleshooting

Solutions infaisables : Vérifiez les contraintes du milieu, les bornes de réaction et la cohérence du modèle Optimisation lente : Essayez différents solveurs (GLPK, CPLEX, Gurobi) via model.solver Solutions non bornées : Vérifiez que les réactions d'échange ont des bornes supérieures appropriées Erreurs d'importation : Assurez-vous du format de fichier correct et des identifiants SBML valides

References

Pour des workflows détaillés et des patterns API, consultez :

  • references/workflows.md - Exemples de workflow complets et détaillés
  • references/api_quick_reference.md - Signatures de fonction courantes et patterns

Documentation officielle : https://cobrapy.readthedocs.io/en/latest/

Skills similaires