fix: give the same result as GECKO where the parity scenarios found differences

CHANGELOG.md

@@ -2,6 +2,68 @@
 
 ## Unreleased
 
+* `apply_kcat_constraints` treats an empty `update_rxns` as all reactions, as MATLAB's
+  `applyKcatConstraints` does with an empty `updateRxns`; it used to do nothing. The callers
+  that build the list from a filter (`set_kcat_for_reactions`, the evotune scoring and reset)
+  skip the call when the list is empty.
+* `apply_custom_kcats` reads a reaction list in double quotes (`"R3,R5"`) without the quotes,
+  as MATLAB's `%q` does, so the quoted list matches R3 and R5 instead of nothing.
+* `make_ec_model` refuses what GECKO's `makeEcModel` refuses, with a `ValueError`: a reaction ID
+  that starts with `usage_prot_` or `prot_pool` or ends in `_REV` or `_EXP_<digits>`, a metabolite
+  ID that starts with `prot_`, an `enzyme_comp` that is not one of the compartment names of the
+  model (a compartment ID that differs from its name no longer passes), and a model without
+  genes. Only a reaction with `lb < 0 < ub` that is not an exchange reaction is split into a
+  `_REV` pair: a reaction with both bounds negative stays as it is instead of raising, and a
+  reaction that only consumes or only produces metabolites, such as `A + B <=>`, is an exchange
+  reaction, as in RAVEN's `getExchangeRxns`.
+* `set_kcat_for_reactions` follows `setKcatForReactions`. `apply` defaults to `False`, so only
+  `ec.kcat` and `ec.source` change; pass `apply=True` to write the S matrix. IDs match without
+  regard to letter case, an ID that matches nothing is ignored, and in a gecko-light model the
+  `###_` prefix is part of the ID, so `"R2"` no longer expands to `001_R2` and `002_R2`.
+  A kcat list with one value per matched row takes its values by position even when another
+  ID matches nothing, as in MATLAB: `["R2", "nonsense"]` with `[1.0, 2.0]` sets the two `R2` rows to 1.0 and 2.0.
+* `flexibilize_enz_concs` with a levelled complex scores a unit by its probes when the solver
+  returns no usable reduced costs (the scipy solver), where a NaN score could pick the first unit.
+* `get_conc_control_coeffs` keeps a reduced cost when raising the bound by 0.01 percent raises
+  the optimum by at least half of the rise the reduced cost predicts, and by more than 1e-12
+  times max(1, |optimum|). A small coefficient on a small bound is kept; the test used to need a
+  rise above 1e-10, which dropped any coefficient with a reduced cost times bound below 1e-6.
+* `fit_sigma(method="bisect")` gives the first grid sigma, `1 / n_sigma_steps`, for a model that
+  cannot grow even at sigma 1, as the grid does; it used to give about 1e-6.
+* `truncate_values` rounds a value exactly halfway between two candidates away from
+  zero, as MATLAB's `round(x, n)` does: 100000.5 gives 100001 and 1234565 gives
+  1234570, where they gave 100000 and 1234560.
+* `save_kcat_source` writes the free text of the `detail` column as
+  `saveKcatSource` does: a tab or line break becomes a space, so a row can no longer
+  be split, and an empty part leaves no leading, trailing or repeated `; `.
+* `fill_kcats_from_isozymes` and `get_kcat_across_isozymes` default to `apply=False`,
+  as `getKcatAcrossIsozymes` only changes `ec.kcat`; pass `apply=True` to rewrite the
+  stoichiometry. Only a kcat of 0 is filled: a NaN kcat is left as it is, and a NaN
+  sibling does not enter the mean. `assign_kcats` passes `apply=True`, as before.
+* `fit_sigma` scans the 0.01 grid of GECKO's `sigmaFitter` by default (`method="bisect"`
+  stays available and still gives a continuous sigma). A target that growth cannot
+  reach gives the first sigma of the growth plateau, where it used to give 1.0. The
+  bounds of the objective reaction are set to `[0, 1000]` before the scan and stay so on
+  the fitted model, and a solve that is not optimal counts as growth 0.
+* `sensitivity_tuning` follows `sensitivityTuning.m`. The result lists the tuned reactions
+  in model order, and reports the kcat and source of the input model, so a reaction an
+  earlier run had tuned no longer reports NaN. `ec.notes` prints the pre-tune kcat as
+  `num2str` does (`preTuneKcat=10`). A `prot_to_ignore` id is matched without regard to
+  case and raises `ValueError` when no usage reaction has it. A tuned kcat has no upper
+  limit. A gecko-light model is tuned, in place of raising `NotImplementedError`.
+* `get_conc_control_coeffs` also takes the one-solve path with the scipy solver. optlang does not expose
+  variable duals there, so the reduced cost is read from the HiGHS bound marginal of the last `linprog`
+  result, after checking that the stored result is the solution `model.optimize` returned. On the tutorial
+  yeast ecModel this takes seconds where the finite-difference loop took several minutes. A scipy result
+  without bound marginals still falls back to the finite-difference loop.
+* `get_conc_control_coeffs` reads the reduced cost of each usage reaction's forward
+  variable and keeps it where a 0.01 percent rise of the bound raises the optimum,
+  in place of the `prot_<id>` shadow price. An enzyme whose own bound does not
+  bind now gets 0 while the protein pool binds, where every used enzyme used to
+  get the pool's price, and a degenerate vertex no longer gives a spurious
+  coefficient. Every nonzero reduced cost is checked, however small, and a usage
+  reaction held at its lower bound gets 0 without a check. Where the solver exposes
+  no variable duals, or returns NaN, the finite-difference loop is used.
 * New: `save_ec_model` records the ecModel format (`gecko_format: "4"`) and
   the release that wrote the file (`generator: geckopy <version>`) in the
   `metaData`, the keys GECKO's `saveEcModel` writes, in place of

