Skip to content

Match MATLAB kcat selection; full PhylDist, ecFVA and light vs full in full_ecModel tutorial - #54

Merged
edkerk merged 6 commits into
developfrom
fix/fuzzy-origin-rank-and-tutorial
Sep 11, 2026
Merged

edkerk merged 6 commits into
developfrom
fix/fuzzy-origin-rank-and-tutorial

Conversation

@edkerk

@edkerk edkerk commented Sep 11, 2026 •

Copy link
Copy Markdown
Collaborator

Running the full_ecModel tutorial of 4.0.0b1 side by side with the MATLAB protocol (GECKO develop4 + RAVEN 3.0.0b1) showed that the Python ecModel diverged from MATLAB already at stage 2: maximum growth before tuning was 0.0013 /hour against 0.116 in MATLAB. This PR fixes the causes, all in kcat gathering, and brings the tutorial up to date. After it, the tutorial tracks MATLAB step by step.

Library fixes

  • fuzzy_kcat_matching: selection among the EC numbers of one reaction. Within one EC the levels are tried in the order 1, 2, 3, 5, 4, 6, as in MATLAB. Across the EC numbers of a reaction, MATLAB compares matches in that same search order, but geckopy compared the output origin numbers, so an any-organism any-substrate kcat (origin 4) beat a correct-organism specific activity (origin 5). New SEARCH_RANK, derived from SEARCH_LEVELS, drives both. Example: r_0005 (EC 2.4.1.34;3.6.5.2) got 0.0005 /s from a human GTPase value for the second EC; MATLAB takes 13.3 /s from the S. cerevisiae specific activity for the first. That one kcat held growth at 0.0013 /hour. This changes kcats of multi-EC reactions, hence a CHANGELOG entry.
  • load_conventional_gem read YAML with cobra.io.load_yaml_model, which drops RAVEN reaction fields. All 2411 eccodes entries of yeast-GEM were lost, so fill_eccodes_from_gem found none. It now uses raven-toolbox's read_yaml_model, which also reads plain cobra YAML.
  • New geckopy.databases.download_phyl_dist fetches RAVEN's keggPhylDist.mat (pinned to the 3.0.0b1 tag; the file is identical on develop3), validating it before it replaces the destination.

Tutorial (tutorials/full_ecModel)

  • data/PhylDist.mat was a two-organism stub, so the closest-organism levels picked arbitrary rows. For example, r_0103_REV got 5e6 /s from Cupriavidus necator where MATLAB takes 3 /s from Aspergillus fumigatus. The stub is removed and the full file is downloaded on first run (git-ignored).
  • STEP 16-17 passes kegg_db, as MATLAB's getECfromDatabase falls back to KEGG. This took EC-string mismatches with MATLAB from 256 to 4.
  • STEP 26 uses merge_kcats instead of the deprecated merge_dlkcat_and_fuzzy_kcats.
  • STEP 43-44 replaces the stale "ABC-SMC is not ported" note with a pointer to evotune (not run, as in MATLAB).
  • STEP 73-75 ecFVA now runs (about 2 min with Gurobi, skipped on GLPK). The commented-out version computed the "bare GEM" column with ec_fva(ec_model_bare, …), duplicating the ecModel column.
  • STEP 76-77 light vs full is ported (code/plot_light_vs_full.py, from plotlightVSfull.m).
  • model_adapter.toml: the max_growth_weight comment said MATLAB weights the flux term, but both weight the max-growth term identically. The same fix is in EvotuneParams.max_growth_weight; the prior_penalty_weight comment is also fixed, since GECKO 4 has the term too.
  • Docs: the three links to a "MATLAB ↔ Python translation guide" on gecko-docs returned 404 (no such page). They now point to the GECKO vs. geckopy page.

Result

MATLAB (develop4 + fixes in SysBioChalmers/GECKO) 4.0.0b1 this PR
STEP 33 max growth 0.116 0.0013 0.124 ¹
sensitivity tuning 7 kcats → 0.427 9 kcats → 0.419 7 kcats → 0.428
STEP 59-63 growth 0.0129 0.0064 0.0129
STEP 69-70 growth / min. pool 0.439 / 123.75 0.430 / 123.75 0.440 / 123.75
STEP 73-75 median FVA range (GEM / ec / ecProt) 3.76 / 0.613 / 0.030 not run 3.76 / 0.624 / 0.031 ²
STEP 76-77 growth light / full 0.0280 / 0.0280 not run 0.0280 / 0.0280

¹ The remaining gap is the pool size: MATLAB sets f from PaxDB at STEP 32 (116.5 vs 125 mg/gDCW); scaled to the same pool, growth is identical.
² From the earlier run of the corrected FVA block.

Of the 4834 shared ec.kcat entries, 839 differed from MATLAB before and 127 differ after. The rest are left for follow-up.

Testing

  • New tests for each library fix, and each fails without its fix: test_multi_ec_org_sa_token_wins_over_any_org_kcat_token, test_yaml_gem_keeps_raven_eccodes, and two download_phyl_dist tests.
  • Full suite: 1422 passed, 1 xfailed; ruff check src tests clean.
  • Tutorial run end to end with Gurobi 13 in 13 minutes, exit status 0.

