End-to-end genome-scale metabolic modeling from genome sequence to flux predictions. Covers automated reconstruction with CarveMe, model validation with memote, FBA/FVA analysis, and gene essentiality prediction. Use when building metabolic models or predicting metabolic phenotypes from genomic data.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
End-to-end genome-scale metabolic modeling from genome sequence to flux predictions. Covers automated reconstruction with CarveMe, model validation with memote, FBA/FVA analysis, and gene essentiality prediction. Use when building metabolic models or predicting metabolic phenotypes from genomic data.
[{"after_reconstruction":"Reactions 1000-2500, growth >0.01 on target media"},{"after_curation":"Memote score >50%, <5% orphan reactions"},{"after_fba":"Realistic growth rate, major pathways active"},{"after_essentiality":"Core essential genes match literature >70%"}]
Metabolic Modeling Pipeline
Complete workflow for genome-scale metabolic modeling: from protein sequences to flux predictions and phenotype analysis.
Workflow Overview
Protein FASTA (genome annotation)
|
v
[1. Reconstruction] --> CarveMe / gapseq / ModelSEED
|
v
[2. Model Curation] --> memote QC, gap-filling
|
| <---- Iterative refinement loop
v
[3. FBA Analysis] --> Growth prediction, flux distribution
|
+-----------------------+
| |
v v
[4a. Gene Essentiality] [4b. Context-Specific]
Single/double KO Tissue-specific models
| |
v v
Essential Gene List Condition-Specific Fluxes
# Basic reconstruction from protein sequences
carve genome.faa -o model_draft.xml
# With gram type specification (improves biomass composition)
carve genome.faa -o model_draft.xml --gram-neg
# Gap-fill for specific media
carve genome.faa -o model_draft.xml --gram-neg --gapfill M9
# Run memote QC
memote run --filename model_draft_report.html model_draft.xml
# Generate snapshot for comparison
memote report snapshot --filename model_snapshot.json model_draft.xml
import json
withopen('model_snapshot.json') as f:
report = json.load(f)
# Key metrics to check
metrics = report['tests']['basic']
print('=== Model QC Metrics ===')
print(f"Reactions with genes: {metrics.get('test_gene_reaction_rule_presence', {}).get('metric', 'N/A')}")
print(f"Metabolites in reactions: {metrics.get('test_metabolites_not_produced', {}).get('metric', 'N/A')}")
print(f"Stoichiometry balance: {metrics.get('test_stoichiometric_consistency', {}).get('result', 'N/A')}")
# Target: memote score >50% for usable model# Score >70% indicates well-curated model
total_score = report.get('score', {}).get('total_score', 0)
print(f'Total memote score: {total_score:.1%}')
Step 3: Model Curation (Iterative)
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
fva = cobra.flux_analysis.flux_variability_analysis(model)
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.medium import minimal_medium
# 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 mediumfor rxn_id in model.exchanges:
rxn = model.reactions.get_by_id(rxn_id)
if rxn_id in target_medium:
rxn.lower_bound = -target_medium[rxn_id]
else:
rxn.lower_bound = 0# 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')
for rxn in gapfill_solution[0]:
model.add_reactions([rxn])
print(f' Added: {rxn.id} - {rxn.name}')
# Verify growth
solution = model.optimize()
print(f'Growth after gap-fill: {solution.objective_value:.4f} h^-1')
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
single_ko = single_gene_deletion(model)
single_ko['growth_ratio'] = single_ko['growth'] / solution.objective_value
# Essential genes: knockout abolishes growth (<10% of WT)# Threshold 0.1 is standard for essentiality
essential = single_ko[single_ko['growth_ratio'] < 0.1]
print(f'Essential genes: {len(essential)} / {len(model.genes)}')
# Export essential genes
essential_list = [list(g)[0].idfor g in essential.index]
withopen('essential_genes.txt', 'w') as f:
f.write('\n'.join(essential_list))
# Compare to experimental data (if available)# Typically expect >70% overlap with experimental essentiality screens
# Double gene knockouts (synthetic lethality)# WARNING: Computationally intensive for large models# Focus on non-essential genes only
non_essential = [g.idfor g in model.genes if g.idnotin essential_list]
# Run pairwise deletions
double_ko = double_gene_deletion(model, gene_list=non_essential[:100]) # Subset for speed# Find synthetic lethal pairs# Synthetic lethality: neither single KO is lethal, but 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
# Build tissue-specific model using expression datadefbuild_context_model(model, expression_data, threshold_percentile=25):
'''Build context-specific model by removing lowly-expressed reactions.
threshold_percentile: reactions below this expression percentile are candidates for removal
'''import numpy as np
context_model = model.copy()
threshold = np.percentile(list(expression_data.values()), threshold_percentile)
reactions_to_remove = []
for rxn in context_model.reactions:
if rxn.gene_reaction_rule:
genes_in_rxn = [g.idfor g in rxn.genes]
# Get expression for genes in reaction
expr_values = [expression_data.get(g, 0) for g in genes_in_rxn]
# If all genes below threshold, candidate for removalifall(e < threshold for e in expr_values):
# Check if removal breaks growthwith context_model:
rxn.knock_out()
sol = context_model.optimize()
if sol.objective_value > 0.01:
reactions_to_remove.append(rxn.id)
for rxn_id in reactions_to_remove:
context_model.reactions.get_by_id(rxn_id).remove_from_model()
return context_model
# Example: liver-specific model
liver_expression = {'gene1': 100, 'gene2': 5, 'gene3': 50} # TPM values
liver_model = build_context_model(model, liver_expression)
print(f'Liver model: {len(liver_model.reactions)} reactions')
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')