docs/gecko_light_status.md

@@ -38,7 +38,7 @@ smaller models (yeast, *E. coli*), full is usually fine.
 
 ```python
 ec_model = make_ec_model(model, adapter, gecko_light=True)
-set_kcat_for_reactions(ec_model, ["R2"], 5.0)   # broadcasts to all isozyme rows
+set_kcat_for_reactions(ec_model, ["001_R2", "002_R2"], 5.0)   # one row per isozyme
 apply_kcat_constraints(ec_model)                 # picks cheapest isozyme per cobra rxn
 set_prot_pool_size(ec_model)
 ec_model.optimize()
@@ -67,7 +67,7 @@ isn't checked out next to geckopy.
 | `src/geckopy/ec_model/make_ec_model.py` | dispatches stage 5 (isozyme split) and stages 9/11 (per-enzyme prot mets + usage rxns) off; calls `allocate_ec_and_coupling_light` instead of the full stages 6/8 |
 | `src/geckopy/ec_model/pipeline/populate_ec.py` | `allocate_ec_and_coupling_light` builds `ec.rxns` (one row per isozyme with the `###_` prefix) and `ec.rxn_enz_mat` in one pass; `split_light_rxn_id` is the inverse |
 | `src/geckopy/ec_model/pipeline/apply_kcat.py` | `_apply_kcat_constraints_light` groups ec rows by cobra reaction, picks the lowest-cost isozyme, writes `-MW_sum / (kcat * 3600)` as the `prot_pool` coefficient; same clear-then-write idempotency contract as full |
-| `src/geckopy/ec_model/pipeline/set_kcat.py` | recognises the `###_` prefix as the light analogue of `_EXP_<N>`; bare base names broadcast to every isozyme row |
+| `src/geckopy/ec_model/pipeline/set_kcat.py` | the `###_` prefix is part of the id in a light model: `001_R2` matches that one row and a bare cobra id such as `R2` matches nothing |
 | `src/geckopy/ec_model/enzyme.py` | `Enzyme.reactions` strips the prefix and dedupes; `Enzyme.mw` setter passes ec.rxns ids (not cobra ids) to `apply_kcat_constraints`; `_repr_html_` catches the light NotImplementedError |
 | `src/geckopy/gather_kcats/fuzzy_kcat_matching.py` | strips the `###_` prefix for cobra lookups |
 | `src/geckopy/gather_kcats/get_standard_kcat.py` | `_ec_rxn_to_cobra_id` handles both forms |