Not addressed here

  • classify_kcat_source (evotune priors) matches ec.source literally, but apply_kcat_list writes BRENDA sources as brenda (wc=0, origin=2). On a geckopy-built model, every BRENDA kcat therefore gets sigma0_log_default instead of the brenda tier. The tutorial TOML's okp group also doesn't match this tutorial's dlkcat kcats. Left alone because fixing it changes evotune results.
  • STEP 32 cannot mirror MATLAB's PaxDB-based f-factor: there is no paxDB.tsv in the tutorial, and calculate_f_factor requires prot_data.
  • Models saved by GECKO 4 MATLAB with RAVEN 3.0.0b1 cannot be read by read_yaml_model. RAVEN's writeYAMLmodel.m writes ec-rxns enzyme entries at the same indentation as their parent key; that fix belongs in RAVEN.

…g MATLAB

Within one EC number the six levels are tried in the order 1, 2, 3, 5,
4, 6 (a correct-organism specific activity before an any-organism
any-substrate kcat), as in fuzzyKcatMatching.m. When a reaction has
several EC numbers, fuzzyKcatMatching.m compares the matches in that
same search order; geckopy compared the output origin numbers instead,
so origin 4 beat origin 5. SEARCH_RANK, derived from SEARCH_LEVELS, now
drives the selection across EC numbers too.

Effect on the yeast-GEM tutorial: r_0005 (EC 2.4.1.34;3.6.5.2) got
0.0005 /s from a human GTPase value for the second EC, where MATLAB
takes 13.3 /s from the S. cerevisiae specific activity for the first;
r_0510 likewise (0.018 vs 603 /s). r_0005 alone held the pre-tuning
growth rate at 0.0013 /hour; with this change it is 0.12 /hour, as in
MATLAB.
cobra.io.load_yaml_model drops RAVEN-specific reaction fields, including
the per-reaction eccodes list that RAVEN YAML uses for EC numbers, so a
yeast-GEM loaded this way had no EC numbers at all and
fill_eccodes_from_gem returned empty codes for every reaction.
read_yaml_model lifts eccodes into annotation['ec-code'] and also reads
plain cobra-style YAML, so it replaces cobra's reader for .yml/.yaml.
fuzzy_kcat_matching needs the KEGG phylogenetic distances to pick the
closest organism when BRENDA has no value for the model's own organism.
MATLAB reads them from RAVEN's keggPhylDist.mat; download_phyl_dist
fetches that file (pinned to RAVEN 3.0.0b1) to a given path, typically
adapter.get_phyl_dist_path(). The file is validated with load_phyl_dist
before it replaces the destination.
…ht vs full

- STEP 16-17: fill_eccodes_from_database gets kegg_db, as MATLAB's
  getECfromDatabase consults KEGG where UniProt has no or only an
  incomplete EC number.
- STEP 18-19: data/PhylDist.mat was a two-organism stub, so the "closest
  organism" levels fell back to arbitrary values (r_0103_REV took
  5e6 /s from Cupriavidus necator where MATLAB takes 3 /s from
  Aspergillus fumigatus). The tutorial now downloads RAVEN's
  keggPhylDist.mat on first run; the file is git-ignored.
- STEP 26: merge_kcats replaces the deprecated
  merge_dlkcat_and_fuzzy_kcats.
- STEP 43-44: the note about the unported ABC-SMC tuning is replaced by
  a pointer to evotune, not run, as in the MATLAB tutorial.
- STEP 73-75: ecFVA runs (skipped on GLPK). The commented-out version
  ran ec_fva on the ecModel for the "bare GEM" column, duplicating the
  ecModel result.
- STEP 76-77: code/plot_light_vs_full.py ports plotlightVSfull.m; light
  and full both reach 0.0280 /hour, as in MATLAB.
- model_adapter.toml: max_growth_weight is applied to the max-growth
  term in both MATLAB and geckopy; the comment said otherwise.

Together with the fuzzy_kcat_matching fix, the tutorial now tracks the
MATLAB protocol: 7 kcats tuned (growth 0.428 vs 0.427 /hour), 0.0129
/hour after flux constraints (MATLAB 0.0129) and 123.75 mg/gDCW minimum
pool usage (identical).
- EvotuneParams.max_growth_weight: MATLAB applies maxGrowthWeight to the
  max-growth term, like geckopy, so no conversion is needed.
- prior_penalty_weight comment: GECKO 4's cmaesKcatTuning has the same
  priorPenaltyWeight term, so reproducing MATLAB does not require 0.
- Migration guide: download_phyl_dist in the kcat example, and YAML
  starting GEMs load through raven-toolbox.
- CHANGELOG: Unreleased section.
The three links to the "MATLAB ↔ Python translation guide" at
gecko-docs.readthedocs.io/en/latest/api/translation/ return 404; that
page does not exist on gecko-docs. The per-function differences they
refer to are on the "GECKO vs. geckopy" page (gecko-to-geckopy.html).
@edkerk
edkerk merged commit 17cc7b0 into develop Sep 11, 2026
5 checks passed
@edkerk
edkerk deleted the fix/fuzzy-origin-rank-and-tutorial branch September 11, 2026 21:09
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant