From 27e39c54943306797baea49e3ba4ffab4e60e2dc Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 23 Sep 2026 01:02:25 +0200 Subject: [PATCH 1/4] Document BRENDA 2026.1 data trap, floor fix, and the DLKcat/floor interaction Investigates every difference between the June-2026 ecYeastGEM build and mainline geckopy today: matcher bug fixes, standard-kcat scope, OKP tiering, BRENDA 2026.1, and prediction source (CatPred/DLKcat/EITLEM). Six changes tune to the same place as June within noise; DLKcat and EITLEM both improve on it. BRENDA 2026.1 breaks tuning outright via a specific-activity/molecular-weight data trap (EC 2.3.1.297), traced to a mouse row surviving a dropped yeast row in the SA-MW join. Two mitigations compared (top_origin_limit=4, and a relative floor on derived kcats); the floor alone gives the best results for three of four predictors but breaks DLKcat via a Crabtree-effect/overflow mismatch, traced by comparing per-condition simulated vs measured ethanol secretion. Combining both mitigations is clean for every predictor tested. --- docs/internal/pipeline_deltas_since_june.md | 173 ++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 docs/internal/pipeline_deltas_since_june.md diff --git a/docs/internal/pipeline_deltas_since_june.md b/docs/internal/pipeline_deltas_since_june.md new file mode 100644 index 0000000..4d635b1 --- /dev/null +++ b/docs/internal/pipeline_deltas_since_june.md @@ -0,0 +1,173 @@ +# What changed in the kcat-assignment pipeline since the June model, and whether it affects tuning + +The ecYeastGEM model tuned in June (`ecYeastGEM_preTune_GECKOderived.yml`) was +built by a snapshot of geckopy that has since moved on. This investigates +every difference between that snapshot and mainline today, one change at a +time, and checks whether CMA-ES tuning (top-112 free parameters, 3 seeds, +budget 6800 evaluations, `max_growth_weight=2`) still lands at the same +place. The June model tunes to a distance of 0.798 ± 0.004 (mean ± sd over +3 seeds) — every number below is judged against that noise floor. + +Everything here is reproduced by re-running the assignment stage (fuzzy +BRENDA matching through `apply_kcat_constraints`) from raw inputs, isolated +one change at a time, on the June ecModel's own reaction set. A byte-exact +replay of every June bug first reproduced the June kcats 100% identically, +so each change below starts from a verified-identical baseline. + +## Summary + +Six pipeline changes tune to the same place as June, within noise: + +| Change | Tuned distance | vs June | +|---|---|---| +| source labels only | 0.798 | +0.0% | +| substring-wildcard fix | 0.796 | −0.1% | +| tie-break fix | 0.792 | −0.7% | +| min-coefficient fix | 0.791 | −0.9% | +| standard kcat (subsystem, not model-wide) | 0.797 | −0.1% | +| OKP predictions in real tiers (not one bloc) | 0.796 | −0.2% | +| all of the above together | 0.794 | −0.5% | + +One change is a clear improvement: + +| Change | Tuned distance | vs June | +|---|---|---| +| DLKcat instead of CatPred | 0.671 | **−15.8%** | +| EITLEM instead of CatPred | 0.769 | −3.6% | + +One change breaks tuning outright, traced to a specific data-quality issue, +with two candidate fixes: + +| Change | Tuned distance | vs June | +|---|---|---| +| BRENDA 2026.1 (unmitigated) | 4.19 | broken | +| BRENDA 2026.1, `top_origin_limit=4` | 0.805–0.850 | same, or +6.5% for the unmitigated combination | +| BRENDA 2026.1, floor filter | 0.731–0.781 | **better**, except paired with DLKcat | +| BRENDA 2026.1, floor + `top_origin_limit=4` | 0.701–0.764 | same or better, every predictor | + +## The BRENDA 2026.1 data trap + +BRENDA's kcat table has one row for EC 2.3.1.297 at a specific activity of +6×10⁻⁷ µmol/min/mg for *Mus musculus*, joined against a molecular weight of +42,000 g/mol — almost certainly a unit or transcription error in the +source data, giving a derived kcat of 4.2×10⁻⁷ 1/s. BRENDA also has a +*Saccharomyces cerevisiae* specific-activity row for the same EC (20.68 +µmol/min/mg, sensible), but no yeast molecular weight for that EC, so the +loader's SA×MW join drops it — the mouse row survives on its own and wins. + +At the mainline default `top_origin_limit=6`, any exact-EC BRENDA match, +including specific-activity-derived ones, outranks a prediction. The +4.2×10⁻⁷ kcat gets assigned, the resulting enzyme-usage coefficient +(MW/(kcat·3600)) is enormous, and CMA-ES cannot recover a sane model in +budget: distance 4.19–4.38 depending on aggregation, vs 0.80 for June. + +This is not unique to this one EC. Across the whole BRENDA 2026.1 dump, +5,363 of 17,829 specific-activity rows (30%) have no same-organism +molecular weight and are silently dropped; 608 kcat-table rows and 1,913 +SA-derived rows sit below 10⁻² 1/s, some of which are likely errors of +the same kind. + +### Two candidate fixes + +**`top_origin_limit=4`** (existing parameter, no code change): demotes +specific-activity-derived matches (origins 5–6) below predictions, so a +prediction is used instead of the mouse row. Fixes the break +(`top_origin_limit=6`: 4.19 → `=4`: 0.81, close to June) but does cost a +little on the unmitigated pipeline (`b26_o4`: +6.5%) and is redundant with +the fix below once that's applied. + +**A relative floor on derived kcats** (`brenda_fix.py`, sandbox-only, not +yet in mainline): drop a specific-activity-derived kcat if it's more than +4 orders of magnitude below the median of the same EC family (widening to +the EC's 3-, 2-, then 1-digit prefix if fewer than 5 reference values are +available). This is the more direct fix — it removes the bad row instead +of demoting a whole tier — and gives the best results seen in this +investigation when combined with `top_origin_limit=4`: + +| Predictor | floor alone (limit=6) | floor + limit=4 | +|---|---|---| +| CatPred | 0.731 (−8.3%) | — | +| EITLEM | 0.718 (−10.0%) | 0.701–0.826 | +| unmitigated pipeline | 0.781 (−2.1%) | — | +| DLKcat | **1.80, broken** | **0.72–0.76, fixed** | + +An earlier version of the floor also carried a molecular-weight fallback +(use the EC family's median MW when the organism has no MW row, recovering +the yeast SA value directly instead of relying on the floor to catch the +alternative). Tested and rejected: combining it with the floor was worse +than either fix alone, and broke the DLKcat pairing (1.13, and 10/112 top +parameters overlapping with June, vs 96+/112 for clean variants) the same +way the unmitigated data does. The floor's own reference pool is built +from the union of the kcat and SA-derived tables; recovering more rows +into that pool before filtering can quietly lower a family's "typical" +value and let a bad entry back in. Removed; only the floor remains. + +### Why the floor alone broke DLKcat, and why `top_origin_limit=4` doesn't + +The floor removes the single obviously-bad row (EC 2.3.1.297), but not +every specific-activity-derived kcat that passes the floor is a *good* +value — some are real BRENDA measurements that are simply less reliable +than a direct kcat. At `top_origin_limit=6`, 365 reactions swap from +whatever the predictor said to one of these floor-surviving +specific-activity values — the same 365 reactions, with the same new +values, regardless of which predictor (CatPred or DLKcat) is paired with +BRENDA, since the swap only depends on BRENDA's own matches. + +Tuned against CatPred, this is fine (`c3_fl`: 0.731). Tuned against +DLKcat, it isn't (`c5_fl`: 1.80). Tracing this by comparing per-condition +RMSE between the two: the failure concentrates in 7 of 33 flux conditions, +all moderate-to-high growth rate on glucose — exactly where *S. cerevisiae* +should switch on the Crabtree effect and secrete ethanol even aerobically, +because respiratory capacity becomes limiting. The DLKcat-tuned model +systematically **under-produces ethanol and over-produces growth** for the +same glucose uptake (e.g. one condition: measured ethanol 9.5 mmol/gDW/h, +simulated 0.0; growth 0.362 vs the CatPred-tuned model's 0.338, closer to +measured). No enzyme is anywhere near its usage cap in these conditions +(all under 2%), so this isn't a protein-budget bottleneck — it's a missed +metabolic switch. DLKcat's own predictions for the reactions it keeps (the +~700+ reactions not part of the 365-reaction swap) apparently combine with +the floor-admitted BRENDA values to make "too respiratory, not enough +overflow" the solution CMA-ES converges to within budget; CatPred's +landscape doesn't have that trap. `top_origin_limit=4` avoids it by +additionally demoting the floor-surviving specific-activity tier below +*any* prediction, so DLKcat's own values are used more broadly instead. + +The exact reaction(s) responsible for the respiration/overflow trade-off +were not pinned down further (would mean comparing the fermentation +pathway — pyruvate decarboxylase, alcohol dehydrogenase — and the +respiratory chain between the working and broken variants for the largest +kcat swaps). The phenomenon and its trigger (floor + DLKcat, without the +origin-limit demotion) are solid; the single-reaction mechanism isn't +required to know **floor + `top_origin_limit=4` together is safe for every +predictor tested**, which is the practical recommendation below. + +## Recommendation + +- **Commit**: nothing yet. The floor fix (`brenda_fix.py`) is sandbox-only + and needs porting into `geckopy.databases.brenda_loader` before it's a + real candidate; it should ship together with `top_origin_limit=4` + becoming the default, not alone, given the DLKcat interaction above. +- **Do not commit**: the molecular-weight fallback. Tested and rejected. +- **Document only**: BRENDA 2026.1's specific-activity-derived tail has a + genuine outlier problem inherited from the source data, independent of + this pipeline. EITLEM (already-supported as an OKP method, just a + different `method=` argument) gives a modest, consistent improvement + over CatPred and is worth considering as an alternative default + predictor; DLKcat's larger improvement is real but its parameter + selection overlaps least with June's (67–83 of 112 top parameters + in common, lowest of any variant tested), meaning the two predictors + are tuning meaningfully different aspects of the model, not just + reproducing the same fit with different labels. + +## Method + +Reassignment driven by `deltas/pipeline.py` (`Cfg` dataclass toggles each +change), using a private copy of the matcher (`deltas/fuzzy_compat/`) for +the three matcher-bug switches and the real mainline `merge_kcats`, +`assign_standard_kcat`, `apply_kcat_constraints` for everything else. +Screening and CMA-ES: `deltas/variant_tune.py`, same procedure as the +June-reproduction runs (screen every free parameter ×2/÷2 with tied +isozyme copies moved as a unit, top-112 by leverage×sigma0, separable +CMA-ES, 3 seeds, budget 6800, `max_growth_weight=2`). Per-condition +decomposition for the DLKcat trace: `deltas/decompose.py`, +`deltas/enzyme_trace2.py`. From 6c67401c7741e69347bd1579a46007fd6e96ec5a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 23 Sep 2026 08:47:01 +0200 Subject: [PATCH 2/4] Note that floor+limit4 is redundant with limit4 alone for CatPred c3_fl4 vs c3_o4 differ in only 34 of 4834 kcats; b26_fl4 vs b26_o4 are byte-identical. Not a bug: top_origin_limit=4 already demotes every origin-5/6 BRENDA match below any prediction, and CatPred's coverage is broad enough that almost every such reaction already has a CatPred prediction to fall back to. The floor only has room to matter where a predictor's coverage is patchier (DLKcat, EITLEM). --- docs/internal/pipeline_deltas_since_june.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/internal/pipeline_deltas_since_june.md b/docs/internal/pipeline_deltas_since_june.md index 4d635b1..353e40a 100644 --- a/docs/internal/pipeline_deltas_since_june.md +++ b/docs/internal/pipeline_deltas_since_june.md @@ -141,6 +141,16 @@ origin-limit demotion) are solid; the single-reaction mechanism isn't required to know **floor + `top_origin_limit=4` together is safe for every predictor tested**, which is the practical recommendation below. +For CatPred specifically, combining the two mitigations makes almost no +difference versus `top_origin_limit=4` alone (`c3_fl4` vs `c3_o4`: 34 of +4,834 kcats differ; `b26_fl4` vs `b26_o4`: byte-identical). This isn't a +bug — `top_origin_limit=4` already demotes every origin-5/6 BRENDA match +below *any* prediction, and CatPred's coverage is broad enough that almost +every reaction with an origin-5/6 match also has a CatPred prediction +available to fall back to. The floor only has room to matter where a +predictor's coverage is patchier (DLKcat, EITLEM) and a reaction's only +alternative to the bad BRENDA value is standard kcat, not a prediction. + ## Recommendation - **Commit**: nothing yet. The floor fix (`brenda_fix.py`) is sandbox-only From 79b75a027d7c002b7a7c87b0650511d0bcc032d5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 23 Sep 2026 09:51:54 +0200 Subject: [PATCH 3/4] Port the BRENDA outlier floor and default top_origin_limit=4 Two changes, meant to ship together (see docs/internal/pipeline_deltas_since_june.md): - load_brenda_data gains filter_outlier_kcats (default True): drops a kcat more than 4 orders of magnitude below the median of its own EC family (widening to a broader EC prefix when fewer than 5 reference values are available), computed from the union of the plain-kcat and SA-derived tables. Catches cases like EC 2.3.1.297, where the only molecular-weight row available is for the wrong organism and the resulting SA x MW join gives a kcat 8 orders of magnitude below a same-EC measurement that has no MW row and would otherwise be silently dropped instead of used. - merge_kcats' and merge_dlkcat_and_fuzzy_kcats' top_origin_limit default changes from 6 to 4, so a specific-activity-derived BRENDA match (origin 5-6) no longer automatically outranks a prediction. Redundant with the filter above when a predictor's coverage is broad (e.g. CatPred), but still needed on its own where coverage is patchier (DLKcat, EITLEM) -- see the doc for the Crabtree-effect mismatch this prevents when paired with the floor alone. Both changes are opt-out (filter_outlier_kcats=False, top_origin_limit=6 explicit) rather than removed, so the pre-existing behaviour stays reachable. Full test suite passes (1435 passed, 1 xfailed, 0 failed); added dedicated tests for the floor filter and updated the one test whose name and assertion encoded the old top_origin_limit default as intentional behaviour. --- src/geckopy/databases/brenda_loader.py | 108 +++++++++++++++++- .../merge_dlkcat_and_fuzzy_kcats.py | 2 +- src/geckopy/gather_kcats/merge_kcats.py | 9 +- tests/test_brenda_loader.py | 79 +++++++++++++ tests/test_merge_dlkcat_and_fuzzy_kcats.py | 21 +++- 5 files changed, 214 insertions(+), 5 deletions(-) diff --git a/src/geckopy/databases/brenda_loader.py b/src/geckopy/databases/brenda_loader.py index 4fe7827..0e005d0 100644 --- a/src/geckopy/databases/brenda_loader.py +++ b/src/geckopy/databases/brenda_loader.py @@ -33,6 +33,7 @@ from pathlib import Path from typing import IO, Literal +import numpy as np import pandas as pd logger = logging.getLogger(__name__) @@ -95,7 +96,9 @@ def sa_for(self, aggregation: Literal["max", "median"]) -> pd.DataFrame: ) -def load_brenda_data(folder: str | Path) -> BrendaData: +def load_brenda_data( + folder: str | Path, *, filter_outlier_kcats: bool = True, +) -> BrendaData: """Load BRENDA kcat, SA, and MW dumps from a folder. Ported from GECKO MATLAB: @@ -129,6 +132,12 @@ def load_brenda_data(folder: str | Path) -> BrendaData: ---------- folder Directory containing all three BRENDA dump files. + filter_outlier_kcats + If True (default), drop a row whose kcat is more than 4 orders + of magnitude below the typical kcat of its EC before it can be + matched to any reaction -- see ``_filter_outlier_rows`` for the + rationale and mechanism. Set False to reproduce the unfiltered + 1.0-era behaviour. Returns ------- @@ -170,6 +179,11 @@ def load_brenda_data(folder: str | Path) -> BrendaData: _sa_view_from_wide(sa_wide, "sa_median"), mw_raw, ) + if filter_outlier_kcats: + kcat_max, kcat_median, sa_max, sa_median = _filter_outlier_rows( + [kcat_max, kcat_median, sa_max, sa_median], + ) + return BrendaData( kcat_max=kcat_max, kcat_median=kcat_median, @@ -371,3 +385,95 @@ def _join_sa_with_mw(sa: pd.DataFrame, mw: pd.DataFrame) -> pd.DataFrame: if not out_rows: return pd.DataFrame(columns=_SA_COLUMNS) return pd.DataFrame(out_rows, columns=_SA_COLUMNS) + + +def _filter_outlier_rows( + tables: list[pd.DataFrame], *, decades: float = 4.0, min_ref: int = 5, +) -> list[pd.DataFrame]: + """Drop a row whose kcat is implausibly low for its EC. + + BRENDA's specific-activity x molecular-weight join can produce a + derived kcat that is many orders of magnitude below anything else + reported for the same EC -- almost always a unit or transcription + error upstream, not a real slow enzyme (a documented case: EC + 2.3.1.297's only molecular-weight row is for a different organism + than its lowest specific-activity row, and the derived kcat this + produces, ~4e-7 1/s, is 8 orders of magnitude below the yeast SA + measurement for the same EC that has no matching MW row and is + silently dropped instead). Once matched to a reaction, a kcat that + small makes the enzyme-usage coefficient enormous and can make a + model untunable. + + This is a *relative* floor, not a fixed threshold: a row is dropped + only if it is more than ``decades`` orders of magnitude below the + median kcat of its own EC. "Its own EC" widens from the full + 4-level code to the 3-, 2-, then 1-level prefix if fewer than + ``min_ref`` values are available at the current level, so a + sparsely-represented EC is still judged against something. The + reference pool for one EC is built from every positive-kcat row of + every table passed in (plain kcat rows and SA-derived rows alike), + computed once up front from the *unfiltered* tables so that + dropping a row in one table cannot shift the reference used to + judge another. + + Real enzymes can be genuinely slow, so this only catches a row far + outside its own EC's spread, not a merely-low one. + + Parameters + ---------- + tables + The four BRENDA views (``kcat_max``, ``kcat_median``, + ``sa_max``, ``sa_median``); each must have ``ec_code`` and + ``kcat`` columns. Order is preserved in the return value. + decades + How many orders of magnitude below the reference counts as an + outlier. + min_ref + Minimum number of reference values required before a level is + trusted; the search widens to a broader EC prefix otherwise. + + Returns + ------- + list[pandas.DataFrame] + The same tables, each with outlier rows removed and its index + reset. + """ + pool: dict[str, list[float]] = {} + for table in tables: + for ec, kcat in zip(table["ec_code"], table["kcat"]): + if kcat > 0 and "-" not in ec: + pool.setdefault(ec, []).append(np.log10(kcat)) + + log_values_by_level: dict[tuple[int, str], list[float]] = {} + for ec, values in pool.items(): + parts = ec.split(".") + for level in (4, 3, 2, 1): + key = (level, ".".join(parts[:level])) + log_values_by_level.setdefault(key, []).extend(values) + + reference_cache: dict[str, float | None] = {} + + def reference_log_kcat(ec: str) -> float | None: + if ec not in reference_cache: + parts = ec.split(".") + reference_cache[ec] = None + for level in (4, 3, 2, 1): + values = log_values_by_level.get((level, ".".join(parts[:level])), []) + if len(values) >= min_ref: + reference_cache[ec] = float(np.median(values)) + break + return reference_cache[ec] + + filtered = [] + for table in tables: + if table.empty: + # An empty boolean mask loses the (already-correct) column + # set on an empty frame; nothing to filter anyway. + filtered.append(table) + continue + keep = [] + for ec, kcat in zip(table["ec_code"], table["kcat"]): + ref = reference_log_kcat(ec) if kcat > 0 else None + keep.append(not (ref is not None and np.log10(kcat) < ref - decades)) + filtered.append(table[pd.Series(keep, index=table.index)].reset_index(drop=True)) + return filtered diff --git a/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py b/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py index 7455d5c..d88db2e 100644 --- a/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py +++ b/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py @@ -22,7 +22,7 @@ def merge_dlkcat_and_fuzzy_kcats( dlkcat_kcats: pd.DataFrame, fuzzy_kcats: pd.DataFrame, *, - top_origin_limit: int = 6, + top_origin_limit: int = 4, bottom_origin_limit: int = 6, wildcard_limit: int = 3, ) -> pd.DataFrame: diff --git a/src/geckopy/gather_kcats/merge_kcats.py b/src/geckopy/gather_kcats/merge_kcats.py index 1443fca..860a89d 100644 --- a/src/geckopy/gather_kcats/merge_kcats.py +++ b/src/geckopy/gather_kcats/merge_kcats.py @@ -94,7 +94,7 @@ def _strip_okp_prefix(normalized_source: str) -> str: def merge_kcats( *kcat_lists: pd.DataFrame, source_priority: Sequence[str], - top_origin_limit: int = 6, + top_origin_limit: int = 4, bottom_origin_limit: int = 6, wildcard_limit: int = 3, database_sources: Iterable[str] = _DEFAULT_DATABASE_SOURCES, @@ -126,6 +126,13 @@ def merge_kcats( source not listed here is dropped (with a warning). top_origin_limit Origin upper bound for ``database_top`` fuzzy rows. ``[1, 6]``. + Default ``4``, not ``6``: origins 5-6 are specific-activity × + molecular-weight derived kcats, which BRENDA's own data quality + does not always support treating as better than a prediction + (a bad SA/MW join can otherwise silently override a good + prediction with an outlier -- see ``load_brenda_data``'s + ``filter_outlier_kcats``, which removes the worst of these + before they reach this function, but does not catch every one). bottom_origin_limit Origin upper bound for ``database_bottom`` fuzzy rows. ``[1, 6]``. wildcard_limit diff --git a/tests/test_brenda_loader.py b/tests/test_brenda_loader.py index 7585f24..266e874 100644 --- a/tests/test_brenda_loader.py +++ b/tests/test_brenda_loader.py @@ -88,6 +88,85 @@ def test_all_empty_files_yield_empty_dataframes(tmp_path): assert result.sa_median.empty +# --------------------------------------------------------------------------- # +# Outlier filtering (filter_outlier_kcats) +# --------------------------------------------------------------------------- # + + +def _kcat_row(ec, organism, kcat): + return f"{ec}\tm1\t{organism}\t{kcat}\t{kcat}\t1\t*\n" + + +def test_outlier_row_dropped_by_default(tmp_path): + # Five real-looking values for the same EC (the reference pool), + # plus one 8-decades-lower outlier for that same EC. + kcat = "".join( + _kcat_row("1.1.1.1", f"org{i}", v) + for i, v in enumerate([10, 12, 8, 15, 11]) + ) + kcat += _kcat_row("1.1.1.1", "org_outlier", 1e-6) + _write_brenda_files(tmp_path, kcat=kcat) + + filtered = load_brenda_data(tmp_path) + assert len(filtered.kcat_max) == 5 + assert 1e-6 not in filtered.kcat_max["kcat"].values + + unfiltered = load_brenda_data(tmp_path, filter_outlier_kcats=False) + assert len(unfiltered.kcat_max) == 6 + assert 1e-6 in unfiltered.kcat_max["kcat"].values + + +def test_reference_pool_widens_to_ec_prefix_when_sparse(tmp_path): + # Only one row for 1.1.1.9 itself (too few for its own reference), + # but four more under the shared 1.1.1.- prefix -- enough (5) once + # widened. The lone 1.1.1.9 row is a real outlier against that + # widened reference and should be dropped. + kcat = "".join( + _kcat_row(ec, f"org{i}", v) + for i, (ec, v) in enumerate([ + ("1.1.1.1", 10), ("1.1.1.2", 12), ("1.1.1.3", 8), ("1.1.1.4", 15), + ]) + ) + kcat += _kcat_row("1.1.1.9", "org_outlier", 1e-6) + _write_brenda_files(tmp_path, kcat=kcat) + + filtered = load_brenda_data(tmp_path) + assert "1.1.1.9" not in filtered.kcat_max["ec_code"].values + assert len(filtered.kcat_max) == 4 + + +def test_outlier_filter_uses_kcat_and_sa_together_as_one_reference_pool(tmp_path): + # The reference pool for one EC draws on both the plain kcat table + # and the SA-derived table; four kcat rows plus one SA row is + # enough (5) to judge a sixth, wildly-off SA row as an outlier. + kcat = "".join( + _kcat_row("2.3.1.297", f"org{i}", v) + for i, v in enumerate([10, 12, 8, 15]) + ) + sa = ( + "2.3.1.297\t*\torg_good\t2.0\t2.0\t1\t*\n" # ~ SA*MW = 20, in range + "2.3.1.297\t*\torg_bad\t1e-8\t1e-8\t1\t*\n" # ~ SA*MW = 1e-7, outlier + ) + mw = ( + "2.3.1.297\t*\torg_good\t10000\t1\t*\n" + "2.3.1.297\t*\torg_bad\t10000\t1\t*\n" + ) + _write_brenda_files(tmp_path, kcat=kcat, sa=sa, mw=mw) + + filtered = load_brenda_data(tmp_path) + assert len(filtered.sa_max) == 1 + assert filtered.sa_max.iloc[0]["organism"] == "org_good" + + +def test_empty_tables_unaffected_by_outlier_filter(tmp_path): + _write_brenda_files(tmp_path) + result = load_brenda_data(tmp_path, filter_outlier_kcats=True) + assert list(result.kcat_max.columns) == [ + "ec_code", "substrate", "organism", "kcat", "n", + ] + assert result.kcat_max.empty + + # --------------------------------------------------------------------------- # # Missing files # --------------------------------------------------------------------------- # diff --git a/tests/test_merge_dlkcat_and_fuzzy_kcats.py b/tests/test_merge_dlkcat_and_fuzzy_kcats.py index f083bde..cd52cc7 100644 --- a/tests/test_merge_dlkcat_and_fuzzy_kcats.py +++ b/tests/test_merge_dlkcat_and_fuzzy_kcats.py @@ -160,13 +160,30 @@ def test_fuzzy_prio1_wins_when_reaction_also_in_dlkcat(): assert row["kcat"] == 5.0 -def test_fuzzy_prio1_uses_default_top_origin_limit_6(): - """With default top_origin_limit=6, origin=6 still qualifies for prio1.""" +def test_fuzzy_prio1_uses_default_top_origin_limit_4(): + """Default top_origin_limit=4: an origin=6 (specific-activity x MW + derived) match does NOT qualify for prio1, so DLKcat wins instead. + Origins 5-6 are demoted by default because BRENDA's SA/MW join can + produce a wildly wrong outlier that would otherwise silently + outrank a good prediction -- see brenda_loader's + filter_outlier_kcats, which this default is meant to pair with.""" out = merge_dlkcat_and_fuzzy_kcats( _dlkcat_df([_dlkcat_row("r1", 100.0)]), _fuzzy_df([_fuzzy_row("r1", 5.0, wildcard_level=0, origin=6)]), ) assert len(out) == 1 + assert out.iloc[0]["source"] == "DLKcat" + + +def test_fuzzy_prio1_origin_6_qualifies_when_limit_explicitly_widened(): + """The pre-4-default behaviour is still reachable by passing + top_origin_limit=6 explicitly.""" + out = merge_dlkcat_and_fuzzy_kcats( + _dlkcat_df([_dlkcat_row("r1", 100.0)]), + _fuzzy_df([_fuzzy_row("r1", 5.0, wildcard_level=0, origin=6)]), + top_origin_limit=6, + ) + assert len(out) == 1 assert out.iloc[0]["source"] == "brenda" From 7e325b28cb4c89b19f6ffe155bea5314fcbe2af1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 23 Sep 2026 17:57:32 +0200 Subject: [PATCH 4/4] Drop the internal pipeline-deltas writeup from this PR Investigation notes, not part of the public change. --- docs/internal/pipeline_deltas_since_june.md | 183 -------------------- 1 file changed, 183 deletions(-) delete mode 100644 docs/internal/pipeline_deltas_since_june.md diff --git a/docs/internal/pipeline_deltas_since_june.md b/docs/internal/pipeline_deltas_since_june.md deleted file mode 100644 index 353e40a..0000000 --- a/docs/internal/pipeline_deltas_since_june.md +++ /dev/null @@ -1,183 +0,0 @@ -# What changed in the kcat-assignment pipeline since the June model, and whether it affects tuning - -The ecYeastGEM model tuned in June (`ecYeastGEM_preTune_GECKOderived.yml`) was -built by a snapshot of geckopy that has since moved on. This investigates -every difference between that snapshot and mainline today, one change at a -time, and checks whether CMA-ES tuning (top-112 free parameters, 3 seeds, -budget 6800 evaluations, `max_growth_weight=2`) still lands at the same -place. The June model tunes to a distance of 0.798 ± 0.004 (mean ± sd over -3 seeds) — every number below is judged against that noise floor. - -Everything here is reproduced by re-running the assignment stage (fuzzy -BRENDA matching through `apply_kcat_constraints`) from raw inputs, isolated -one change at a time, on the June ecModel's own reaction set. A byte-exact -replay of every June bug first reproduced the June kcats 100% identically, -so each change below starts from a verified-identical baseline. - -## Summary - -Six pipeline changes tune to the same place as June, within noise: - -| Change | Tuned distance | vs June | -|---|---|---| -| source labels only | 0.798 | +0.0% | -| substring-wildcard fix | 0.796 | −0.1% | -| tie-break fix | 0.792 | −0.7% | -| min-coefficient fix | 0.791 | −0.9% | -| standard kcat (subsystem, not model-wide) | 0.797 | −0.1% | -| OKP predictions in real tiers (not one bloc) | 0.796 | −0.2% | -| all of the above together | 0.794 | −0.5% | - -One change is a clear improvement: - -| Change | Tuned distance | vs June | -|---|---|---| -| DLKcat instead of CatPred | 0.671 | **−15.8%** | -| EITLEM instead of CatPred | 0.769 | −3.6% | - -One change breaks tuning outright, traced to a specific data-quality issue, -with two candidate fixes: - -| Change | Tuned distance | vs June | -|---|---|---| -| BRENDA 2026.1 (unmitigated) | 4.19 | broken | -| BRENDA 2026.1, `top_origin_limit=4` | 0.805–0.850 | same, or +6.5% for the unmitigated combination | -| BRENDA 2026.1, floor filter | 0.731–0.781 | **better**, except paired with DLKcat | -| BRENDA 2026.1, floor + `top_origin_limit=4` | 0.701–0.764 | same or better, every predictor | - -## The BRENDA 2026.1 data trap - -BRENDA's kcat table has one row for EC 2.3.1.297 at a specific activity of -6×10⁻⁷ µmol/min/mg for *Mus musculus*, joined against a molecular weight of -42,000 g/mol — almost certainly a unit or transcription error in the -source data, giving a derived kcat of 4.2×10⁻⁷ 1/s. BRENDA also has a -*Saccharomyces cerevisiae* specific-activity row for the same EC (20.68 -µmol/min/mg, sensible), but no yeast molecular weight for that EC, so the -loader's SA×MW join drops it — the mouse row survives on its own and wins. - -At the mainline default `top_origin_limit=6`, any exact-EC BRENDA match, -including specific-activity-derived ones, outranks a prediction. The -4.2×10⁻⁷ kcat gets assigned, the resulting enzyme-usage coefficient -(MW/(kcat·3600)) is enormous, and CMA-ES cannot recover a sane model in -budget: distance 4.19–4.38 depending on aggregation, vs 0.80 for June. - -This is not unique to this one EC. Across the whole BRENDA 2026.1 dump, -5,363 of 17,829 specific-activity rows (30%) have no same-organism -molecular weight and are silently dropped; 608 kcat-table rows and 1,913 -SA-derived rows sit below 10⁻² 1/s, some of which are likely errors of -the same kind. - -### Two candidate fixes - -**`top_origin_limit=4`** (existing parameter, no code change): demotes -specific-activity-derived matches (origins 5–6) below predictions, so a -prediction is used instead of the mouse row. Fixes the break -(`top_origin_limit=6`: 4.19 → `=4`: 0.81, close to June) but does cost a -little on the unmitigated pipeline (`b26_o4`: +6.5%) and is redundant with -the fix below once that's applied. - -**A relative floor on derived kcats** (`brenda_fix.py`, sandbox-only, not -yet in mainline): drop a specific-activity-derived kcat if it's more than -4 orders of magnitude below the median of the same EC family (widening to -the EC's 3-, 2-, then 1-digit prefix if fewer than 5 reference values are -available). This is the more direct fix — it removes the bad row instead -of demoting a whole tier — and gives the best results seen in this -investigation when combined with `top_origin_limit=4`: - -| Predictor | floor alone (limit=6) | floor + limit=4 | -|---|---|---| -| CatPred | 0.731 (−8.3%) | — | -| EITLEM | 0.718 (−10.0%) | 0.701–0.826 | -| unmitigated pipeline | 0.781 (−2.1%) | — | -| DLKcat | **1.80, broken** | **0.72–0.76, fixed** | - -An earlier version of the floor also carried a molecular-weight fallback -(use the EC family's median MW when the organism has no MW row, recovering -the yeast SA value directly instead of relying on the floor to catch the -alternative). Tested and rejected: combining it with the floor was worse -than either fix alone, and broke the DLKcat pairing (1.13, and 10/112 top -parameters overlapping with June, vs 96+/112 for clean variants) the same -way the unmitigated data does. The floor's own reference pool is built -from the union of the kcat and SA-derived tables; recovering more rows -into that pool before filtering can quietly lower a family's "typical" -value and let a bad entry back in. Removed; only the floor remains. - -### Why the floor alone broke DLKcat, and why `top_origin_limit=4` doesn't - -The floor removes the single obviously-bad row (EC 2.3.1.297), but not -every specific-activity-derived kcat that passes the floor is a *good* -value — some are real BRENDA measurements that are simply less reliable -than a direct kcat. At `top_origin_limit=6`, 365 reactions swap from -whatever the predictor said to one of these floor-surviving -specific-activity values — the same 365 reactions, with the same new -values, regardless of which predictor (CatPred or DLKcat) is paired with -BRENDA, since the swap only depends on BRENDA's own matches. - -Tuned against CatPred, this is fine (`c3_fl`: 0.731). Tuned against -DLKcat, it isn't (`c5_fl`: 1.80). Tracing this by comparing per-condition -RMSE between the two: the failure concentrates in 7 of 33 flux conditions, -all moderate-to-high growth rate on glucose — exactly where *S. cerevisiae* -should switch on the Crabtree effect and secrete ethanol even aerobically, -because respiratory capacity becomes limiting. The DLKcat-tuned model -systematically **under-produces ethanol and over-produces growth** for the -same glucose uptake (e.g. one condition: measured ethanol 9.5 mmol/gDW/h, -simulated 0.0; growth 0.362 vs the CatPred-tuned model's 0.338, closer to -measured). No enzyme is anywhere near its usage cap in these conditions -(all under 2%), so this isn't a protein-budget bottleneck — it's a missed -metabolic switch. DLKcat's own predictions for the reactions it keeps (the -~700+ reactions not part of the 365-reaction swap) apparently combine with -the floor-admitted BRENDA values to make "too respiratory, not enough -overflow" the solution CMA-ES converges to within budget; CatPred's -landscape doesn't have that trap. `top_origin_limit=4` avoids it by -additionally demoting the floor-surviving specific-activity tier below -*any* prediction, so DLKcat's own values are used more broadly instead. - -The exact reaction(s) responsible for the respiration/overflow trade-off -were not pinned down further (would mean comparing the fermentation -pathway — pyruvate decarboxylase, alcohol dehydrogenase — and the -respiratory chain between the working and broken variants for the largest -kcat swaps). The phenomenon and its trigger (floor + DLKcat, without the -origin-limit demotion) are solid; the single-reaction mechanism isn't -required to know **floor + `top_origin_limit=4` together is safe for every -predictor tested**, which is the practical recommendation below. - -For CatPred specifically, combining the two mitigations makes almost no -difference versus `top_origin_limit=4` alone (`c3_fl4` vs `c3_o4`: 34 of -4,834 kcats differ; `b26_fl4` vs `b26_o4`: byte-identical). This isn't a -bug — `top_origin_limit=4` already demotes every origin-5/6 BRENDA match -below *any* prediction, and CatPred's coverage is broad enough that almost -every reaction with an origin-5/6 match also has a CatPred prediction -available to fall back to. The floor only has room to matter where a -predictor's coverage is patchier (DLKcat, EITLEM) and a reaction's only -alternative to the bad BRENDA value is standard kcat, not a prediction. - -## Recommendation - -- **Commit**: nothing yet. The floor fix (`brenda_fix.py`) is sandbox-only - and needs porting into `geckopy.databases.brenda_loader` before it's a - real candidate; it should ship together with `top_origin_limit=4` - becoming the default, not alone, given the DLKcat interaction above. -- **Do not commit**: the molecular-weight fallback. Tested and rejected. -- **Document only**: BRENDA 2026.1's specific-activity-derived tail has a - genuine outlier problem inherited from the source data, independent of - this pipeline. EITLEM (already-supported as an OKP method, just a - different `method=` argument) gives a modest, consistent improvement - over CatPred and is worth considering as an alternative default - predictor; DLKcat's larger improvement is real but its parameter - selection overlaps least with June's (67–83 of 112 top parameters - in common, lowest of any variant tested), meaning the two predictors - are tuning meaningfully different aspects of the model, not just - reproducing the same fit with different labels. - -## Method - -Reassignment driven by `deltas/pipeline.py` (`Cfg` dataclass toggles each -change), using a private copy of the matcher (`deltas/fuzzy_compat/`) for -the three matcher-bug switches and the real mainline `merge_kcats`, -`assign_standard_kcat`, `apply_kcat_constraints` for everything else. -Screening and CMA-ES: `deltas/variant_tune.py`, same procedure as the -June-reproduction runs (screen every free parameter ×2/÷2 with tied -isozyme copies moved as a unit, top-112 by leverage×sigma0, separable -CMA-ES, 3 seeds, budget 6800, `max_growth_weight=2`). Per-condition -decomposition for the DLKcat trace: `deltas/decompose.py`, -`deltas/enzyme_trace2.py`.