docs/migrating_from_gecko_matlab.md

@@ -228,7 +228,7 @@ Each row maps a MATLAB function to its geckopy equivalent. The
 | `getComplexData`             | `get_complex_data`            | —                                                                                                                                                                              |
 | `getKcatAcrossIsozymes`      | `fill_kcats_from_isozymes`    | **`[Py]`** Renamed; `get_kcat_across_isozymes` kept as deprecated alias. Aggregation strategy (`max` / `median` / `mean`) is adapter-configurable via `kcat_aggregate_isozymes`. |
 | `getReactionsFromEnzyme`     | `get_reactions_from_enzyme`   | **`[Py]`** Returns one `DataFrame` instead of five parallel cells.                                                                                                            |
-| `setKcatForReactions`        | `set_kcat_for_reactions`      | **`[Py]`** Recognises both `_EXP_<N>` (full) and `###_` (light) suffixes; understands base names that broadcast across isozymes.                                              |
+| `setKcatForReactions`        | `set_kcat_for_reactions`      | Matches IDs without regard to letter case, ignores an unknown ID and does not apply the constraints unless `apply=True`, as `setKcatForReactions`; a light-model row is addressed by its `###_` ID.                    |
 | `setProtPoolSize`            | `set_prot_pool_size`          | **`[3→4]`** Sets the **upper** bound of `prot_pool_exchange` (forward direction).                                                                                              |
 
 ### 5.2 Gather kcats (`src/geckomat/gather_kcats`)
@@ -300,7 +300,7 @@ Each row maps a MATLAB function to its geckopy equivalent. The
 | `constrainFluxData`   | `apply_flux_data_constraints`        | **`[Py]`** Renamed; `constrain_flux_data` kept as deprecated alias. Takes pre-loaded `FluxData`. 0-indexed `condition` (default `0`, MATLAB `1`).                                                                                                                                                              |
 | `fillEnzConcs`        | `fill_enz_concs`                     | —                                                                                                                                                                                                                                                                                                              |
 | `flexibilizeEnzConcs` | `flexibilize_enz_concs`              | **`[Py]`** Returns a `FlexEnzResult` dataclass. Raises `NotImplementedError` on gecko-light.                                                                                                                                                                                                                  |
-| `getConcControlCoeffs` | `get_conc_control_coeffs`            | **`[3→4]`** Probe direction reads the forward `usage_prot_*` upper bound. **`[Py]`** Default path uses the prot-met shadow price (one LP solve, *n* times fewer LPs than the finite-difference probe); falls back to finite-difference when the solver doesn't expose LP duals. Gates on `sol.status == "optimal"`. |
+| `getConcControlCoeffs` | `get_conc_control_coeffs`            | **`[3→4]`** Probe direction reads the forward `usage_prot_*` upper bound. **`[Py]`** Default path reads the reduced cost of each usage reaction's forward variable (one LP solve plus one check solve per nonzero value, against one solve per enzyme for the finite-difference probe), so an enzyme whose own bound does not bind gets 0 also while the protein pool binds; falls back to finite-difference when the solver exposes no variable duals or returns NaN. Gates on `sol.status == "optimal"`. |
 | `updateProtPool`      | *(not ported)*                       | Obsolete since GECKO 3.2.0. Use `set_prot_pool_size` instead.                                                                                                                                                                                                                                                  |
 | —                     | `relax_proteomics_greedy` *(new)*    | **`[Py]`** Alternative to `flexibilize_enz_concs`: relaxes a fixed top-*k* of measured concentrations in a greedy loop until a target growth is reached. Returns a `GreedyRelaxResult`.                                                                                                                       |
 

docs/relaxation_methods.md

@@ -14,8 +14,8 @@ the per-complex options for all of them are described at the end.
 | aspect | `flexibilize_enz_concs` | `relax_proteomics_greedy` |
 |---|---|---|
 | Origin | GECKO MATLAB `flexibilizeEnzConcs.m` (faithful port) | legacy geckopy, Carrasco et al. 2023 |
