executing-plans
executing-plans
Use when executing an implementation plan in the current session as the implementer yourself — your human partner chose inline execution, or no subagent tool is available
bio-systems-biology-metabolic-reconstructionBuilds draft genome-scale metabolic models from an annotated genome using CarveMe (top-down carving of a BiGG universal model) or gapseq (bottom-up pathway-evidence reconstruction), then loads and sanity-checks the draft in COBRApy. Use when creating a model for an organism without one, choosing between CarveMe and gapseq, gap-filling to a target medium, un…
Independent GitHub snapshot · Automated classification · Original documentation governs
Reference examples tested with: CarveMe 1.6+, gapseq 1.2+, COBRApy 0.29+, DIAMOND 2.1+, Python 3.10+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signatures<tool> --version then <tool> --help to confirm flagsIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Note: CarveMe needs an LP solver (academic CPLEX/Gurobi; SCIP is a slow open-source fallback) and a DIAMOND install, and its universal model is BiGG-derived so the BiGG-model release matters. gapseq is cloned from GitHub (not pip-installable), emits ModelSEED-namespace models, and its reference DB version matters. Model predictions are only comparable within the same tool, DB, and namespace.
"Build a metabolic model from my organism's genome" -> Map the annotated proteome/genome to reactions in a reference database, assemble a draft network with a biomass reaction, and gap-fill so it can grow on a chosen medium.
carve genome.faa -o model.xml (CarveMe, top-down); gapseq doall genome.fna (gapseq, bottom-up)Automated reconstruction produces a DRAFT, not a finished model. The single most misleading signal is growth: CarveMe and gapseq GAP-FILL specifically to force biomass production on a chosen medium, so a draft that grows proves nothing biological - it was made to grow. Consequences:
| Goal | Tool | Why / trade-off |
|---|---|---|
| Fast draft(s) for well-studied bacteria; batch/community | CarveMe (carve) | top-down MILP carving of a curated BiGG universe; minutes; universe is simulation-ready but BiGG-centric; universal biomass; weak transporters |
| Non-model/environmental clade; carbon-source & fermentation phenotypes | gapseq | bottom-up homology + pathway-completeness; slower, more transparent; better SCFA/carbon-use recovery; ModelSEED namespace complicates merging |
| Fully-automated web pipeline (RAST annotation) | ModelSEED/KBase | template-based; convenient; template biomass and aggressive gap-fill can force implausible reactions |
| Eukaryotes / fungi / actinomycetes | RAVEN (MATLAB) | KEGG/MetaCyc-based, template or de novo; MATLAB license; the eukaryote-capable option |
Do NOT treat "CarveMe and gapseq do the same thing, pick the faster one" as true: different philosophies, namespaces (BiGG vs ModelSEED), and failure modes. The choice is scientific. No single tool dominates - which is why consensus/ensemble reconstruction exists.
pip install carveme # also needs DIAMOND and an LP solver (CPLEX/Gurobi; SCIP fallback)
# Draft from a PROTEIN FASTA (default input). Raw/GenBank genomes are NOT accepted.
carve genome.faa -o model.xml
# Gram type and universe are VALUES of -u/--universe, NOT --grampos/--gramneg flags.
carve genome.faa -o model.xml -u grampos # {bacteria (default), grampos, gramneg, archaea, cyanobacteria}
# Gap-fill to force growth on a medium (opt-in; records what was added for that medium).
carve genome.faa -o model.xml --gapfill M9
carve genome.faa -o model.xml -u gramneg --gapfill M9,LB # multiple media
# Nucleotide input instead of protein, or download by accession:
carve genome.fna --dna -o model.xml
Community reconstruction uses a SEPARATE merge_community command (not carve); see systems-biology/community-metabolic-modeling.
git clone https://github.com/jotech/gapseq && cd gapseq && ./gapseq test # cloned, not pip; check deps
# One-shot: find + find-transport + draft + fill
./gapseq doall genome.fna
# Or the explicit steps (note find-transport is its OWN subcommand, not `find -t`):
./gapseq find -p all genome.fna # -> genome-all-Reactions.tbl, genome-all-Pathways.tbl
./gapseq find-transport genome.fna # -> genome-Transporter.tbl (singular)
./gapseq draft -r genome-all-Reactions.tbl -t genome-Transporter.tbl \
-p genome-all-Pathways.tbl -c genome.fna # -> genome-draft.RDS, genome-rxnWeights.RDS
./gapseq fill -m genome-draft.RDS -n dat/media/M9.csv \
-c genome-rxnWeights.RDS -g genome-rxnXgenes.RDS # -> genome.xml / genome.RDS
Goal: Read the draft, confirm it grows on the gap-fill medium, and inventory the parts most likely to be wrong.
Approach: Load the SBML into COBRApy, report network size and gene coverage, test growth, and count orphan (gene-less) reactions and exchanges - the draft's soft spots before curation.
import cobra
model = cobra.io.read_sbml_model('model.xml')
print(f'reactions={len(model.reactions)} metabolites={len(model.metabolites)} genes={len(model.genes)}')
print(f'grows on gap-fill medium: {model.slim_optimize() > 1e-3}') # true by construction if gap-filled
orphans = [r for r in model.reactions if not r.genes] # no GPR: gap-filled, spontaneous, or transport
print(f'orphan (gene-less) reactions: {len(orphans)} exchanges: {len(model.exchanges)}')
# Typical bacterial draft: ~1000-2500 reactions. Far outside that range flags an annotation problem.
# Reaction/metabolite IDs come from the tool's reference DB: CarveMe = BiGG, gapseq/ModelSEED =
# ModelSEED (seed.*), RAVEN = KEGG/MetaCyc. Two models in different namespaces cannot be merged or
# compared directly. Reconcile through MetaNetX/MNXref (MNXM* metabolites, MNXR* reactions) BEFORE
# any cross-tool merge or community build. This BiGG-vs-ModelSEED split is exactly why community
# modeling of CarveMe + gapseq outputs breaks without reconciliation.
| Symptom | Cause | Fix |
|---|---|---|
carve errors on a genome file | GenBank/nucleotide passed where protein FASTA expected | supply a protein FASTA, or add --dna for nucleotide |
--grampos/--gramneg not recognized | those are -u/--universe VALUES, not flags | carve ... -u grampos |
| Draft cannot grow at all | no gap-filling requested, or wrong medium | add --gapfill <medium>; confirm the medium supplies biomass precursors |
| Draft grows on everything / implausibly | gap-fill forced reactions for the chosen medium | flag gap-filled reactions low-confidence; re-gap-fill on the correct medium; curate |
| Two models will not merge / IDs mismatch | different namespaces (BiGG vs ModelSEED) | reconcile via MetaNetX/MNXref before merging |
gapseq find -t fails | transport is the find-transport subcommand | use ./gapseq find-transport genome.fna |
| Very few genes / tiny network | poor annotation or wrong input file | check the proteome/annotation; verify gene IDs |