Orchestrates genome-scale metabolic modeling from a protein FASTA to flux predictions, chaining CarveMe/gapseq reconstruction, memote QC, gap-filling, media-constrained FBA/FVA, gene essentiality, and context-specific models. Use when committing the reconstruction tool (which locks the identifier NAMESPACE forever - BiGG vs ModelSEED vs KEGG, no automatic translation), setting the medium BEFORE FBA (the exchange bounds ARE the medium; essentiality and gap-fill are computed relative to it), curating iteratively (stoichiometric-consistency first, then mass/charge, then directionality, then GPR) with energy-generating-cycle removal, and reading a MEMOTE score as well-formedness NOT correctness. Hands mechanism to the systems-biology component skills; not a re-teach of any single step.
Orchestrates genome-scale metabolic modeling from a protein FASTA to flux predictions, chaining CarveMe/gapseq reconstruction, memote QC, gap-filling, media-constrained FBA/FVA, gene essentiality, and context-specific models. Use when committing the reconstruction tool (which locks the identifier NAMESPACE forever - BiGG vs ModelSEED vs KEGG, no automatic translation), setting the medium BEFORE FBA (the exchange bounds ARE the medium; essentiality and gap-fill are computed relative to it), curating iteratively (stoichiometric-consistency first, then mass/charge, then directionality, then GPR) with energy-generating-cycle removal, and reading a MEMOTE score as well-formedness NOT correctness. Hands mechanism to the systems-biology component skills; not a re-teach of any single step.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Metabolic Modeling Pipeline
"Build and analyze a metabolic model for my organism" -> Orchestrate CarveMe reconstruction, memote quality scoring, gap-filling, FBA/FVA flux analysis, gene essentiality prediction, and context-specific model building from expression data.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
The governing principle
An automated reconstruction is a HYPOTHESIS about metabolism, not a finished model — the tool does only stage 1 of a four-stage process (Thiele & Palsson) and makes a 30-second draft look finished. Three commitments are made before flux ever means anything, and a high MEMOTE score does NOT certify any of them.
The reconstruction tool locks the identifier NAMESPACE forever. D-glucose is glc__D_c (BiGG, CarveMe) = cpd00027_c0 (ModelSEED, gapseq) = C00031 (KEGG). Mixing namespaces (merging a BiGG CarveMe model with a ModelSEED gapseq model) is a silent, error-free failure: metabolites fail to string-match, producing duplicated metabolites, disconnected reactions, and wrong growth. MetaNetX/MNXref is the required, non-automatic reconciliation layer before ANY cross-tool merge.
The medium IS the model, and it is set before FBA means anything. The exchange-reaction bounds (EX_*_e lower bounds = uptake) define what the cell can do; FBA on a wrong/default medium invalidates everything, and gap-fill + essentiality are computed RELATIVE to the medium (a biosynthetic gene is essential in minimal, non-essential in rich). The in-silico medium must match the wet-lab condition being validated against.
The biomass objective function drives every flux and every knockout call. The BOF is a gene-less pseudo-reaction whose flux = growth rate; copying E. coli biomass into a non-model organism carries the wrong physiology. Two ATP-maintenance terms both matter (GAM inside biomass; NGAM as a separate ATPM with a fixed floor — ATPM lower bound = 0 is a red flag).
MEMOTE 95% is well-formedness, not correctness. It scores stoichiometric consistency + annotation coverage + SBO, and is partly gameable; the predictive tests (energy-generating-cycle sweep, essentiality, Biolog) are NOT in the scored total. A 95% model can still make ATP from nothing.
Every metabolite ID; a cross-namespace merge silently duplicates metabolites and breaks growth
Input = protein FASTA with genes split (CarveMe) vs genome FASTA (gapseq)
Whether the reconstruction runs at all; feeding a raw assembly to CarveMe is the classic trap
Medium (exchange bounds)
What the cell can do; gap-fill and essentiality are computed relative to it
Biomass objective function + GAM/NGAM
Every flux and knockout call; a foreign or zero-ATPM BOF invalidates essentiality
Workflow Overview
Protein FASTA (genome annotation)
|
v
[1. Reconstruction] --> CarveMe / gapseq / ModelSEED
|
v
[2. Model Validation] --> memote QC (snapshot report)
|
v
[3. Model Curation] --> gap-filling, mass/charge balance
|
| <---- Iterative refinement loop
v
[4. FBA Analysis] --> Growth prediction, flux distribution
|
+-----------------------+
| |
v v
[5a. Gene Essentiality] [5b. Context-Specific]
Single/double KO Tissue-specific models
| |
v v
Essential Gene List Condition-Specific Fluxes
CarveMe's MILP carving needs a commercial solver (CPLEX or Gurobi, both free for academics) via reframed (CarveMe 1.5+ depends on reframed, the refactored framed) — the #1 install pain; the open SCIP fallback is far slower. gapseq (GLPK default) is the solver-free reconstruction alternative.
Required data:
Protein FASTA file from genome annotation
BiGG universal model (downloaded by CarveMe)
Primary Path: Bacterial Model from Genome
Step 1: Automated Reconstruction with CarveMe
# Basic reconstruction from a PROTEIN FASTA (CarveMe rejects raw/GenBank genomes)
carve genome.faa -o model_draft.xml
# Gram type / universe is a VALUE of -u/--universe, NOT a --gram-neg flag
carve genome.faa -o model_draft.xml -u gramneg
# Gap-fill for a specific medium (opt-in; the medium determines what gets added)
carve genome.faa -o model_draft.xml -u gramneg --gapfill M9
# Run the memote test suite (results stored as JSON when --filename is given)
memote run --filename model_result.json.gz model_draft.xml
# Generate the human-readable HTML snapshot report
memote report snapshot --filename model_report.html model_draft.xml
# The memote SCORE measures consistency and annotation (well-formedness), NOT biological# correctness -- a model can score high and mispredict every knockout. Read WHICH tests fail# (stoichiometric consistency, mass/charge balance, energy-generating cycles) in the HTML report,# and validate predictions separately (see systems-biology/model-curation). Programmatic access:from memote.suite.api import test_model
code, result = test_model(model, results=True) # result is a MemoteResult of the raw outcomes
Step 3: Model Curation (Iterative)
Curation order matters and the steps interact: fix stoichiometric consistency FIRST (unconserved metabolites poison everything), then per-reaction mass/charge balance, then directionality, then GPR. Gap-filling can break a balance and an energy-generating-cycle fix frequently unmasks a SECOND nested EGC, so re-run memote + the EGC sweep + a growth check after every change. Mechanism lives in systems-biology/model-curation.
import cobra
from cobra.flux_analysis import gapfill
model = cobra.io.read_sbml_model('model_draft.xml')
# Check for common issuesdefdiagnose_model(model):
issues = []
# Dead-end metabolites (produced but not consumed, or vice versa)for met in model.metabolites:
producing = [r for r in met.reactions if met in r.products]
consuming = [r for r in met.reactions if met in r.reactants]
iflen(producing) > 0andlen(consuming) == 0:
issues.append(f'Dead-end (not consumed): {met.id}')
eliflen(producing) == 0andlen(consuming) > 0:
issues.append(f'Dead-end (not produced): {met.id}')
# Blocked reactions: reactions that CANNOT carry flux under any feasible state. fraction_of_optimum# defaults to 1.0, which instead pins growth at the optimum and reports reactions unused by that# particular optimal solution -- a different, much larger set. Pass 0 (or use find_blocked_reactions).
fva = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0)
blocked = fva[(fva['minimum'] == 0) & (fva['maximum'] == 0)]
iflen(blocked) > 0:
issues.append(f'Blocked reactions: {len(blocked)}')
return issues
issues = diagnose_model(model)
print(f'Found {len(issues)} issues')
for issue in issues[:10]:
print(f' {issue}')
# Gap-filling for growth on specific mediafrom cobra.flux_analysis import gapfill
# Load universal reaction database for gap-filling
universal = cobra.io.read_sbml_model('universal_model.xml')
# Define target medium (e.g., glucose minimal)
target_medium = {
'EX_glc__D_e': 10, # Glucose uptake'EX_o2_e': 20, # Oxygen'EX_nh4_e': 100, # Ammonium'EX_pi_e': 100, # Phosphate'EX_so4_e': 100, # Sulfate
}
# Apply medium (model.exchanges yields Reaction objects, not id strings). This is a FULLY DEFINED# medium: every exchange absent from target_medium is closed, so target_medium must also list the# trace metals and cofactors biomass requires (Fe, K, Mg, Ca, Zn, ...) or growth is zero. To vary only# the carbon source instead, layer overrides onto `model.medium` rather than replacing it.for rxn in model.exchanges:
rxn.lower_bound = -target_medium[rxn.id] if rxn.idin target_medium else0# block other uptakes# Gap-fill to enable growth# Gap-filling adds minimal reactions from universal model to enable growth
gapfill_solution = gapfill(model, universal, demand_reactions=False)
print(f'Gap-fill added {len(gapfill_solution[0])} reactions')
# Gap-filled reactions are the LEAST-evidenced part of the model (added to force growth on this# medium, not because homology supports them) -- flag them low-confidence, do not treat as validated.for rxn in gapfill_solution[0]:
model.add_reactions([rxn])
print(f' Added (low-confidence): {rxn.id} - {rxn.name}')
# Verify growth
solution = model.optimize()
print(f'Growth after gap-fill: {solution.objective_value:.4f} h^-1')
# Persist the curated model -- downstream steps read model_curated.xml, not the draft
cobra.io.write_sbml_model(model, 'model_curated.xml')
Step 4: Flux Balance Analysis
import cobra
import pandas as pd
import matplotlib.pyplot as plt
model = cobra.io.read_sbml_model('model_curated.xml')
# Basic FBA
solution = model.optimize()
print(f'Objective (growth): {solution.objective_value:.4f} h^-1')
print(f'Status: {solution.status}')
# Get active fluxes
fluxes = solution.fluxes
active_fluxes = fluxes[abs(fluxes) > 1e-6]
print(f'Active reactions: {len(active_fluxes)} / {len(model.reactions)}')
# Key exchange fluxes (uptake/secretion)
exchange_fluxes = fluxes[[r.idfor r in model.exchanges]]
significant_exchanges = exchange_fluxes[abs(exchange_fluxes) > 0.1]
print('\nSignificant exchanges:')
print(significant_exchanges.sort_values())
from cobra.flux_analysis import single_gene_deletion, double_gene_deletion
# Single gene knockouts. Result columns: ids (a SET of gene-id strings), growth, status.
single_ko = single_gene_deletion(model)
single_ko['growth_ratio'] = single_ko['growth'] / solution.objective_value
# Essential genes: knockout drops growth below the cutoff (a policy, not a library default).# Report and sweep the cutoff; match the medium to any experiment being compared. Essentiality# is model- and medium-relative (see systems-biology/gene-essentiality).
essential = single_ko[single_ko['growth_ratio'] < 0.1]
print(f'Essential genes: {len(essential)} / {len(model.genes)} on this medium')
# ids elements are gene-id STRINGS (a set), so list(s)[0] gives the id -- there is no .id attribute.
essential_list = [list(s)[0] for s in essential['ids']]
withopen('essential_genes.txt', 'w') as f:
f.write('\n'.join(essential_list))
# Double gene knockouts (synthetic lethality)# WARNING: Computationally intensive for large models# Focus on non-essential genes only (a synthetic lethal needs both singles viable)
non_essential = [g.idfor g in model.genes if g.idnotin essential_list]
# Run pairwise deletions (positional gene_list1/gene_list2; cap the O(n^2) sweep)
double_ko = double_gene_deletion(model, non_essential[:100], non_essential[:100])
# Synthetic lethality: neither single KO is lethal, but the double KO is
synthetic_lethal = double_ko[double_ko['growth'] < 0.01]
print(f'Synthetic lethal pairs: {len(synthetic_lethal)}')
Step 5b: Context-Specific Models
Use a validated extraction method rather than ad-hoc pruning. COBRApy has no native GIMME/iMAT/INIT; the real Python options are troppo and corda, and the threshold/method choice dominates the result more than the data does. See systems-biology/context-specific-models for the method decision table and the threshold-sensitivity discipline.
# corda is the most turnkey native-Python extraction method.from corda import CORDA, reaction_confidence
# Translate expression into CORDA confidence classes (-1 absent, 0 unknown, 1 low, 2 med, 3 high)# through the GPR, then build the context model.
gene_conf = {g.id: 2for g in model.genes} # derive from expression quantiles
rxn_conf = {r.id: reaction_confidence(r, gene_conf) for r in model.reactions} # corda 0.5+ takes the Reaction object (reads r.gpr); the old GPR-string form was removed
opt = CORDA(model, rxn_conf)
opt.build()
context_model = opt.cobra_model('tissue')
print(f'Context model: {len(context_model.reactions)} reactions')
# Rebuild at 2-3 thresholds and report which reactions are threshold-dependent (hypotheses).
Visualization with Escher
import escher
# Load model and solution
model = cobra.io.read_sbml_model('model_curated.xml')
solution = model.optimize()
# Create Escher map
builder = escher.Builder(
map_name='e_coli_core.Core metabolism',
model=model,
reaction_data=solution.fluxes.to_dict()
)
builder.save_html('flux_map.html')
Reconcile via MetaNetX/MNXref before any merge; never string-join across namespaces
Every yield and essentiality inflated
Energy-generating cycle (free ATP from over-permissive reversibility + gap-fill)
Close exchanges + add ATP demand, confirm max=0; constrain directionality (eQuilibrator dG); EGCs nest, re-test
Unrealistic growth rate
Unbounded/wrong uptake
Audit the full boundary set; set model.medium (POSITIVE magnitudes) to mirror the assay
MEMOTE 95% but predictions wrong
MEMOTE scores well-formedness, not correctness
Treat MEMOTE as a hygiene floor; add EGC + essentiality/Biolog validation (not in the scored total)
Extra "essential" genes vs literature
Essentiality reported without stating the medium
In-silico medium = wet-lab medium; report and SWEEP the cutoff (1/2/5/10% WT growth)
Gap-filled reactions trusted as real
Gap-fill adds reactions to force growth, not from evidence
Flag low-confidence, keep distinguishable, prioritize for experimental verification
No growth
Missing reactions on THIS medium
Gap-fill FOR the specified medium; do not gap-fill to one medium then predict on another
References
Machado D, Andrejev S, Tramontano M, Patil KR (2018) Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Research 46:7542-7553. DOI 10.1093/nar/gky537. (CarveMe.)
Zimmermann J, Kaleta C, Waschina S (2021) gapseq: informed prediction of bacterial metabolic pathways and reconstruction of accurate metabolic models. Genome Biology 22:81. DOI 10.1186/s13059-021-02295-1.
Lieven C, Beber ME, Olivier BG, et al (2020) MEMOTE for standardized genome-scale metabolic model testing. Nature Biotechnology 38:272-276. DOI 10.1038/s41587-020-0446-y.
Thiele I, Palsson BO (2010) A protocol for generating a high-quality genome-scale metabolic reconstruction. Nature Protocols 5:93-121. DOI 10.1038/nprot.2009.203. (four-stage reconstruction.)
Fritzemeier CJ, Hartleb D, Szappanos B, Papp B, Lercher MJ (2017) Erroneous energy-generating cycles in published genome scale metabolic networks: identification and removal. PLoS Computational Biology 13:e1005494. DOI 10.1371/journal.pcbi.1005494.
Output Files
File
Description
model_draft.xml
Initial reconstruction (SBML)
model_curated.xml
Gap-filled and validated model
model_report.html
Memote QC report
essential_genes.txt
Predicted essential genes
gene_essentiality.tsv
Full single-gene-deletion table (growth ratio per gene)
fba_fluxes.tsv
Optimal flux distribution
fva_results.tsv
Flux variability ranges
model_analysis_summary.pdf / .png
Summary figure (exchanges, FVA, essentiality)
flux_map.html
Escher visualization
Extensions
Beyond the core genome-to-flux path, the model feeds two further analyses: build a multi-species community from several reconstructions (systems-biology/community-metabolic-modeling), or design growth-coupled knockouts to overproduce a target chemical (systems-biology/strain-design).
