GPTomics/bio-systems-biology-metabolic-reconstruction

Version Compatibility

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:

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

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.

Metabolic Reconstruction

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

• CLI: carve genome.faa -o model.xml (CarveMe, top-down); gapseq doall genome.fna (gapseq, bottom-up)

The governing principle: a draft is a hypothesis, and "it grows" is guaranteed by construction

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:

• Gap-filled reactions are the least-evidenced part of the model (added to close a hole, not because homology supports them), and the gap-fill medium determines what gets added. Gap-filling on the wrong medium bakes in the wrong reactions. Flag gap-filled reactions as low-confidence and record the medium.
• Draft quality is bounded by the annotation and the reference database. CarveMe can only ever include reactions in the BiGG universe (biased toward well-studied organisms); gapseq's homology thresholds and pathway logic set its floor. Peripheral/novel metabolism is systematically underrepresented.
• The draft is the START of curation, not the end. Different tools produce markedly different models from the same genome, and successive models of the same organism (iJR904 -> iJO1366 -> iML1515) give different predictions. Reconstruction feeds model-curation, never bypasses it.

Decision: CarveMe vs gapseq vs ModelSEED

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

CarveMe (top-down)

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.

gapseq (bottom-up, pathway-evidence)

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

Load and Sanity-Check the Draft

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.

Namespaces (the silent killer of model comparison)

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

Common Errors

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

Related Skills

• systems-biology/model-curation - Curate, gap-fill deliberately, and validate the draft (the required next step)
• systems-biology/flux-balance-analysis - Predict growth/flux once the model is trustworthy
• systems-biology/community-metabolic-modeling - Combine reconstructions into a community model
• genome-annotation/prokaryotic-annotation - Produce the annotated protein FASTA CarveMe/gapseq consume
• database-access/ncbi-datasets-cli - Fetch genome/proteome inputs

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 Res 46(15):7542-7553. (CarveMe)
• Zimmermann J, Kaleta C, Waschina S. 2021. gapseq: informed prediction of bacterial metabolic pathways and reconstruction of accurate metabolic models. Genome Biol 22(1):81.
• Henry CS, DeJongh M, Best AA, et al. 2010. High-throughput generation, optimization and analysis of genome-scale metabolic models. Nat Biotechnol 28(9):977-982. (ModelSEED)
• Wang H, Marcisauskas S, Sanchez BJ, et al. 2018. RAVEN 2.0: a versatile toolbox for metabolic network reconstruction. PLoS Comput Biol 14(10):e1006541.
• Thiele I, Palsson BO. 2010. A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat Protoc 5(1):93-121.
• Mendoza SN, Olivier BG, Molenaar D, Teusink B. 2019. A systematic assessment of current genome-scale metabolic reconstruction tools. Genome Biol 20(1):158.
• Moretti S, Tran VDT, Mehl F, et al. 2021. MetaNetX/MNXref: unified namespace for metabolites and biochemical reactions. Nucleic Acids Res 49(D1):D570-D574.
• Feist AM, Palsson BO. 2010. The biomass objective function. Curr Opin Microbiol 13(3):344-349.
• Monk JM, Lloyd CJ, Brunk E, et al. 2017. iML1515, a knowledgebase that computes Escherichia coli traits. Nat Biotechnol 35(10):904-908.