-| Which enzyme to relax | largest **control coefficient** (finite-difference: how much each measured enzyme's bound limits growth, from `get_conc_control_coeffs`) | largest **\|shadow price\|** (LP dual of the usage upper bound) |
-| Solves per selection | re-solves the LP **once per candidate protein** (expensive) | a **single** `slim_optimize` exposes all duals (cheap) |
+| Which enzyme to relax | largest **control coefficient** (how much each measured enzyme's bound limits growth, from `get_conc_control_coeffs`: the reduced cost of its usage bound, or a finite difference where the solver exposes no duals) | largest **\|shadow price\|** (LP dual of the usage upper bound) |
+| Solves per selection | one LP solve plus one check solve per enzyme with a nonzero coefficient (one solve per candidate protein where the solver exposes no duals) | a **single** `slim_optimize` exposes all duals (cheap) |
 | Relaxation step | **gradual**: usage `ub ← conc·(1 + fold_change·k)`; the same enzyme can be picked repeatedly and loosened a bit more each time | **all-or-nothing**: usage `ub ← default_upper_bound` (≈ unconstrained) in one shot |
 | Tighten-back | yes — a refinement pass minimises the protein pool at `exp_growth` and **drops** enzymes whose true usage is below their original conc (un-relaxes over-relaxed ones) | none — each picked enzyme stays fully unconstrained |
 | Result stringency | **tight** (minimal, proteomics-faithful relaxation) | **loose** (tends to over-relax) |

src/geckopy/ec_model/make_ec_model.py

@@ -23,6 +23,7 @@
 from __future__ import annotations
 
 import logging
+import re
 from typing import TYPE_CHECKING, Optional
 
 import cobra
@@ -52,6 +53,50 @@
 
 logger = logging.getLogger(__name__)
 
+_EXP_SUFFIX = re.compile(r"_EXP_\d+$")
+
+
+def _check_inputs(model: cobra.Model, adapter: "ModelAdapter") -> None:
+    """Refuse inputs that GECKO's ``makeEcModel`` refuses."""
+    enzyme_comp = adapter.params.enzyme_comp
+    if enzyme_comp not in set(model.compartments.values()):
+        raise ValueError(
+            f"Compartment {enzyme_comp} (specified in params.enzyme_comp) "
+            "cannot be found in the compartment names of the model."
+        )
+    if any(met.id.startswith("prot_") for met in model.metabolites):
+        raise ValueError(
+            "The identifiers in model.metabolites are not allowed to start "
+            "with 'prot_'."
+        )
+    if any(
+        rxn.id.startswith(("usage_prot_", "prot_pool"))
+        or rxn.id.endswith("_REV")
+        or _EXP_SUFFIX.search(rxn.id)
+        for rxn in model.reactions
+    ):
+        raise ValueError(
+            "The identifiers in model.reactions are not allowed to start "
+            "with 'usage_prot' or 'prot_pool', or end with '_REV' or "
+            "'_EXP_[digit]'."
+        )
+    if not model.genes:
+        raise ValueError(
+            "None of the proteins in uniprot.tsv or kegg.tsv match the genes "
+            "in the model: the model has no genes."
+        )
+
+
+def _is_exchange(rxn: cobra.Reaction) -> bool:
+    """Whether a reaction only consumes or only produces metabolites.
+
+    RAVEN's ``getExchangeRxns`` counts every reaction without a positive or
+    without a negative stoichiometric coefficient, so ``A + B <=>`` is an
+    exchange reaction as well as ``A <=>``.
+    """
+    coefficients = list(rxn.metabolites.values())
+    return all(c <= 0 for c in coefficients) or all(c >= 0 for c in coefficients)
+
 
 def make_ec_model(
     model: cobra.Model,
@@ -131,9 +176,15 @@ def make_ec_model(
         If ``uniprot_db`` is None and no ``uniprot.tsv`` is found
         under the adapter's data folder.
     ValueError
-        Propagated from individual pipeline stages on
-        inconsistencies (e.g. genes referenced in GPRs that don't
-        exist in the model).
+        As ``makeEcModel`` does, if the compartment named by
+        ``adapter.params.enzyme_comp`` is not one of the model's
+        compartment names (a compartment ID that differs from its name is
+        refused), if a metabolite ID starts with ``prot_``, if a reaction ID
+        starts with ``usage_prot_`` or ``prot_pool`` or ends in ``_REV`` or
+        ``_EXP_<digits>``, or if the model has no genes or none of them is
+        found in the databases. Also propagated from individual pipeline
+        stages on inconsistencies (e.g. genes referenced in GPRs that
+        don't exist in the model).
 
     Genes referenced in the model's GPRs but not found in UniProt are
     left enzyme-unconstrained: they are reported as a single warning
@@ -146,6 +197,8 @@ def make_ec_model(
             "populated ec substructure. Run it only on a conventional GEM."
         )
 
+    _check_inputs(model, adapter)
+
     # Work on a copy: the pipeline stages below mutate the model in place
     # (splitting reactions, adding pseudometabolites, ...), but the caller's
     # input GEM must be left untouched.
@@ -156,11 +209,17 @@ def make_ec_model(
     # them per-row in ec instead.
     remove_pseudoreaction_gprs(model, adapter)
     invert_backwards_only_reactions(model)
-    # Exchange reactions are left reversible: MATLAB's makeEcModel.m
-    # excludes them here too, passing convertToIrrev only nonExchRxns
-    # (model.rxns minus getExchangeRxns(model)).
-    non_exchange_rxns = [rxn.id for rxn in model.reactions if not rxn.boundary]
-    convert_to_irreversible(model, rxns=non_exchange_rxns)
+    # Only reactions that carry flux in both directions (lb < 0 < ub) and
+    # are not exchange reactions are split, as MATLAB's makeEcModel.m does:
+    # it passes convertToIrrev the reversible reactions without
+    # getExchangeRxns(model). A reaction with both bounds negative stays as
+    # it is.
+    reversible_rxns = [
+        rxn.id
+        for rxn in model.reactions
+        if rxn.lower_bound < 0 < rxn.upper_bound and not _is_exchange(rxn)
+    ]
+    convert_to_irreversible(model, rxns=reversible_rxns)
     if not gecko_light:
         expand_model(model)
         # Sort reactions, so that reversible and isozymic reactions are kept

src/geckopy/ec_model/pipeline/apply_custom_kcats.py

@@ -41,8 +41,10 @@ def apply_custom_kcats(
     Mode A (rxns only, proteins empty)
         Apply the kcat to every entry of ``ec.rxns`` whose ID, after
         stripping any ``_EXP_<n>`` suffix, equals one of the listed
-        reaction IDs (comma-separated). Reverse direction is honored:
-        ``"R2"`` matches ``R2_EXP_*`` but not ``R2_REV_*``.
+        reaction IDs (comma-separated). A list in double quotes, as a
+        spreadsheet exports it, is read without the quotes. Reverse
+        direction is honored: ``"R2"`` matches ``R2_EXP_*`` but not
+        ``R2_REV_*``.
 
     Mode B (proteins only, rxns empty)
         Candidate reactions are those catalyzed by any listed protein.
@@ -266,13 +268,24 @@ def _parse_custom_kcats_tsv(path: Path) -> list[dict]:
                 "index": line_no - 1,  # 1-based row index excluding header
                 "proteins": parts[0].strip(),
                 "kcat": kcat,
-                "rxns": parts[4].strip(),
+                "rxns": _unquote(parts[4].strip()),
                 "notes": parts[5].strip(),
                 "stoicho": parts[6].strip(),
             })
     return rows
 
 
+def _unquote(field: str) -> str:
+    """Strip one pair of double quotes around a field, as MATLAB's ``%q`` does.
+
+    A spreadsheet exports a reaction list such as ``R3,R5`` as ``"R3,R5"``;
+    a doubled quote inside the field stands for one quote.
+    """
+    if len(field) >= 2 and field.startswith('"') and field.endswith('"'):
+        return field[1:-1].replace('""', '"')
+    return field
+
+
 def _resolve_rxn_ids(
     rxns_field: str, no_suffix: list[str], ec_rxn_ids: list[str],
 ) -> list[int]: