diff --git a/CHANGELOG.md b/CHANGELOG.md index cadf5dd..d2a9f5b 100644 --- a/CHANGELOG.md +++ b/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_`, 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_` 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 `) in the `metaData`, the keys GECKO's `saveEcModel` writes, in place of diff --git a/docs/gecko_light_status.md b/docs/gecko_light_status.md index 2536b01..1a7a2e3 100644 --- a/docs/gecko_light_status.md +++ b/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_`; 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 | diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index b4cefa4..541777f 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/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_` (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`. | diff --git a/docs/relaxation_methods.md b/docs/relaxation_methods.md index 3b823d5..b49781e 100644 --- a/docs/relaxation_methods.md +++ b/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) | diff --git a/src/geckopy/ec_model/make_ec_model.py b/src/geckopy/ec_model/make_ec_model.py index 5d68e54..5734201 100644 --- a/src/geckopy/ec_model/make_ec_model.py +++ b/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_``, 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 diff --git a/src/geckopy/ec_model/pipeline/apply_custom_kcats.py b/src/geckopy/ec_model/pipeline/apply_custom_kcats.py index 272dcbc..736472a 100644 --- a/src/geckopy/ec_model/pipeline/apply_custom_kcats.py +++ b/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_`` 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]: diff --git a/src/geckopy/ec_model/pipeline/apply_kcat.py b/src/geckopy/ec_model/pipeline/apply_kcat.py index f0378c5..12b0557 100644 --- a/src/geckopy/ec_model/pipeline/apply_kcat.py +++ b/src/geckopy/ec_model/pipeline/apply_kcat.py @@ -59,8 +59,10 @@ def apply_kcat_constraints( An EcModel produced by ``make_ec_model``, with at least some kcat values assigned to ``model.ec.kcat``. Mutated in place. update_rxns - Reaction IDs (from ``model.ec.rxns``) to update. If None, - updates all entries. + Reaction IDs (from ``model.ec.rxns``) to update. ``None`` or an + empty sequence updates all entries, as ``applyKcatConstraints`` + does with an empty ``updateRxns``; a caller that builds the list + from a filter must skip the call when the filter finds nothing. Raises ------ @@ -80,7 +82,10 @@ def apply_kcat_constraints( ec_rxn_ids = model.ec.rxns # Resolve update_rxns to a numpy array of ec-row indices. - if update_rxns is None: + if update_rxns is not None: + update_rxns = list(update_rxns) + + if not update_rxns: selected_idx = np.arange(len(ec_rxn_ids)) else: id_to_idx = {rxn_id: i for i, rxn_id in enumerate(ec_rxn_ids)} @@ -194,13 +199,17 @@ def _apply_kcat_constraints_light( model A gecko-light EcModel with ``ec.gecko_light is True``. update_rxns - ec.rxns IDs (with the ``###_`` prefix) to update. ``None`` - means update every cobra reaction backed by at least one - ec.rxns row. + ec.rxns IDs (with the ``###_`` prefix) to update. ``None`` or an + empty sequence means update every cobra reaction backed by at + least one ec.rxns row. A reaction is updated as a whole: naming one of its + isozyme rows recomputes the cheapest cost over all its rows. """ ec_rxn_ids = model.ec.rxns - if update_rxns is None: + if update_rxns is not None: + update_rxns = list(update_rxns) + + if not update_rxns: selected_idx = np.arange(len(ec_rxn_ids)) else: id_to_idx = {rxn_id: i for i, rxn_id in enumerate(ec_rxn_ids)} @@ -217,12 +226,17 @@ def _apply_kcat_constraints_light( if selected_idx.size == 0: return - # Group the selected ec rows by the cobra reaction they belong to; - # the cobra reaction is what carries the LP constraint. + # The cobra reaction carries the LP constraint, so selecting one isozyme + # row updates the reaction as a whole: every ec row of that reaction + # competes for the cheapest cost, selected or not. + affected = { + split_light_rxn_id(ec_rxn_ids[idx])[1] for idx in selected_idx + } rows_by_cobra: dict[str, list[int]] = defaultdict(list) - for idx in selected_idx: - _, cobra_id = split_light_rxn_id(ec_rxn_ids[idx]) - rows_by_cobra[cobra_id].append(int(idx)) + for idx, ec_rxn_id in enumerate(ec_rxn_ids): + _, cobra_id = split_light_rxn_id(ec_rxn_id) + if cobra_id in affected: + rows_by_cobra[cobra_id].append(idx) pool_met = model.metabolites.get_by_id(POOL_ID) rxn_enz_mat_csr = model.ec.rxn_enz_mat.tocsr() diff --git a/src/geckopy/ec_model/pipeline/fill_kcats.py b/src/geckopy/ec_model/pipeline/fill_kcats.py index c16df8a..8af9f1a 100644 --- a/src/geckopy/ec_model/pipeline/fill_kcats.py +++ b/src/geckopy/ec_model/pipeline/fill_kcats.py @@ -32,7 +32,7 @@ def fill_kcats_from_isozymes( model: "EcModel", *, - apply: bool = True, + apply: bool = False, aggregate: Optional[str] = None, ) -> None: """Fill in missing kcats by combining known kcats across sibling isozymes. @@ -45,7 +45,10 @@ def fill_kcats_from_isozymes( For each reaction with no kcat assigned (``ec.kcat == 0``), find its siblings whose kcat is set (``> 0``). If any exist, take their mean and assign it. Set ``ec.source`` to ``"isozymes"`` for the - filled entries. + filled entries. Only a kcat of 0 is filled, as in ``getKcatAcrossIsozymes``: + a NaN kcat is not a missing kcat and is left as it is. A NaN sibling + does not contribute to the mean here; in ``getKcatAcrossIsozymes`` it + does, so there the mean is NaN and the row is not filled. Reactions whose siblings all also lack a kcat stay at 0. Reactions with a single isozyme (no siblings beyond themselves) also stay at 0. @@ -58,8 +61,9 @@ def fill_kcats_from_isozymes( An EcModel produced by ``make_ec_model``. Mutated in place. Must not be a gecko-light model. apply - If True (default), call ``apply_kcat_constraints`` after - updating ec.kcat so the new values reflect in the S matrix. + If True, call ``apply_kcat_constraints`` after updating ec.kcat + so the new values reflect in the S matrix. Default False, as + in ``getKcatAcrossIsozymes``, which only changes ``ec.kcat``. aggregate How to combine sibling kcats: ``"mean"`` (matches MATLAB GECKO), ``"median"`` (more robust for log-distributed rate constants), or @@ -88,6 +92,7 @@ def fill_kcats_from_isozymes( kcat = model.ec.kcat known_mask = kcat > 0 + missing_mask = kcat == 0 if kcat.size == 0 or not known_mask.any(): logger.warning( "fill_kcats_from_isozymes: ec.kcat has no known values to " @@ -113,8 +118,8 @@ def fill_kcats_from_isozymes( # For each missing entry, look up its base group and average if available. filled_indices: list[int] = [] - for i, (bid, known) in enumerate(zip(base_ids, known_mask)): - if known: + for i, (bid, missing) in enumerate(zip(base_ids, missing_mask)): + if not missing: continue siblings = known_by_base.get(bid) if not siblings: @@ -140,7 +145,7 @@ def fill_kcats_from_isozymes( apply_kcat_constraints(model, update_rxns=filled_rxn_ids) -def get_kcat_across_isozymes(model: "EcModel", *, apply: bool = True) -> None: +def get_kcat_across_isozymes(model: "EcModel", *, apply: bool = False) -> None: """Deprecated alias for :func:`fill_kcats_from_isozymes`. Kept for backward compatibility with the original MATLAB name. diff --git a/src/geckopy/ec_model/pipeline/protein_pool.py b/src/geckopy/ec_model/pipeline/protein_pool.py index e7bc8cc..14047a4 100644 --- a/src/geckopy/ec_model/pipeline/protein_pool.py +++ b/src/geckopy/ec_model/pipeline/protein_pool.py @@ -39,9 +39,9 @@ def _resolve_enzyme_compartment_id(model: "EcModel", requested: str) -> str: """Resolve ``adapter.params.enzyme_comp`` to a compartment ID. Matches ``requested`` first against compartment names (the typical - case, e.g. "cytoplasm"), then falls back to compartment IDs (e.g. - "c" when the model happens to use short IDs as names, as in - ecTestGEM). + case, e.g. "cytoplasm"), then falls back to compartment IDs. This + resolver accepts the ID fallback; ``make_ec_model`` itself requires + ``enzyme_comp`` to be a compartment name. Raises ValueError if neither matches. """ diff --git a/src/geckopy/ec_model/pipeline/set_kcat.py b/src/geckopy/ec_model/pipeline/set_kcat.py index 62be5c6..d20db97 100644 --- a/src/geckopy/ec_model/pipeline/set_kcat.py +++ b/src/geckopy/ec_model/pipeline/set_kcat.py @@ -9,14 +9,12 @@ from .apply_kcat import apply_kcat_constraints -from .populate_ec import split_light_rxn_id if TYPE_CHECKING: from ..ec_model import EcModel _EXP_SUFFIX_PATTERN = re.compile(r"_EXP_\d+$") -_LIGHT_PREFIX_PATTERN = re.compile(r"^\d{3}_") _SOURCE_TAG = "setKcatForReactions" @@ -25,26 +23,30 @@ def set_kcat_for_reactions( rxn_ids: Sequence[str], kcat: Union[float, Sequence[float]], *, - apply: bool = True, + apply: bool = False, ) -> list[str]: """Set kcat values for one or more reactions in ec.kcat. - Each ID in ``rxn_ids`` is interpreted as follows: + As GECKO's ``setKcatForReactions``, each ID in ``rxn_ids`` is matched + against ``ec.rxns`` without regard to letter case: - - If the ID ends in ``_EXP_`` (full layout) or starts with - ``###_`` (gecko-light), only that exact reaction is matched. - - Otherwise, the ID is treated as a base name and matches every - ec.rxns entry whose ID, after stripping the layout-specific - isozyme marker, equals the base name. So in a full model, - ``"R2"`` matches ``R2``, ``R2_EXP_1``, ``R2_EXP_2``; in a light - model, ``"R2"`` matches ``001_R2``, ``002_R2``. + - An ID that ends in ``_EXP_`` matches the ec.rxns entry with + exactly that ID. + - Any other ID matches every ec.rxns entry whose ID, after removing a + trailing ``_EXP_``, equals it. So in a full model ``"R2"`` + matches ``R2``, ``R2_EXP_1`` and ``R2_EXP_2``. In a gecko-light + model the ec.rxns IDs carry a ``###_`` prefix (``001_R2``), which + is part of the ID: ``"R2"`` matches nothing there and + ``"001_R2"`` matches that one row. + - An ID that matches nothing is ignored. The ``kcat`` argument follows numpy-style broadcasting: a single - float applies to every matched reaction; a sequence of length - ``len(rxn_ids)`` gives one value per input ID (broadcast across - that ID's own matches); a sequence of length equal to the total - number of matched reactions gives one value per resolved match, in - ``ec.rxns`` order. + float applies to every matched reaction; a sequence of length equal + to the total number of matched reactions gives one value per + resolved match, in ``ec.rxns`` order; a sequence of length + ``len(rxn_ids)`` otherwise gives one value per input ID (broadcast + across that ID's own matches). When both lengths fit, the per-match + form applies, so an ID that matches nothing never shifts the values. That last form lets one un-suffixed ``rxn_id`` that expands to multiple isozymes take a list of per-isozyme values, matching @@ -69,10 +71,11 @@ def set_kcat_for_reactions( kcat New kcat value(s) in 1/s. Scalar or sequence (see above). apply - If True (default), call ``apply_kcat_constraints`` after - updating ec.kcat / ec.source so the change is reflected in the - S matrix immediately. If False, the caller must invoke - ``apply_kcat_constraints`` themselves. + If True, call ``apply_kcat_constraints`` for the changed + reactions after updating ec.kcat / ec.source so the change is + reflected in the S matrix immediately. If False (default, as in + MATLAB, which only changes ec.kcat and ec.source), the caller + invokes ``apply_kcat_constraints`` afterwards. Returns ------- @@ -83,37 +86,28 @@ def set_kcat_for_reactions( Raises ------ ValueError - If any rxn_id matches zero reactions, or if kcat has a length - other than 1, ``len(rxn_ids)``, or the total number of resolved - matches. + If kcat has a length other than 1, ``len(rxn_ids)``, or the + total number of resolved matches. """ rxn_ids = list(rxn_ids) if not rxn_ids: return [] ec_rxns = model.ec.rxns - if model.ec.gecko_light: - # Light: strip the leading ``###_`` counter so the cobra reaction - # id is what the base-name lookup sees. - base_ec_rxns = [split_light_rxn_id(r)[1] for r in ec_rxns] - explicit_pattern = _LIGHT_PREFIX_PATTERN - else: - base_ec_rxns = [_EXP_SUFFIX_PATTERN.sub("", r) for r in ec_rxns] - explicit_pattern = _EXP_SUFFIX_PATTERN + ec_rxns_lower = [r.lower() for r in ec_rxns] + base_ec_rxns_lower = [_EXP_SUFFIX_PATTERN.sub("", r).lower() for r in ec_rxns] # Resolve every input ID to a list of ec-row indices. matches_per_input: list[list[int]] = [] for rxn_id in rxn_ids: - if explicit_pattern.search(rxn_id): - indices = [i for i, r in enumerate(ec_rxns) if r == rxn_id] + key = rxn_id.lower() + if _EXP_SUFFIX_PATTERN.search(rxn_id): + indices = [i for i, r in enumerate(ec_rxns_lower) if r == key] else: + base_key = _EXP_SUFFIX_PATTERN.sub("", key) indices = [ - i for i, r in enumerate(base_ec_rxns) if r == rxn_id + i for i, r in enumerate(base_ec_rxns_lower) if r == base_key ] - if not indices: - raise ValueError( - f"rxn_id '{rxn_id}' matched no entries in ec.rxns." - ) matches_per_input.append(indices) total_matches = sum(len(m) for m in matches_per_input) @@ -123,24 +117,24 @@ def set_kcat_for_reactions( kcat_per_index = [float(kcat)] * total_matches else: kcat_seq = list(kcat) - # kcat_seq has length len(rxn_ids) (one value per input ID, - # broadcast across that ID's own matches), or length - # total_matches (one value per resolved match, in ec.rxns - # order -- MATLAB setKcatForReactions's positional semantics; - # see the docstring's fragility note). When every input matches - # exactly one row the two are the same list, so there's no - # ambiguity between them. - if len(kcat_seq) == len(rxn_ids): + # kcat_seq has length total_matches (one value per resolved match, + # in input order -- MATLAB setKcatForReactions's positional + # semantics), or length len(rxn_ids) (one value per input ID, + # broadcast across that ID's own matches). The positional form is + # tried first, so an input ID without a match never shifts the + # values the way the per-input form would; when every input + # matches exactly one row the two forms are the same list. + if len(kcat_seq) == total_matches: + kcat_per_index = [float(v) for v in kcat_seq] + elif len(kcat_seq) == len(rxn_ids): kcat_per_index = [] for value, indices in zip(kcat_seq, matches_per_input): kcat_per_index.extend([float(value)] * len(indices)) - elif len(kcat_seq) == total_matches: - kcat_per_index = [float(v) for v in kcat_seq] else: raise ValueError( f"kcat has length {len(kcat_seq)}; expected a scalar, " - f"length {len(rxn_ids)} (one per input rxn_id), or length " - f"{total_matches} (one per resolved match)." + f"length {total_matches} (one per resolved match) or " + f"length {len(rxn_ids)} (one per input rxn_id)." ) # Apply to ec.kcat and ec.source. @@ -154,7 +148,7 @@ def set_kcat_for_reactions( updated_ids = [ec_rxns[i] for i in flat_indices] - if apply: + if apply and updated_ids: apply_kcat_constraints(model, update_rxns=updated_ids) return updated_ids diff --git a/src/geckopy/gather_kcats/kcat_assignment.py b/src/geckopy/gather_kcats/kcat_assignment.py index 4b4548e..99aa80a 100644 --- a/src/geckopy/gather_kcats/kcat_assignment.py +++ b/src/geckopy/gather_kcats/kcat_assignment.py @@ -252,7 +252,7 @@ def load(name: str) -> Optional[pd.DataFrame]: stage[i], candidates[i], spread[i] = "custom", 1, 0.0 if config.isozymes and not ec.gecko_light: before = ec.kcat == 0 - fill_kcats_from_isozymes(model) + fill_kcats_from_isozymes(model, apply=True) for i in np.flatnonzero(before & (ec.kcat > 0)): stage[i] = "isozymes" if config.standard: diff --git a/src/geckopy/gather_kcats/kcat_sources.py b/src/geckopy/gather_kcats/kcat_sources.py index d9c987e..41ef74c 100644 --- a/src/geckopy/gather_kcats/kcat_sources.py +++ b/src/geckopy/gather_kcats/kcat_sources.py @@ -207,6 +207,23 @@ def _clean(value) -> str: return re.sub(r"[\t\r\n]", " ", str(value)) +def _detail_text(wildcard_level, origin, extra) -> str: + """The detail column: ``wc=``, ``origin=`` and the free text, joined by ``; ``. + + Tabs and line breaks of the free text become spaces, so a row can never be split. Empty + parts leave no ``; `` behind: a leading or trailing ``; `` is dropped and a repeated + ``; `` is written once, as ``saveKcatSource`` does. + """ + parts = [ + f"wc={int(wildcard_level)}" if wildcard_level is not None and pd.notna(wildcard_level) else "", + f"origin={int(origin)}" if origin is not None and pd.notna(origin) else "", + _clean(extra).strip(), + ] + text = "; ".join(parts) + text = re.sub(r"^(; )+|(; )+$", "", text) + return re.sub(r"(; ){2,}", "; ", text) + + def save_kcat_source( model: "EcModel", kcat_list: pd.DataFrame, @@ -284,16 +301,9 @@ def col(column): if not enzymes[i]: no_enzymes += 1 continue - detail = [] - if wc is not None and pd.notna(wc): - detail.append(f"wc={int(wc)}") - if origin is not None and pd.notna(origin): - detail.append(f"origin={int(origin)}") - if extra is not None and not (isinstance(extra, float) and math.isnan(extra)) and str(extra).strip(): - detail.append(str(extra).strip()) rows.append([ _clean(base[i]), _clean(enzymes[i]), _clean(subs), format_number(kcat), - format_number(sd), _clean(src), "; ".join(detail), _clean(ec_code), "", + format_number(sd), _clean(src), _detail_text(wc, origin, extra), _clean(ec_code), "", ]) if no_enzymes: logger.info("%d kcats of reactions without enzymes are not written to %s.", no_enzymes, path) diff --git a/src/geckopy/get_enzyme_data/find_ec_in_db.py b/src/geckopy/get_enzyme_data/find_ec_in_db.py index f97a0ab..595cbbf 100644 --- a/src/geckopy/get_enzyme_data/find_ec_in_db.py +++ b/src/geckopy/get_enzyme_data/find_ec_in_db.py @@ -6,12 +6,26 @@ from __future__ import annotations import logging +import re import numpy as np logger = logging.getLogger(__name__) +_EC_SEPARATORS = re.compile(r"[;\s]+") + + +def _split_ec_cell(cell: str) -> list[str]: + """Split a database EC cell into its EC tokens. + + UniProt separates codes with ``"; "`` and KEGG with a space, so a cell is + split on semicolons and whitespace alike and empty pieces are dropped, + as GECKO MATLAB's ``getECstring`` does. + """ + return [token for token in _EC_SEPARATORS.split(cell) if token] + + def _split_levels(token: str) -> list[str]: return token.split(".") @@ -122,8 +136,9 @@ def find_ec_in_db( cross-gene UNION / INTERSECTION ordering. db_eccodes EC string per protein row in the database. Empty string means - "no EC for this protein". Multi-EC strings (``;``-joined) are - permitted. + "no EC for this protein". A cell may hold several codes, separated + by semicolons, spaces or both (``"1.1.1.1; 2.2.2.2"`` from UniProt, + ``"1.1.1.4 1.1.1.-"`` from KEGG). db_mw Molecular weight per protein row, in Da. NaN is tolerated. gene_to_protein_indices @@ -186,7 +201,7 @@ def find_ec_in_db( min_idx = int(np.argmin(finite_mws)) chosen = db_eccodes[matched[min_idx]] - per_gene_tokens.append(chosen.split(";")) + per_gene_tokens.append(_split_ec_cell(chosen)) per_gene_tokens = [t for t in per_gene_tokens if t] diff --git a/src/geckopy/kcat_tuning/evotune/tuning.py b/src/geckopy/kcat_tuning/evotune/tuning.py index e646d19..6c448f2 100644 --- a/src/geckopy/kcat_tuning/evotune/tuning.py +++ b/src/geckopy/kcat_tuning/evotune/tuning.py @@ -330,7 +330,8 @@ def screen_kcat_leverage( # scores its own copy instead, so this is a no-op there. Either # way, `model` must come out holding the prior it went in with. model.ec.kcat[tunable_idx] = kcat0 - apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) + if ec_rxn_ids_tunable: + apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) base_rmse = rmses[0] rows = [] @@ -578,7 +579,8 @@ def _score_batch(kcat_matrix, pool): break model.ec.kcat[tunable_idx] = best_vec - apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) + if ec_rxn_ids_tunable: + apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) return EvotuneResult( rxns=ec_rxn_ids_tunable, old_kcat=kcat0.copy(), new_kcat=best_vec, @@ -832,7 +834,8 @@ def _score_kcat_vector( across every candidate it's asked to score. """ model.ec.kcat[tunable_idx] = kcat_vec - apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) + if ec_rxn_ids_tunable: + apply_kcat_constraints(model, update_rxns=ec_rxn_ids_tunable) _reset_solver_basis(model) flux_sims = None diff --git a/src/geckopy/kcat_tuning/sensitivity_tuning.py b/src/geckopy/kcat_tuning/sensitivity_tuning.py index 386567d..da903a5 100644 --- a/src/geckopy/kcat_tuning/sensitivity_tuning.py +++ b/src/geckopy/kcat_tuning/sensitivity_tuning.py @@ -25,16 +25,19 @@ from __future__ import annotations import logging +import math from dataclasses import dataclass, field from typing import TYPE_CHECKING, Optional, Sequence import numpy as np from ..ec_model.constants import ( + POOL_ID, PROT_PREFIX, USAGE_PREFIX, ) from ..ec_model.pipeline.apply_kcat import apply_kcat_constraints +from ..ec_model.pipeline.populate_ec import split_light_rxn_id from .corrections import Correction, corrections as build_corrections if TYPE_CHECKING: @@ -49,9 +52,20 @@ # than this (absolute, 1/h). Using a tolerance instead of exact equality # avoids spinning on sub-numerical-noise increases. _GROWTH_STALL_TOL = 1e-9 -# Physical upper bound for a tuned kcat (Bar-Even et al. 2011); tuning never -# pushes a kcat above the diffusion limit, no matter how limiting it looks. -_MAX_TUNED_KCAT = 1e7 + + +def _num2str(value: float) -> str: + """A scalar as MATLAB's ``num2str`` prints it: integers as such, other + values with ``max(floor(log10(abs(x))) + 5, 5)`` significant digits.""" + value = float(value) + if math.isnan(value): + return "NaN" + if math.isinf(value): + return "Inf" if value > 0 else "-Inf" + if value == int(value): + return str(int(value)) + digits = max(math.floor(math.log10(abs(value))) + 5, 5) + return f"{value:.{digits}g}" @dataclass @@ -61,7 +75,8 @@ class TunedKcatsResult: Attributes ---------- rxns - IDs of the reactions whose kcat was bumped, in tuning order. + IDs of the reactions whose kcat was bumped, in model order (the + order of ``model.reactions``), as ``sensitivityTuning.m`` lists them. rxn_names Reaction names parallel to ``rxns`` (falls back to the reaction ID when it has no name). @@ -69,12 +84,13 @@ class TunedKcatsResult: For each tuned reaction, the catalysing enzyme ID(s), ``;``-joined when more than one enzyme is involved. old_kcat - kcat value before tuning, parallel to ``rxns``. + kcat value of the input model, parallel to ``rxns``. A reaction that + an earlier run had already tuned reports that run's value. new_kcat kcat value after tuning, parallel to ``rxns``. source - Original kcat provenance (``model.ec.source``) before - tuning, parallel to ``rxns``. + kcat provenance (``model.ec.source``) of the input model, + parallel to ``rxns``. """ rxns: list[str] = field(default_factory=list) @@ -135,11 +151,17 @@ def sensitivity_tuning( is in ``prot_to_ignore``. 5. Finds the catalysed reaction consuming the most of that enzyme's pseudometabolite and multiplies that reaction's - ``ec.kcat`` by ``fold_change``. - 6. Annotates ``ec.notes`` with the original kcat and source - (unless the source is already ``"sensitivityTuning"``). + ``ec.kcat`` by ``fold_change``. A tuned kcat has no upper limit. + 6. Annotates ``ec.notes`` with the original kcat (printed as MATLAB's + ``num2str`` does) and source, unless the source is already + ``"sensitivityTuning"``. 7. Re-applies kcat constraints for that single reaction. + On a gecko-light model each iteration instead picks the reaction with + the most flux through the protein pool, multiplies the kcat of all its + isozyme rows by ``fold_change`` and re-applies the constraint of that + reaction; ``ec.notes`` and ``ec.source`` are left as they are. + Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/sensitivityTuning.m. @@ -155,10 +177,13 @@ def sensitivity_tuning( Per-iteration kcat multiplier. prot_to_ignore UniProt IDs whose usage rxns are excluded from being picked - as the limiting enzyme. + as the limiting enzyme, matched without regard to case. An ID with + no usage reaction raises ``ValueError`` on a full model; on a + gecko-light model it is ignored. max_iterations - Hard cap on tuning iterations. A backstop against a target - growth rate that is unreachable by raising kcats alone. + Hard cap on tuning iterations (``sensitivityTuning.m`` has none). A + backstop against a target growth rate that is unreachable by + raising kcats alone. bio_rxn Biomass reaction id to set as the objective. Defaults to ``params.bio_rxn`` from the adapter when one is attached; if @@ -173,16 +198,11 @@ def sensitivity_tuning( Raises ------ ValueError - If the initial LP is infeasible, or if - ``desired_growth_rate`` is needed but no adapter / value is - available. - NotImplementedError - If ``model.ec.gecko_light`` is True. + If the initial LP is infeasible, if ``desired_growth_rate`` is + needed but no adapter / value is available, or if a + ``prot_to_ignore`` ID has no usage reaction (full model). """ - if model.ec.gecko_light: - raise NotImplementedError( - "sensitivity_tuning does not yet support gecko-light models." - ) + light = bool(model.ec.gecko_light) if desired_growth_rate is None: from ..adapter import resolve_adapter @@ -218,17 +238,25 @@ def sensitivity_tuning( "exchange constraints." ) + # The result reports kcat and source as the input model has them. + input_kcat = np.array(model.ec.kcat, dtype=float) + input_source = list(model.ec.source) + usage_rxn_ids = [ r.id for r in model.reactions if r.id.startswith(USAGE_PREFIX) ] - ignore_usage_ids = {f"{USAGE_PREFIX}{p}" for p in prot_to_ignore} + if light: + ignore_usage_ids: set[str] = set() + ignore_rxn_ids = _light_rxns_of_enzymes(model, prot_to_ignore) + else: + ignore_usage_ids = _usage_ids_to_ignore(usage_rxn_ids, prot_to_ignore) + ignore_rxn_ids = set() last_growth = float("nan") - tuned_ec_indices: list[int] = [] # indices into model.ec.rxns - # ec_idx -> (pre-tune kcat, pre-tune source), captured the first time a - # reaction is tuned so the result doesn't have to re-parse ec.notes. - pretune: dict[int, tuple[float, str]] = {} + # Indices into model.ec.rxns; on a light model the first isozyme row of + # each tuned reaction. + tuned_ec_indices: list[int] = [] iteration = 1 while iteration <= max_iterations: @@ -252,6 +280,22 @@ def sensitivity_tuning( iteration += 1 + if light: + step = _light_step(model, sol, ignore_rxn_ids) + if step is None: + logger.warning( + "sensitivity_tuning: no reaction draws on the protein " + "pool; cannot continue tuning." + ) + break + ec_rows = step + for i in ec_rows: + model.ec.kcat[i] = float(model.ec.kcat[i]) * fold_change + apply_kcat_constraints(model, update_rxns=[model.ec.rxns[ec_rows[0]]]) + if ec_rows[0] not in tuned_ec_indices: + tuned_ec_indices.append(ec_rows[0]) + continue + # Pick usage rxn with max flux (= most produced enzyme), # excluding any in prot_to_ignore. best_usage_id: Optional[str] = None @@ -318,27 +362,16 @@ def sensitivity_tuning( # Record the pre-tune kcat/source the first time this reaction is # tuned, and annotate notes (idempotently). if model.ec.source[ec_idx] != _TUNING_SOURCE: - pretune.setdefault( - ec_idx, - (float(model.ec.kcat[ec_idx]), model.ec.source[ec_idx]), - ) old_note = model.ec.notes[ec_idx] new_note = ( - f"preTuneKcat={model.ec.kcat[ec_idx]} | " + f"preTuneKcat={_num2str(model.ec.kcat[ec_idx])} | " f"source:{model.ec.source[ec_idx]}" ) model.ec.notes[ec_idx] = ( f"{old_note}; {new_note}" if old_note else new_note ) - target_kcat = float(model.ec.kcat[ec_idx]) * fold_change - if target_kcat > _MAX_TUNED_KCAT: - target_kcat = _MAX_TUNED_KCAT - logger.warning( - "sensitivity_tuning: kcat for %s capped at the physical " - "ceiling %g /s.", target_rxn.id, _MAX_TUNED_KCAT, - ) - model.ec.kcat[ec_idx] = target_kcat + model.ec.kcat[ec_idx] = float(model.ec.kcat[ec_idx]) * fold_change model.ec.source[ec_idx] = _TUNING_SOURCE apply_kcat_constraints(model, update_rxns=[target_rxn.id]) @@ -351,33 +384,95 @@ def sensitivity_tuning( max_iterations, desired_growth_rate, ) - return _build_result(model, tuned_ec_indices, pretune) + return _build_result(model, tuned_ec_indices, input_kcat, input_source) + + +def _usage_ids_to_ignore( + usage_rxn_ids: list[str], prot_to_ignore: list[str], +) -> set[str]: + """Usage reaction ids of ``prot_to_ignore``, matched without regard to + case; an ID with no usage reaction raises, as MATLAB's ``find`` does.""" + ignore: set[str] = set() + for prot in prot_to_ignore: + wanted = f"{USAGE_PREFIX}{prot}".lower() + hits = [uid for uid in usage_rxn_ids if uid.lower() == wanted] + if not hits: + raise ValueError( + f"prot_to_ignore: no usage reaction for protein {prot!r}." + ) + ignore.update(hits) + return ignore + + +def _light_rxns_of_enzymes( + model: "EcModel", prot_to_ignore: list[str], +) -> set[str]: + """Cobra reaction ids that any of the ``prot_to_ignore`` enzymes catalyse + (gecko-light), matched without regard to case; unknown IDs are ignored.""" + ignore: set[str] = set() + if not prot_to_ignore: + return ignore + mat = model.ec.rxn_enz_mat.tocsc() + wanted = {p.lower() for p in prot_to_ignore} + for j, enzyme in enumerate(model.ec.enzymes): + if enzyme.lower() not in wanted: + continue + for row in mat[:, j].nonzero()[0]: + ignore.add(split_light_rxn_id(model.ec.rxns[int(row)])[1]) + return ignore + + +def _light_step( + model: "EcModel", sol, ignore_rxn_ids: set[str], +) -> Optional[list[int]]: + """The ec rows of the reaction with the most flux through the protein + pool (gecko-light), or ``None`` when no reaction draws on it.""" + pool = model.metabolites.get_by_id(POOL_ID) + best_id: Optional[str] = None + best_value = 0.0 + for rxn in model.reactions: + if rxn.id in ignore_rxn_ids or pool not in rxn.metabolites: + continue + value = rxn.metabolites[pool] * float(sol.fluxes.get(rxn.id, 0.0)) + if value < best_value: + best_value = value + best_id = rxn.id + if best_id is None: + return None + rows = [ + i for i, ec_id in enumerate(model.ec.rxns) + if split_light_rxn_id(ec_id)[1] == best_id + ] + return rows or None def _build_result( model: "EcModel", ec_indices: list[int], - pretune: dict[int, tuple[float, str]], + input_kcat: np.ndarray, + input_source: list[str], ) -> TunedKcatsResult: - """Snapshot pre/post kcat values for tuned reactions.""" + """Snapshot input and tuned kcat values of the tuned reactions, listed in + model order.""" if not ec_indices: return TunedKcatsResult() - rxns = [model.ec.rxns[i] for i in ec_indices] - rxn_names = [ - model.reactions.get_by_id(r).name or r for r in rxns + light = bool(model.ec.gecko_light) + cobra_ids = [ + split_light_rxn_id(model.ec.rxns[i])[1] if light else model.ec.rxns[i] + for i in ec_indices ] + order = sorted( + range(len(ec_indices)), key=lambda k: model.reactions.index(cobra_ids[k]), + ) + ec_indices = [ec_indices[k] for k in order] + rxns = [cobra_ids[k] for k in order] + rxn_names = [model.reactions.get_by_id(r).name or r for r in rxns] new_kcat = np.array( [float(model.ec.kcat[i]) for i in ec_indices], dtype=float, ) - - # Pre-tune kcat/source come from the dict captured during tuning. - old_kcat: list[float] = [] - sources: list[str] = [] - for i in ec_indices: - old, src = pretune.get(i, (float("nan"), "")) - old_kcat.append(old) - sources.append(src) + old_kcat = np.array([float(input_kcat[i]) for i in ec_indices], dtype=float) + sources = [input_source[i] for i in ec_indices] enzymes: list[str] = [] rxn_enz_dense = model.ec.rxn_enz_mat.toarray() @@ -389,7 +484,7 @@ def _build_result( rxns=rxns, rxn_names=rxn_names, enzymes=enzymes, - old_kcat=np.array(old_kcat, dtype=float), + old_kcat=old_kcat, new_kcat=new_kcat, source=sources, ) diff --git a/src/geckopy/kcat_tuning/sigma_fitter.py b/src/geckopy/kcat_tuning/sigma_fitter.py index 566778f..d7fe490 100644 --- a/src/geckopy/kcat_tuning/sigma_fitter.py +++ b/src/geckopy/kcat_tuning/sigma_fitter.py @@ -9,7 +9,9 @@ The algorithm: try a range of sigma values, set the protein-pool bound at each, solve the model, compare the predicted growth rate to the experimental one, pick the sigma that minimises the -difference. +difference. The default scan is the 0.01 grid of GECKO's +``sigmaFitter``; the objective reaction's bounds are opened to +``[0, 1000]`` first, and stay open on the fitted model. Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/sigmaFitter.m. @@ -29,6 +31,9 @@ # Sigma resolution the bisection converges to (fine enough to resolve a steep # growth curve; ~20 iterations from the [0, 1] bracket). _SIGMA_BISECT_TOL = 1e-6 +# Relative tolerance within which growth counts as being on its plateau, so +# solver noise does not move the start of the plateau. +_PLATEAU_RTOL = 1e-9 @dataclass @@ -63,11 +68,18 @@ class SigmaFitterResult: def _growth_at_sigma( model: "EcModel", p_tot: float, f: float, sigma: float, ) -> float: - """Solve the model with the pool sized for ``sigma``; revert after.""" + """Solve the model with the pool sized for ``sigma``; revert after. + + A solve that is not optimal counts as growth 0, as an empty solution + does in ``sigmaFitter.m``. + """ with model: set_prot_pool_size(model, p_tot=p_tot, f=f, sigma=float(sigma)) sol = model.optimize() - return float(sol.objective_value or 0.0) + if sol.status != "optimal" or sol.objective_value is None: + return 0.0 + value = float(sol.objective_value) + return value if np.isfinite(value) else 0.0 def _sigma_error_pct(growth: float, growth_rate: float) -> float: @@ -83,17 +95,28 @@ def fit_sigma( p_tot: Optional[float] = None, f: Optional[float] = None, n_sigma_steps: int = 100, - method: Literal["grid", "bisect"] = "bisect", + method: Literal["grid", "bisect"] = "grid", ) -> SigmaFitterResult: """Sweep sigma and pick the one matching ``growth_rate`` best. - Growth is monotone non-decreasing in sigma, so ``method="bisect"`` - (default) finds the best sigma in about ``log2(n)`` solves; the - returned grids then contain only the sigmas actually evaluated. - ``method="grid"`` keeps the legacy ``[1/n, 2/n, ..., 1.0]`` sweep - for diagnostic plotting (every grid point evaluated). - - The best sigma is re-applied to ``model`` before return. + ``method="grid"`` (default) is the sweep of GECKO's ``sigmaFitter``: + sigma ``i / n`` for ``i = 1..n`` (every grid point evaluated), and the + first grid point with the smallest error wins. For a target that growth + cannot reach this is the first sigma of the growth plateau when the + plateau values are identical; values that differ by solver noise can move + the pick to a later grid point, as in ``sigmaFitter.m``. A model that + cannot grow at any sigma gives the first grid point, ``1 / n_sigma_steps``. + ``method="bisect"`` exploits that growth is monotone non-decreasing in + sigma and finds the best sigma, to a resolution of 1e-6 and not on the + grid, in about ``log2(n)`` solves; the returned grids then contain only + the sigmas actually evaluated. For an unreachable target it finds where + the growth plateau starts (growth within a relative 1e-9 of its value at + sigma 1); a model that cannot grow at sigma 1 gives ``1 / n_sigma_steps``, + as the grid does. + + The bounds of the objective reaction are set to ``[0, 1000]`` before the + scan, as ``sigmaFitter.m`` does, and stay so on ``model``. The best sigma + is re-applied to ``model`` before return. Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/sigmaFitter.m. @@ -103,8 +126,9 @@ def fit_sigma( model EcModel with the protein pool machinery installed and the objective set (typically the biomass reaction). Mutated in - place: the optimal sigma is applied via - ``set_prot_pool_size`` before return. + place: the objective reaction's bounds become ``[0, 1000]`` and + the optimal sigma is applied via ``set_prot_pool_size`` before + return. growth_rate Experimental growth rate to match. Defaults to ``model.adapter.params.gr_exp``. @@ -120,9 +144,9 @@ def fit_sigma( it caps the number of bisection iterations (the default 100 is far more than the ~20 needed to converge). Default 100. method - ``"bisect"`` (default) finds the best sigma in ~``log2(n)`` - solves using the monotonicity of growth in sigma; ``"grid"`` - sweeps the full grid. + ``"grid"`` (default) sweeps the full grid as ``sigmaFitter.m`` + does; ``"bisect"`` finds the best sigma in ~``log2(n)`` solves using + the monotonicity of growth in sigma. Returns ------- @@ -149,6 +173,11 @@ def fit_sigma( ) params = model.adapter.params if model.adapter is not None else None + # sigmaFitter.m opens the bounds of the objective for good. + for rxn in model.reactions: + if rxn.objective_coefficient != 0: + rxn.bounds = (0.0, 1000.0) + if growth_rate is None: growth_rate = float(params.gr_exp) if p_tot is None: @@ -156,6 +185,7 @@ def fit_sigma( if f is None: f = float(params.f) + plateau_sigma: Optional[float] = None if method == "grid": sigma_grid = np.array( [(i + 1) / n_sigma_steps for i in range(n_sigma_steps)], @@ -174,23 +204,34 @@ def growth_at(sigma: float) -> float: evaluated[sigma] = _growth_at_sigma(model, p_tot, f, sigma) return evaluated[sigma] - if growth_at(1.0) < growth_rate: - # Target unreachable even at the full budget; closest is sigma=1. - best_sigma = 1.0 + g_max = growth_at(1.0) + if g_max < growth_rate: + # Target unreachable even at the full budget: the closest growth + # is the plateau, so look for the smallest sigma that reaches it + # (to solver noise). + reach = g_max - _PLATEAU_RTOL * abs(g_max) + else: + reach = growth_rate + # Bisect for the smallest sigma reaching ``reach``. The tolerance is + # fine (not 1/n) so a steep growth curve still resolves the + # crossing; n_sigma_steps just caps iterations. + lo, hi = 0.0, 1.0 + if g_max <= 0: + # No growth at any sigma: every error is equal, and the grid takes + # its first point. + hi = 1.0 / n_sigma_steps else: - # Bisect for the smallest sigma reaching the target growth. The - # tolerance is fine (not 1/n) so a steep growth curve still - # resolves the crossing; n_sigma_steps just caps iterations. - lo, hi = 0.0, 1.0 for _ in range(n_sigma_steps): if hi - lo <= _SIGMA_BISECT_TOL: break mid = 0.5 * (lo + hi) - if growth_at(mid) >= growth_rate: + if growth_at(mid) >= reach: hi = mid else: lo = mid - best_sigma = hi + best_sigma = hi + if g_max < growth_rate: + plateau_sigma = hi growth_at(best_sigma) sigma_grid = np.array(sorted(evaluated), dtype=float) growth_grid = np.array( @@ -200,7 +241,10 @@ def growth_at(sigma: float) -> float: error_grid = np.array( [_sigma_error_pct(g, growth_rate) for g in growth_grid], dtype=float, ) - best_sigma = float(sigma_grid[int(np.argmin(error_grid))]) + if plateau_sigma is not None: + best_sigma = float(plateau_sigma) + else: + best_sigma = float(sigma_grid[int(np.argmin(error_grid))]) # Apply the best sigma to the model permanently (geckopy divergence). set_prot_pool_size(model, p_tot=p_tot, f=f, sigma=best_sigma) diff --git a/src/geckopy/kcat_tuning/truncate_values.py b/src/geckopy/kcat_tuning/truncate_values.py index cbc33d0..42ff155 100644 --- a/src/geckopy/kcat_tuning/truncate_values.py +++ b/src/geckopy/kcat_tuning/truncate_values.py @@ -5,6 +5,7 @@ """ from __future__ import annotations +from decimal import ROUND_HALF_UP, Decimal from typing import Optional import numpy as np @@ -21,6 +22,10 @@ def truncate_values( order_magn = max(ceil(log10(|v|)), 0) result = round(v, 6 - order_magn) + A value exactly halfway between two candidates rounds away from + zero, as MATLAB's ``round(x, n)`` does (``100000.5`` gives + ``100001`` and ``1234565`` gives ``1234570``). The value is taken + at its shortest decimal representation, so ``2.675`` is a tie. Zero, NaN, and infinite values are returned unchanged. For 2D input, only the columns selected by ``columns`` are @@ -69,9 +74,11 @@ def truncate_values( def _truncate_scalar(v: float) -> float: """Round ``v`` to ~6 significant figures. - Zero, NaN, and infinite values are returned unchanged. + Ties round away from zero. Zero, NaN, and infinite values are + returned unchanged. """ if v == 0 or np.isnan(v) or np.isinf(v): return v order_magn = max(int(np.ceil(np.log10(abs(v)))), 0) - return round(float(v), 6 - order_magn) + quantum = Decimal(1).scaleb(order_magn - 6) + return float(Decimal(repr(float(v))).quantize(quantum, rounding=ROUND_HALF_UP)) diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 4905278..14a7398 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -282,6 +282,10 @@ def flexibilize_enz_concs( ) if pool_pred_growth is not None: pred_growth = pool_pred_growth + # No enzyme limits growth any more: the highest growth reached is the + # target of the final pass and of the reported model, whether or not + # exp_growth was reached. + exp_growth = pred_growth break if families is not None: @@ -492,11 +496,15 @@ def _unit_scores( candidates.append((k, reactions, active)) if _solver_supports_duals(model): reduced = sol.reduced_costs - for k, _, active in candidates: - scores[k] = sum(abs(float(reduced[r.id])) * r.upper_bound for r in active) - scores[scores < 1e-12] = 0.0 - if (scores > 0).any(): - return scores + needed = [r.id for _, _, active in candidates for r in active] + # cobra fills reduced_costs with NaN when the solver exposes no duals + if np.isfinite(np.asarray([float(reduced[rid]) for rid in needed])).all(): + for k, _, active in candidates: + scores[k] = sum(abs(float(reduced[r.id])) * r.upper_bound for r in active) + scores[scores < 1e-12] = 0.0 + if (scores > 0).any(): + return scores + scores = np.zeros(len(units)) # No dual singles a unit out (none exposed, or the protein pool binds as # well and carries the whole price): probe each candidate with the pool # lifted, so the probe finds the caps that limit once the pool is relaxed. diff --git a/src/geckopy/limit_proteins/get_conc_control_coeffs.py b/src/geckopy/limit_proteins/get_conc_control_coeffs.py index 1608c30..012ac44 100644 --- a/src/geckopy/limit_proteins/get_conc_control_coeffs.py +++ b/src/geckopy/limit_proteins/get_conc_control_coeffs.py @@ -20,9 +20,22 @@ _GROWTH_DELTA_THRESHOLD = 1e-10 -# Solvers we have explicitly confirmed do not expose LP duals via optlang. -# Anything else is assumed to (true for glpk, glpk_exact, gurobi, cplex). -_NO_DUALS_SOLVER_TAGS = ("scipy",) +# Relative raise of a usage upper bound that confirms a nonzero reduced cost (0.01 percent). +_VERIFY_STEP = 1e-4 + +# A reduced cost is kept when the optimum rises by at least this fraction of the rise the +# reduced cost predicts for the step (reduced cost times the bound times ``_VERIFY_STEP``). +_VERIFY_MIN_FRACTION = 0.5 + +# Smallest rise of the optimum that counts, relative to max(1, |optimum|); below it the +# two solves do not differ beyond solver noise. +_VERIFY_GROWTH_FLOOR = 1e-12 + +# Solvers we have explicitly confirmed do not expose LP duals at all. Empty: optlang's scipy +# backend raises on ``dual`` but keeps HiGHS's bound marginals (see +# :func:`_scipy_upper_marginal`), and any solver whose dual read fails or returns NaN falls +# back to the finite-difference loop anyway. +_NO_DUALS_SOLVER_TAGS: tuple[str, ...] = () # One-shot info reporting when the shadow-price path is unavailable. _FALLBACK_REPORTED: set[str] = set() @@ -42,13 +55,26 @@ def get_conc_control_coeffs( ``limit`` times their upper bound are skipped (no sensitivity to expect there). - The default path uses the ``prot_`` metabolite **shadow price** - from a single LP solve — equivalent to the per-protein - finite-difference probe at the margin, but ``n`` times fewer LPs. - When the active solver is one that does not expose LP duals (notably - optlang's scipy backend), the function automatically falls back to - the per-protein finite-difference loop and logs once that it has - done so. The fold-change probe size is only used in that fallback. + The default path reads, from a single LP solve, the reduced cost of + each analysed usage reaction's forward variable. That is the local + marginal d(objective)/d(usage_ub) of the enzyme's own bound: it is + zero for an enzyme whose bound does not bind, also when the shared + protein pool does, because the pool's price is part of the + ``prot_`` metabolite dual and cancels out of the reduced cost + (``|dual(prot_) - dual(prot_pool)|``). A nonzero reduced cost is + kept only if raising that bound by 0.01 percent raises the optimum + by at least half of the rise the reduced cost predicts (one extra + solve per nonzero value), because the dual overestimates the effect + at a degenerate vertex. The test is relative to the predicted rise, + so a small coefficient on a small bound is kept as readily as a large + one; the rise must also exceed 1e-12 times max(1, |optimum|). A zero + reduced cost needs no check. With optlang's scipy backend (HiGHS) the reduced cost is + read from the bound marginal of the last ``linprog`` result, since optlang does not + expose variable duals there. When the active solver gives no usable variable duals + (a missing or NaN value, or a scipy result without bound marginals) the function falls + back to the per-protein finite-difference loop, which needs one solve per protein + (several minutes on a genome-scale ecModel), and logs once that it has done so. The + fold-change probe size is only used in that fallback. Ported from GECKO MATLAB: src/geckomat/limit_proteins/getConcControlCoeffs.m. @@ -86,15 +112,20 @@ def get_conc_control_coeffs( return np.zeros(0, dtype=bool), np.zeros(0, dtype=float) if _solver_supports_duals(model): - return _shadow_price_coeffs(model, proteins, limit) - _report_finite_difference_fallback(model) + return _shadow_price_coeffs(model, proteins, limit, fold_change) + _report_finite_difference_fallback(model, "does not expose LP duals") return _finite_difference_coeffs(model, proteins, fold_change, limit) def _shadow_price_coeffs( model: "EcModel", proteins: list[str], limit: float, + fold_change: float = 2.0, ) -> tuple[np.ndarray, np.ndarray]: - """One-solve control coefficients via the prot_ metabolite duals.""" + """One-solve control coefficients from the usage reactions' reduced costs. + + Falls back to :func:`_finite_difference_coeffs` when the solver gives + no usable duals. + """ n = len(proteins) enz_mask = np.zeros(n, dtype=bool) coeffs = np.zeros(n, dtype=float) @@ -102,9 +133,12 @@ def _shadow_price_coeffs( sol = model.optimize() if not _solution_is_optimal(sol): return enz_mask, coeffs + initial_growth = sol.objective_value fluxes = sol.fluxes cobra_rxn_ids = {r.id for r in model.reactions} + reduced_costs: dict[int, float] = {} + pinned_at_lower: set[int] = set() for i, protein in enumerate(proteins): rxn_id = f"{USAGE_PREFIX}{protein}" if rxn_id not in cobra_rxn_ids: @@ -117,16 +151,88 @@ def _shadow_price_coeffs( if usage_flux / prev_ub <= limit: continue enz_mask[i] = True - # |shadow price| of the prot_ metabolite constraint equals the - # local marginal d(objective)/d(usage_ub). Verified in tests across - # several scenarios and solvers (glpk, glpk_exact, gurobi). - sp = abs(model.enzymes.get_by_id(protein).shadow_price) - if sp > _GROWTH_DELTA_THRESHOLD: - coeffs[i] = sp + rc = _forward_reduced_cost(rxn, usage_flux) + if rc is None: + _report_finite_difference_fallback(model, "returns no usable variable duals") + return _finite_difference_coeffs(model, proteins, fold_change, limit) + reduced_costs[i] = rc + if rxn.lower_bound > 0 and usage_flux <= rxn.lower_bound * (1.0 + 1e-9): + # Held at its lower bound: raising the upper bound changes nothing. + pinned_at_lower.add(i) + + for i, rc in reduced_costs.items(): + if rc == 0.0 or i in pinned_at_lower: + continue + rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{proteins[i]}") + predicted_gain = rc * rxn.upper_bound * _VERIFY_STEP + with model: + rxn.upper_bound = rxn.upper_bound * (1.0 + _VERIFY_STEP) + new_sol = model.optimize() + if not _solution_is_optimal(new_sol): + continue + gain = new_sol.objective_value - initial_growth + if gain >= max( + _VERIFY_MIN_FRACTION * predicted_gain, + _VERIFY_GROWTH_FLOOR * max(1.0, abs(initial_growth)), + ): + coeffs[i] = rc return enz_mask, coeffs +def _forward_reduced_cost(rxn, flux: Optional[float] = None) -> Optional[float]: + """|reduced cost| of the reaction's forward variable, or None if unavailable. + + ``flux`` is the reaction's flux in the solution the reduced cost must belong to; the + scipy path needs it to confirm that the stored result is that solution. + + ``reaction.reduced_cost`` is the forward variable's value minus the reverse + variable's, and the reverse variable of a usage reaction carries the opposite + sign, so it is twice the sensitivity to the upper bound; the forward variable + alone is that sensitivity. + """ + try: + value = float(rxn.forward_variable.dual) + except NotImplementedError: + value = _scipy_upper_marginal(rxn, flux) + if value is None: + return None + except (AttributeError, TypeError, ValueError, RuntimeError): + return None + if not np.isfinite(value): + return None + return abs(value) + + +def _scipy_upper_marginal(rxn, flux: Optional[float]) -> Optional[float]: + """Marginal of the forward variable's upper bound in optlang's scipy backend, or None. + + optlang raises ``NotImplementedError`` for ``Variable.dual`` with the scipy solver but + keeps the whole ``linprog`` result in ``problem._solution``. For HiGHS that result holds + ``upper.marginals``: d(objective)/d(upper bound) of the minimised problem, in the order + of the variables handed to ``linprog``. ``cobra.Model.optimize`` resets optlang's own copy + of the primal values after the solve (it restores the objective direction), so the stored + result is checked against ``flux``, the flux cobra returned for this reaction: the forward + minus the reverse variable of the stored result must equal it. Without ``flux``, or when + the values differ, the result is not used (None) and the finite-difference loop runs. + """ + if flux is None: + return None + try: + problem = rxn.model.solver.problem + solution = problem._solution + marginals = solution.upper.marginals + names = list(problem.bounds) + forward = names.index(rxn.forward_variable.name) + reverse = names.index(rxn.reverse_variable.name) + stored = float(solution.x[forward]) - float(solution.x[reverse]) + if abs(stored - flux) > 1e-9 * max(1.0, abs(flux)): + return None + return float(marginals[forward]) + except (AttributeError, IndexError, ValueError, TypeError, KeyError): + return None + + def _finite_difference_coeffs( model: "EcModel", proteins: list[str], @@ -135,8 +241,8 @@ def _finite_difference_coeffs( ) -> tuple[np.ndarray, np.ndarray]: """Per-protein finite-difference fallback (one solve per protein). - Same contract as :func:`get_conc_control_coeffs`. Kept as the - fallback for solvers that do not expose LP duals. + Same contract as :func:`get_conc_control_coeffs`. Used when the solver + does not expose usable variable duals. """ n = len(proteins) enz_mask = np.zeros(n, dtype=bool) @@ -193,26 +299,26 @@ def _solution_is_optimal(sol) -> bool: def _solver_supports_duals(model: "EcModel") -> bool: - """True iff the model's LP solver exposes constraint duals. + """True unless the solver is listed in ``_NO_DUALS_SOLVER_TAGS``. - optlang's scipy backend explicitly raises ``NotImplementedError`` on - ``constraint.dual``; glpk/glpk_exact/gurobi/cplex all expose duals - for LP optima. We detect by the solver module name (cheap; no probe - solve) and assume "supports duals" unless we know otherwise. + glpk/glpk_exact/gurobi/cplex expose variable duals for LP optima and optlang's scipy + backend exposes them through :func:`_scipy_upper_marginal`. A solver whose dual read + fails is detected per value in :func:`_forward_reduced_cost`, which sends the call to the + finite-difference loop. """ module = type(model.solver).__module__.lower() return not any(tag in module for tag in _NO_DUALS_SOLVER_TAGS) -def _report_finite_difference_fallback(model: "EcModel") -> None: +def _report_finite_difference_fallback(model: "EcModel", reason: str) -> None: """Log once per solver that we fell back to the per-protein loop.""" module = type(model.solver).__module__ if module in _FALLBACK_REPORTED: return _FALLBACK_REPORTED.add(module) logger.info( - "get_conc_control_coeffs: solver %r does not expose LP duals; " + "get_conc_control_coeffs: solver %r %s; " "using the per-protein finite-difference fallback (one solve per " - "protein, %sx slower). This is reported once per process.", - module, "n", + "protein, nx slower). This is reported once per process.", + module, reason, ) diff --git a/src/geckopy/utilities/enzyme_usage.py b/src/geckopy/utilities/enzyme_usage.py index 5f0c60f..328d462 100644 --- a/src/geckopy/utilities/enzyme_usage.py +++ b/src/geckopy/utilities/enzyme_usage.py @@ -65,7 +65,8 @@ def enzyme_usage( * ``abs_usage = |fluxes[usage_prot_]|`` (mg/gDCW). * ``ub`` = the usage reaction's upper bound. - * ``cap_usage = abs_usage / ub`` (or 0 if ``ub == 0``). + * ``cap_usage = abs_usage / ub`` (NaN if ``ub == 0``, as + :attr:`Enzyme.cap_usage` and GECKO MATLAB's ``getEnzymeUsage`` give). Ported from GECKO MATLAB: src/geckomat/utilities/enzymeUsage.m. @@ -111,7 +112,7 @@ def enzyme_usage( ub_arr = np.array(ub_list, dtype=float) with np.errstate(divide="ignore", invalid="ignore"): cap_arr = np.where( - ub_arr > 0, abs_usage_arr / ub_arr, 0.0, + ub_arr != 0, abs_usage_arr / ub_arr, np.nan, ) if not include_zero: diff --git a/tests/test_apply_kcat.py b/tests/test_apply_kcat.py index 5a0c7f8..73e87f5 100644 --- a/tests/test_apply_kcat.py +++ b/tests/test_apply_kcat.py @@ -182,14 +182,14 @@ def test_update_rxns_unknown_id_raises(): apply_kcat_constraints(ec_model, update_rxns=["not_a_rxn"]) -def test_update_rxns_empty_list_is_noop(): +def test_update_rxns_empty_list_updates_all(): ec_model = _ectestgem_ec_model() ec_model.ec.kcat[ec_model.ec.rxns.index("R3")] = 10.0 apply_kcat_constraints(ec_model, update_rxns=[]) r3_coef = _get_s_coef(ec_model.reactions.get_by_id("R3"), "prot_P4") - assert r3_coef == 0.0 + assert r3_coef < 0.0 def test_setting_kcat_to_zero_and_reapplying_clears_old_constraint(): diff --git a/tests/test_conc_control_coeffs_shadow.py b/tests/test_conc_control_coeffs_shadow.py new file mode 100644 index 0000000..2eca72e --- /dev/null +++ b/tests/test_conc_control_coeffs_shadow.py @@ -0,0 +1,127 @@ +"""get_conc_control_coeffs: reduced-cost path across solvers, and its fallbacks. + +The coefficient of an enzyme is the reduced cost of its usage reaction's forward variable, +kept only where a 0.01 percent rise of the bound raises the optimum by at least half of +the rise the reduced cost predicts. An enzyme whose own bound does not bind has coefficient +0 also while the shared protein pool binds. +""" +import sys +from types import SimpleNamespace + +import numpy as np +import pytest + +from geckopy.limit_proteins import get_conc_control_coeffs +from geckopy.limit_proteins.get_conc_control_coeffs import ( + _finite_difference_coeffs, + _forward_reduced_cost, +) + +from test_get_conc_control_coeffs import _build_toy_ec_model + +# The toy models have finite bounds only. glpk_exact aborts the process on a model with an +# infinite bound, which every ecModel built by geckopy has, whatever this function does. +SOLVERS = ["glpk", "glpk_exact", "gurobi", "scipy"] + + +def _with_solver(model, solver): + try: + model.solver = solver + model.optimize() + except Exception as exc: # solver not installed or not licensed + pytest.skip(f"solver {solver} unavailable: {exc}") + return model + + +def _pool_binding_model(solver): + """Enzyme bound 5 does not bind; the protein pool caps growth at 3.""" + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + model.reactions.get_by_id("prot_pool_exchange").upper_bound = 3.0 + return model + + +def _degenerate_model(solver): + """The enzyme bound and R_AB's own bound both equal 5: raising the enzyme bound alone raises nothing.""" + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + model.reactions.get_by_id("R_AB").upper_bound = 5.0 + return model + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_pool_binding_gives_zero_for_an_enzyme_whose_bound_does_not_bind(solver): + model = _pool_binding_model(solver) + assert model.optimize().objective_value == pytest.approx(3.0) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == 0.0 + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_binding_enzyme_bound_gives_its_marginal_growth(solver): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_shadow_price_equals_the_local_finite_difference(solver): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0, enzyme_coeff_in_R=-2.0), solver) + _, coeffs = get_conc_control_coeffs(model) + _, local = _finite_difference_coeffs(model, ["E"], 1.001, 0.0) + assert coeffs[0] == pytest.approx(local[0], rel=1e-6) + assert coeffs[0] == pytest.approx(0.5, rel=1e-6) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_degenerate_vertex_gives_zero(solver): + model = _degenerate_model(solver) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == 0.0 + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_model_is_unchanged_after_the_call(solver): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + usage = model.reactions.get_by_id("usage_prot_E") + before = (usage.lower_bound, usage.upper_bound) + get_conc_control_coeffs(model) + assert (usage.lower_bound, usage.upper_bound) == before + + +# --------------------------------------------------------------------------- # +# Duals that are missing or not finite +# --------------------------------------------------------------------------- # + +class _Raises: + @property + def dual(self): + raise NotImplementedError + + +def test_forward_reduced_cost_is_none_without_a_usable_dual(): + assert _forward_reduced_cost(SimpleNamespace(forward_variable=_Raises())) is None + assert _forward_reduced_cost(SimpleNamespace(forward_variable=SimpleNamespace(dual=float("nan")))) is None + assert _forward_reduced_cost(SimpleNamespace(forward_variable=SimpleNamespace(dual=float("inf")))) is None + + +def test_forward_reduced_cost_is_the_absolute_value(): + assert _forward_reduced_cost(SimpleNamespace(forward_variable=SimpleNamespace(dual=-2.5))) == 2.5 + + +@pytest.mark.parametrize("solver", ["glpk", "glpk_exact", "gurobi"]) +def test_unusable_duals_fall_back_to_the_finite_difference(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + expected = _finite_difference_coeffs(model, ["E"], 2.0, 0.0) + module = sys.modules["geckopy.limit_proteins.get_conc_control_coeffs"] + monkeypatch.setattr(module, "_forward_reduced_cost", lambda rxn, flux=None: None) + enz, coeffs = get_conc_control_coeffs(model) + assert enz.tolist() == expected[0].tolist() + assert coeffs.tolist() == expected[1].tolist() + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) + + +def test_scipy_reads_the_bound_marginal(): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), "scipy") + _, coeffs = get_conc_control_coeffs(model, fold_change=2.0) + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) diff --git a/tests/test_enzyme_usage.py b/tests/test_enzyme_usage.py index bd84f4b..b6b2888 100644 --- a/tests/test_enzyme_usage.py +++ b/tests/test_enzyme_usage.py @@ -167,16 +167,16 @@ def test_missing_rxn_in_fluxes_treated_as_zero(): # --------------------------------------------------------------------------- # -# Zero upper bound -> cap_usage = 0 +# Zero upper bound -> cap_usage = NaN # --------------------------------------------------------------------------- # -def test_zero_ub_cap_usage_is_zero_not_inf(): +def test_zero_ub_cap_usage_is_nan(): model = _build_toy(enzyme_ub=0.0) sol = model.optimize() result = enzyme_usage(model, sol.fluxes) - # ub=0 -> no growth -> abs_usage=0; cap_usage must be 0, not inf or NaN. - assert result.cap_usage[0] == 0.0 - assert np.isfinite(result.cap_usage[0]) + # ub=0 -> no growth -> abs_usage=0; the ratio is undefined, as for Enzyme.cap_usage. + assert result.abs_usage[0] == 0.0 + assert np.isnan(result.cap_usage[0]) # --------------------------------------------------------------------------- # diff --git a/tests/test_get_conc_control_coeffs.py b/tests/test_get_conc_control_coeffs.py index dbae1b8..05454bb 100644 --- a/tests/test_get_conc_control_coeffs.py +++ b/tests/test_get_conc_control_coeffs.py @@ -245,7 +245,7 @@ def test_shadow_price_matches_finite_difference_coeff(ub, coeff_in_R, expected): # Both implementations live in the production module after the swap; tests # compare them as ground-truth helpers (the dispatcher uses shadow-price by -# default and falls back to finite-difference on scipy). +# default and falls back to finite-difference where a solver gives no usable duals). from geckopy.limit_proteins.get_conc_control_coeffs import ( # noqa: E402 _finite_difference_coeffs as _fd_coeffs, _shadow_price_coeffs as _sp_coeffs, @@ -392,24 +392,15 @@ def test_equiv_infeasible_model_both_zero(): [ "glpk", "glpk_exact", - # optlang's scipy backend explicitly does NOT implement duals - # (`constraint.dual` raises NotImplementedError), so the - # shadow-price approach is unavailable there. The finite-difference - # function still works on scipy. - pytest.param( - "scipy", - marks=pytest.mark.xfail( - raises=NotImplementedError, - reason="optlang scipy interface does not expose LP duals", - strict=True, - ), - ), + # optlang's scipy backend raises on variable duals; the reduced cost is read from + # HiGHS's bound marginal instead. + "scipy", "gurobi", ], ) def test_equiv_across_solvers(solver): - """The shadow-price/finite-difference equivalence holds for every LP - solver that exposes duals. Uses the single-enzyme linear-regime model + """The reduced-cost/finite-difference equivalence holds for every LP + solver. Uses the single-enzyme linear-regime model so the LP is non-degenerate / non-breakpoint (where duals would be non-unique).""" from cobra.util.solver import solvers @@ -444,7 +435,7 @@ def test_finite_diff_underestimates_when_step_crosses_breakpoint(): # --------------------------------------------------------------------------- # # Dispatcher: the public function defaults to shadow-price and falls back to -# finite-difference on solvers (scipy) that do not expose duals. +# finite-difference where the solver gives no usable variable duals. # --------------------------------------------------------------------------- # def test_dispatcher_uses_shadow_price_on_default_solver(): @@ -463,9 +454,10 @@ def test_dispatcher_uses_shadow_price_on_default_solver(): assert pc[0] != pytest.approx(fc[0]) -def test_dispatcher_falls_back_to_finite_diff_on_scipy(caplog): - """With the scipy solver (no LP duals), the dispatcher uses the - finite-difference path and reports the fallback once.""" +def test_dispatcher_falls_back_to_finite_diff_without_usable_duals(caplog, monkeypatch): + """When the solver gives no usable variable dual (here: the scipy bound + marginal is unavailable), the dispatcher uses the finite-difference path + and reports the fallback once.""" from cobra.util.solver import solvers if "scipy" not in solvers: pytest.skip("scipy not installed") @@ -478,6 +470,7 @@ def test_dispatcher_falls_back_to_finite_diff_on_scipy(caplog): ) # Reset the one-shot guard so the info log fires for this test. mod_ccc._FALLBACK_REPORTED.clear() + monkeypatch.setattr(mod_ccc, "_scipy_upper_marginal", lambda rxn, flux: None) model = _build_toy_ec_model(enzyme_ub=8.0) model.reactions.get_by_id("EX_A").lower_bound = -10.0 @@ -488,4 +481,4 @@ def test_dispatcher_falls_back_to_finite_diff_on_scipy(caplog): fm, fc = _fd_helper(model, fold_change=2.0) np.testing.assert_array_equal(pm, fm) np.testing.assert_allclose(pc, fc, rtol=1e-6, atol=1e-12) - assert "does not expose LP duals" in caplog.text + assert "returns no usable variable duals" in caplog.text diff --git a/tests/test_make_ec_model.py b/tests/test_make_ec_model.py index f4d06bb..88a2601 100644 --- a/tests/test_make_ec_model.py +++ b/tests/test_make_ec_model.py @@ -273,6 +273,7 @@ def test_make_ec_model_missing_uniprot_tsv_raises(tmp_path): adapter = ModelAdapter.from_folder(tmp_path) model = cobra.Model("m") + model.compartments = {"c": "c"} m = cobra.Metabolite("A", compartment="c") model.add_metabolites([m]) r = cobra.Reaction("r1") @@ -306,6 +307,7 @@ def test_make_ec_model_unmatched_genes_logged_and_annotated(tmp_path, caplog): adapter = ModelAdapter.from_folder(tmp_path) model = cobra.Model("m") + model.compartments = {"c": "c"} m_a = cobra.Metabolite("A", compartment="c") m_b = cobra.Metabolite("B", compartment="c") model.add_metabolites([m_a, m_b]) diff --git a/tests/test_parity_decisions_t2.py b/tests/test_parity_decisions_t2.py new file mode 100644 index 0000000..817002c --- /dev/null +++ b/tests/test_parity_decisions_t2.py @@ -0,0 +1,174 @@ +"""make_ec_model refuses what GECKO's makeEcModel refuses and treats reversibility and +exchange reactions as it does; set_kcat_for_reactions follows setKcatForReactions.""" +from pathlib import Path + +import cobra +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.ec_model.pipeline import set_kcat_for_reactions + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + + +def _inputs(): + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + model = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + return model, adapter + + +def _build(edit=None, *, light=False, adapter_edit=None): + model, adapter = _inputs() + if edit is not None: + edit(model) + if adapter_edit is not None: + adapter_edit(adapter) + return make_ec_model(model, adapter, gecko_light=light) + + +# --------------------------------------------------------------------------- # +# Refusals +# --------------------------------------------------------------------------- # + +def test_the_unedited_model_builds(): + assert _build().ec.n_rxns > 0 + assert _build(light=True).ec.n_rxns > 0 + + +def _rename_rxn(old, new): + return lambda model: setattr(model.reactions.get_by_id(old), "id", new) + + +@pytest.mark.parametrize("light", [False, True]) +@pytest.mark.parametrize("old,new", [ + ("R1", "R1_REV"), + ("R1", "R1_EXP_1"), + ("R5", "usage_prot_P5"), + ("R5", "prot_pool_R5"), +]) +def test_reaction_id_conflicts_are_refused(old, new, light): + with pytest.raises(ValueError, match="model.reactions"): + _build(_rename_rxn(old, new), light=light) + + +def test_a_reaction_id_with_the_suffix_text_inside_is_accepted(): + ec_model = _build(_rename_rxn("R1", "R_REVX")) + assert "R_REVX" in ec_model.reactions + + +@pytest.mark.parametrize("light", [False, True]) +def test_metabolite_id_starting_with_prot_is_refused(light): + def edit(model): + model.metabolites.get_by_id("m2c").id = "prot_m2" + with pytest.raises(ValueError, match="prot_"): + _build(edit, light=light) + + +def test_a_compartment_id_that_differs_from_its_name_is_refused(): + def edit(model): + model.compartments = {"c": "cytosol", "e": "extra"} + with pytest.raises(ValueError, match="compartment names"): + _build(edit) + + +def test_an_unknown_enzyme_compartment_is_refused(): + def adapter_edit(adapter): + adapter.params.enzyme_comp = "zz" + with pytest.raises(ValueError, match="compartment names"): + _build(adapter_edit=adapter_edit) + + +@pytest.mark.parametrize("light", [False, True]) +def test_a_model_without_genes_is_refused(light): + def edit(model): + for rxn in model.reactions: + rxn.gene_reaction_rule = "" + cobra.manipulation.delete.remove_genes( + model, list(model.genes), remove_reactions=False + ) + with pytest.raises(ValueError, match="no genes"): + _build(edit, light=light) + + +# --------------------------------------------------------------------------- # +# Reversibility and exchange reactions +# --------------------------------------------------------------------------- # + +@pytest.mark.parametrize("light", [False, True]) +def test_a_reaction_with_both_bounds_negative_is_not_split(light): + def edit(model): + model.reactions.get_by_id("R3").bounds = (-10.0, -1.0) + ec_model = _build(edit, light=light) + assert "R3_REV" not in ec_model.reactions + assert ec_model.reactions.get_by_id("R3").bounds == (-10.0, -1.0) + + +def _add_sink(model): + sink = cobra.Reaction("SINK2", lower_bound=-1000.0, upper_bound=1000.0) + model.add_reactions([sink]) + sink.add_metabolites({ + model.metabolites.get_by_id("m1c"): -1.0, + model.metabolites.get_by_id("m2c"): -1.0, + }) + + +@pytest.mark.parametrize("light", [False, True]) +def test_a_multi_metabolite_sink_is_an_exchange_reaction_and_stays_reversible(light): + ec_model = _build(_add_sink, light=light) + assert "SINK2_REV" not in ec_model.reactions + assert ec_model.reactions.get_by_id("SINK2").bounds == (-1000.0, 1000.0) + + +def test_a_reversible_internal_reaction_is_still_split(): + ec_model = _build() + assert "R2_REV_EXP_1" in ec_model.reactions + + +# --------------------------------------------------------------------------- # +# set_kcat_for_reactions +# --------------------------------------------------------------------------- # + +def _ec_model(light=False): + return _build(light=light) + + +def _kcat(ec_model, rxn_id): + return float(ec_model.ec.kcat[ec_model.ec.rxns.index(rxn_id)]) + + +def test_set_kcat_does_not_apply_the_constraints_by_default(): + ec_model = _ec_model() + set_kcat_for_reactions(ec_model, ["R3"], 10.0) + coefficients = {m.id: c for m, c in ec_model.reactions.get_by_id("R3").metabolites.items()} + assert "prot_P4" not in coefficients + assert _kcat(ec_model, "R3") == 10.0 + + +def test_set_kcat_applies_the_constraints_when_asked(): + ec_model = _ec_model() + set_kcat_for_reactions(ec_model, ["R3"], 10.0, apply=True) + coefficients = {m.id: c for m, c in ec_model.reactions.get_by_id("R3").metabolites.items()} + assert coefficients["prot_P4"] == pytest.approx(-(40000.0 / (10.0 * 3600.0))) + + +def test_set_kcat_matches_ids_without_regard_to_letter_case(): + ec_model = _ec_model() + assert set_kcat_for_reactions(ec_model, ["r3"], 7.0) == ["R3"] + assert _kcat(ec_model, "R3") == 7.0 + assert set_kcat_for_reactions(ec_model, ["r2_EXP_1"], 8.0) == ["R2_EXP_1"] + assert _kcat(ec_model, "R2_EXP_1") == 8.0 + + +def test_set_kcat_ignores_an_unknown_id(): + ec_model = _ec_model() + assert set_kcat_for_reactions(ec_model, ["nonsense"], 1.0) == [] + assert set_kcat_for_reactions(ec_model, ["R3", "nonsense"], 1.0) == ["R3"] + assert _kcat(ec_model, "R3") == 1.0 + + +def test_set_kcat_on_a_light_model_does_not_expand_a_bare_cobra_id(): + ec_model = _ec_model(light=True) + assert set_kcat_for_reactions(ec_model, ["R2"], 5.0) == [] + assert _kcat(ec_model, "001_R2") == 0.0 + assert set_kcat_for_reactions(ec_model, ["001_R2", "002_R2"], 5.0) == ["001_R2", "002_R2"] + assert _kcat(ec_model, "002_R2") == 5.0 diff --git a/tests/test_parity_fixes_f3.py b/tests/test_parity_fixes_f3.py new file mode 100644 index 0000000..bbc509d --- /dev/null +++ b/tests/test_parity_fixes_f3.py @@ -0,0 +1,236 @@ +"""Behaviour shared with GECKO MATLAB: EC cell separators, zero-bound capacity, +the target growth of flexibilize_enz_concs and isozyme competition in gecko-light.""" +import copy +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter, make_ec_model +from geckopy.ec_model.ec_data import EcData +from geckopy.ec_model.pipeline import apply_kcat_constraints +from geckopy.get_enzyme_data import find_ec_in_db +from geckopy.limit_proteins import flexibilize_enz_concs +from geckopy.utilities import enzyme_usage, report_enzyme_usage + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + + +# --------------------------------------------------------------------------- # +# find_ec_in_db: a database cell may separate its codes with ';', spaces or both +# --------------------------------------------------------------------------- # + +def _db(rows): + codes = [ec for _, ec, _ in rows] + mw = np.array([m for _, _, m in rows], dtype=float) + index: dict[str, list[int]] = {} + for i, (gene, _, _) in enumerate(rows): + index.setdefault(gene, []).append(i) + return codes, mw, index + + +def test_semicolon_space_cell_intersects_with_partner_gene(): + codes, mw, index = _db([("g1", "1.1.1.1; 2.2.2.2", 10.0), ("g2", "2.2.2.2", 10.0)]) + assert find_ec_in_db(["g1", "g2"], codes, mw, index) == "2.2.2.2" + + +def test_space_separated_cell_is_split(): + codes, mw, index = _db([("g1", "1.1.1.4 1.1.1.303", 10.0)]) + assert find_ec_in_db(["g1"], codes, mw, index) == "1.1.1.4;1.1.1.303" + + +def test_semicolon_without_space_is_split(): + codes, mw, index = _db([("g1", "1.1.1.4;1.1.1.303", 10.0)]) + assert find_ec_in_db(["g1"], codes, mw, index) == "1.1.1.4;1.1.1.303" + + +def test_space_separated_wildcard_cell_drops_the_parent_wildcard(): + codes, mw, index = _db([("g1", "1.1.1.4 1.1.1.- 1.1.1.303", 10.0)]) + assert find_ec_in_db(["g1"], codes, mw, index) == "1.1.1.4;1.1.1.303" + + +# --------------------------------------------------------------------------- # +# enzyme_usage: a zero usage upper bound gives NaN, never 0 +# --------------------------------------------------------------------------- # + +def _toy_usage(enzyme_ub: float) -> EcModel: + model = EcModel("toy") + a_e = cobra.Metabolite("A_e", compartment="e") + a_c = cobra.Metabolite("A_c", compartment="c") + b_c = cobra.Metabolite("B_c", compartment="c") + pool = cobra.Metabolite("prot_pool", compartment="c") + prot_e = cobra.Metabolite("prot_E", compartment="c") + model.add_metabolites([a_e, a_c, b_c, pool, prot_e]) + ex_a = cobra.Reaction("EX_A") + ex_a.add_metabolites({a_e: -1.0}) + ex_a.bounds = (-1000.0, 0.0) + tr_a = cobra.Reaction("TR_A") + tr_a.add_metabolites({a_e: -1.0, a_c: 1.0}) + tr_a.bounds = (0.0, 1000.0) + r_ab = cobra.Reaction("R_AB") + r_ab.add_metabolites({a_c: -1.0, b_c: 1.0, prot_e: -1.0}) + r_ab.bounds = (0.0, 1000.0) + sk_b = cobra.Reaction("SK_B") + sk_b.add_metabolites({b_c: -1.0}) + sk_b.bounds = (0.0, 1000.0) + pool_ex = cobra.Reaction("prot_pool_exchange") + pool_ex.add_metabolites({pool: 1.0}) + pool_ex.bounds = (0.0, 1000.0) + usage = cobra.Reaction("usage_prot_E") + usage.add_metabolites({pool: -1.0, prot_e: 1.0}) + usage.bounds = (0.0, enzyme_ub) + model.add_reactions([ex_a, tr_a, r_ab, sk_b, pool_ex, usage]) + model.objective = "SK_B" + model.ec = EcData( + rxns=["R_AB"], kcat=np.array([1.0]), source=[""], notes=[""], eccodes=[""], + genes=["g_E"], enzymes=["E"], mw=np.array([100.0]), sequence=[""], + concs=np.array([float(enzyme_ub)]), rxn_enz_mat=sparse.csr_matrix([[1.0]]), + ) + return model + + +def test_zero_upper_bound_gives_nan_capacity(): + model = _toy_usage(0.0) + result = enzyme_usage(model, model.optimize().fluxes) + assert np.isnan(result.cap_usage[0]) + assert result.abs_usage[0] == 0.0 + + +def test_positive_upper_bound_gives_the_ratio(): + model = _toy_usage(5.0) + result = enzyme_usage(model, model.optimize().fluxes) + assert result.cap_usage[0] == pytest.approx(1.0) + + +def test_report_lists_no_zero_bound_enzyme_as_high_capacity(): + model = _toy_usage(0.0) + usage = enzyme_usage(model, model.optimize().fluxes) + report = report_enzyme_usage(model, usage) + assert len(report.high_cap_usage) == 0 + + +# --------------------------------------------------------------------------- # +# flexibilize_enz_concs: when no enzyme limits growth any more, the growth reached +# is the target of the final pass +# --------------------------------------------------------------------------- # + +def _flex_toy(tmp_path: Path, *, uptake_cap: float) -> EcModel: + (tmp_path / "model_adapter.toml").write_text( + 'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = 5.0\nbio_rxn = "biomass"\n' + ) + adapter = ModelAdapter.from_folder(tmp_path) + model = EcModel("toy", adapter=adapter) + a_e = cobra.Metabolite("A_e", compartment="e") + a_c = cobra.Metabolite("A_c", compartment="c") + b_c = cobra.Metabolite("B_c", compartment="c") + pool = cobra.Metabolite("prot_pool", compartment="c") + prot_e = cobra.Metabolite("prot_E", compartment="c") + model.add_metabolites([a_e, a_c, b_c, pool, prot_e]) + ex_a = cobra.Reaction("EX_A") + ex_a.add_metabolites({a_e: -1.0}) + ex_a.bounds = (-1000.0, 0.0) + tr_a = cobra.Reaction("TR_A") + tr_a.add_metabolites({a_e: -1.0, a_c: 1.0}) + tr_a.bounds = (0.0, uptake_cap) + r_ab = cobra.Reaction("R_AB") + r_ab.add_metabolites({a_c: -1.0, b_c: 1.0, prot_e: -1.0}) + r_ab.bounds = (0.0, 1000.0) + bio = cobra.Reaction("biomass") + bio.add_metabolites({b_c: -1.0}) + bio.bounds = (0.0, 1000.0) + pool_ex = cobra.Reaction("prot_pool_exchange") + pool_ex.add_metabolites({pool: 1.0}) + pool_ex.bounds = (0.0, 1000.0) + usage = cobra.Reaction("usage_prot_E") + usage.add_metabolites({pool: -1.0, prot_e: 1.0}) + usage.bounds = (0.0, 1.0) + model.add_reactions([ex_a, tr_a, r_ab, bio, pool_ex, usage]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R_AB"], kcat=np.array([1.0]), source=[""], notes=[""], eccodes=[""], + genes=["g_E"], enzymes=["E"], mw=np.array([100.0]), sequence=[""], + concs=np.array([1.0]), rxn_enz_mat=sparse.csr_matrix([[1.0]]), + ) + return model + + +def test_growth_capped_below_target_is_the_final_target(tmp_path): + model = _flex_toy(tmp_path, uptake_cap=2.5) + result = flexibilize_enz_concs(model) + assert result.uniprot_ids == ["E"] + upper = model.reactions.get_by_id("usage_prot_E").upper_bound + # The band's low edge below the reached growth 2.5, not an infeasible 5.0. + assert np.isfinite(upper) + assert upper == pytest.approx(2.5 * (1 - 0.25 / 100), rel=1e-3) + assert result.flex_concs[0] == pytest.approx(upper) + + +# --------------------------------------------------------------------------- # +# apply_kcat_constraints, gecko-light: naming one isozyme row updates the whole +# reaction, so every isozyme of it competes for the cheapest cost +# --------------------------------------------------------------------------- # + +_LIGHT_CACHE: EcModel | None = None + + +def _light_model() -> EcModel: + global _LIGHT_CACHE + if _LIGHT_CACHE is None: + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + gem = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + _LIGHT_CACHE = make_ec_model(gem, adapter, gecko_light=True) + return copy.deepcopy(_LIGHT_CACHE) + + +def _pool_coef(model: EcModel, rxn_id: str) -> float: + rxn = model.reactions.get_by_id(rxn_id) + for met, coef in rxn.metabolites.items(): + if met.id == "prot_pool": + return float(coef) + return 0.0 + + +def _row_cost(model: EcModel, row: int) -> float: + csr = model.ec.rxn_enz_mat.tocsr().getrow(row) + return float(np.sum(csr.data * model.ec.mw[csr.indices])) / float(model.ec.kcat[row]) + + +def test_naming_the_worse_isozyme_keeps_the_cheapest_cost_of_the_reaction(): + model = _light_model() + complex_row = model.ec.rxns.index("001_R2") + single_row = model.ec.rxns.index("002_R2") + model.ec.kcat[complex_row] = 50.0 + model.ec.kcat[single_row] = 50.0 + apply_kcat_constraints(model) + # The single-subunit isozyme is the cheaper one; make it worse and name only it. + model.ec.kcat[single_row] = 0.5 + cheapest = min(_row_cost(model, complex_row), _row_cost(model, single_row)) + assert _row_cost(model, complex_row) < _row_cost(model, single_row) + + apply_kcat_constraints(model, update_rxns=["002_R2"]) + + assert _pool_coef(model, "R2") == pytest.approx(-cheapest / 3600.0, rel=1e-12) + + +def test_naming_one_isozyme_leaves_other_reactions_untouched(): + model = _light_model() + for rxn_id in ("001_R2", "002_R2", "001_R3"): + model.ec.kcat[model.ec.rxns.index(rxn_id)] = 10.0 + apply_kcat_constraints(model, update_rxns=["001_R3"]) + assert _pool_coef(model, "R2") == 0.0 + assert _pool_coef(model, "R3") < 0.0 + + +def test_isozyme_row_without_enzymes_is_skipped(): + model = _light_model() + complex_row = model.ec.rxns.index("001_R2") + single_row = model.ec.rxns.index("002_R2") + model.ec.kcat[complex_row] = 50.0 + model.ec.kcat[single_row] = 50.0 + matrix = model.ec.rxn_enz_mat.tolil() + matrix[single_row, :] = 0.0 + model.ec.rxn_enz_mat = matrix.tocsr() + apply_kcat_constraints(model, update_rxns=["001_R2"]) + assert _pool_coef(model, "R2") == pytest.approx(-_row_cost(model, complex_row) / 3600.0, rel=1e-12) diff --git a/tests/test_parity_fixes_g1a.py b/tests/test_parity_fixes_g1a.py new file mode 100644 index 0000000..4a1a5b6 --- /dev/null +++ b/tests/test_parity_fixes_g1a.py @@ -0,0 +1,92 @@ +"""Behaviour that follows GECKO: truncate_values ties, the detail text of a kcat source file, +and the defaults of fill_kcats_from_isozymes.""" +import numpy as np +import pandas as pd +import pytest + +from geckopy.ec_model.pipeline import fill_kcats_from_isozymes, set_kcat_for_reactions +from geckopy.gather_kcats import read_kcat_source, save_kcat_source +from geckopy.kcat_tuning import truncate_values +from geckopy.kcat_tuning.truncate_values import _truncate_scalar + +from test_fill_kcats import _ectestgem_ec_model, _isozyme_model, _kcat +from test_kcat_sources import _adapter, _kcat_list, _model + + +# --------------------------------------------------------------------------- # +# truncate_values: a tie rounds away from zero, as MATLAB round(x, n) does +# --------------------------------------------------------------------------- # + +@pytest.mark.parametrize("value, expected", [ + (100000.5, 100001.0), + (1234565.0, 1234570.0), + (-100000.5, -100001.0), + (0.1234565, 0.123457), +]) +def test_truncate_tie_rounds_away_from_zero(value, expected): + assert _truncate_scalar(value) == expected + + +def test_truncate_array_ties_and_unchanged_values(): + out = truncate_values(np.array([100000.5, 1234565.0, 0.0, np.inf, np.nan, 12.3456789])) + assert out[0] == 100001.0 and out[1] == 1234570.0 + assert out[2] == 0.0 and np.isinf(out[3]) and np.isnan(out[4]) + assert out[5] == pytest.approx(12.3457) + + +# --------------------------------------------------------------------------- # +# save_kcat_source: the detail column can not split a row +# --------------------------------------------------------------------------- # + +def _saved_details(tmp_path, extras): + adapter = _adapter(tmp_path) + model = _model(adapter, light=False) + kcats = _kcat_list(model) + kcats["detail"] = [extras[i % len(extras)] for i in range(len(kcats))] + path = save_kcat_source(model, kcats, "brenda") + data_lines = [line for line in path.read_text(encoding="utf-8").split("\n") + if line and not line.startswith("#")] + return data_lines, read_kcat_source(path)[0] + + +def test_detail_with_tab_and_newline_keeps_one_row_per_kcat(tmp_path): + lines, records = _saved_details(tmp_path, ["a\tb", "c\nd"]) + assert {len(line.split("\t")) for line in lines} == {9} + assert len(records) == len(lines) - 1 + assert set(records["detail"]) <= {"wc=0; origin=%d; a b" % o for o in range(1, 7)} | { + "wc=0; origin=%d; c d" % o for o in range(1, 7)} + + +def test_detail_leaves_no_empty_parts(tmp_path): + _, records = _saved_details(tmp_path, ["; x", "p; ; q"]) + details = set(records["detail"]) + assert all("; ;" not in d and not d.endswith("; ") for d in details) + assert "wc=0; origin=1; x" in details + assert "wc=0; origin=2; p; q" in details + + +# --------------------------------------------------------------------------- # +# fill_kcats_from_isozymes: only a kcat of 0 is missing; apply is off by default +# --------------------------------------------------------------------------- # + +def test_nan_kcat_is_not_filled_from_siblings(): + model = _isozyme_model([5.0, float("nan"), 0.0]) + fill_kcats_from_isozymes(model) + assert np.isnan(model.ec.kcat[1]) + assert model.ec.source[1] == "x" + assert model.ec.kcat[2] == 5.0 and model.ec.source[2] == "isozymes" + + +def test_nan_sibling_does_not_enter_the_mean(): + model = _isozyme_model([4.0, float("nan"), 8.0, 0.0]) + fill_kcats_from_isozymes(model) + assert model.ec.kcat[3] == 6.0 + + +def test_apply_is_off_by_default(): + ec_model = _ectestgem_ec_model() + set_kcat_for_reactions(ec_model, ["R2_EXP_1"], 100.0, apply=False) + fill_kcats_from_isozymes(ec_model) + assert _kcat(ec_model, "R2_EXP_2") == 100.0 + r = ec_model.reactions.get_by_id("R2_EXP_2") + assert next((c for m, c in r.metabolites.items() if m.id == "prot_P3"), 0.0) == 0.0 diff --git a/tests/test_parity_fixes_g1b.py b/tests/test_parity_fixes_g1b.py new file mode 100644 index 0000000..355b64e --- /dev/null +++ b/tests/test_parity_fixes_g1b.py @@ -0,0 +1,239 @@ +"""fit_sigma and sensitivity_tuning follow GECKO's sigmaFitter and sensitivityTuning.""" +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter, fit_sigma, make_ec_model, sensitivity_tuning +from geckopy.ec_model.ec_data import EcData +from geckopy.ec_model.pipeline import apply_kcat_constraints, set_prot_pool_size + +from test_sigma_fitter import _adapter as _sigma_adapter +from test_sigma_fitter import _build_pool_only_model + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + + +# --------------------------------------------------------------------------- # +# fit_sigma +# --------------------------------------------------------------------------- # + +def test_default_scan_is_the_grid_of_sigma_fitter(tmp_path): + """growth = 1000 * sigma; a target of 10.5 lies between the grid points 0.01 and 0.02.""" + adapter = _sigma_adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=10.5) + model = _build_pool_only_model(adapter) + result = fit_sigma(model) + assert len(result.sigma_grid) == 100 + assert result.sigma == pytest.approx(0.01, abs=1e-12) + + +def test_bisect_is_still_selectable(tmp_path): + adapter = _sigma_adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=10.5) + model = _build_pool_only_model(adapter) + result = fit_sigma(model, method="bisect") + assert result.sigma == pytest.approx(0.0105, abs=1e-5) + assert len(result.sigma_grid) < 30 + + +def _capped_model(tmp_path, *, gr_exp, cap=30.0): + """Growth is min(1000 * sigma, cap): a second substrate, limited to ``cap``, + is consumed by the biomass reaction too.""" + adapter = _sigma_adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=gr_exp) + model = _build_pool_only_model(adapter) + substrate = cobra.Metabolite("cap_met", compartment="c") + supply = cobra.Reaction("EX_cap") + supply.add_metabolites({substrate: 1.0}) + supply.bounds = (0.0, cap) + model.add_reactions([supply]) + model.reactions.get_by_id("biomass").add_metabolites({substrate: -1.0}) + return model + + +def test_unreachable_target_gives_the_start_of_the_growth_plateau(tmp_path): + """Growth is min(1000 * sigma, 30); a target of 100 is out of reach, so the + first grid point with the plateau's growth wins, as in sigmaFitter.m.""" + model = _capped_model(tmp_path, gr_exp=100.0) + result = fit_sigma(model) + assert result.sigma == pytest.approx(0.03, abs=1e-12) + + +def test_bisect_finds_the_start_of_the_plateau_for_an_unreachable_target(tmp_path): + model = _capped_model(tmp_path, gr_exp=100.0) + result = fit_sigma(model, method="bisect") + assert result.sigma == pytest.approx(0.03, abs=1e-5) + + +def test_objective_bounds_are_opened_and_stay_open(tmp_path): + """sigmaFitter.m sets lb 0 and ub 1000 on the objective reaction for good.""" + adapter = _sigma_adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=10.0) + model = _build_pool_only_model(adapter) + biomass = model.reactions.get_by_id("biomass") + biomass.bounds = (5.0, 5.0) + result = fit_sigma(model) + assert biomass.bounds == (0.0, 1000.0) + assert result.sigma == pytest.approx(0.01, abs=1e-12) + + +# --------------------------------------------------------------------------- # +# sensitivity_tuning, full model +# --------------------------------------------------------------------------- # + +def _toy_adapter(tmp_path: Path) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text( + 'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = 0.5\nbio_rxn = "biomass"\n' + ) + return ModelAdapter.from_folder(tmp_path) + + +def _chain(adapter, kcats, mws=(3600.0, 3600.0), pool=1.0) -> EcModel: + """A -> B (R1, enzyme E1) -> biomass (R2, enzyme E2), both drawing on one pool. + + The pool cost of one unit of flux through Ri is mw_i / (kcat_i * 3600).""" + model = EcModel("chain", adapter=adapter) + mets = {name: cobra.Metabolite(name, compartment="c") for name in + ("A_e", "A_c", "B_c", "pool", "prot_E1", "prot_E2", "bio_met")} + mets["pool"].id = "prot_pool" + model.add_metabolites(list(mets.values())) + + def rxn(rid, stoich, lb, ub): + r = cobra.Reaction(rid) + r.add_metabolites({mets[k]: v for k, v in stoich.items()}) + r.bounds = (lb, ub) + return r + + c1 = mws[0] / (kcats[0] * 3600.0) + c2 = mws[1] / (kcats[1] * 3600.0) + model.add_reactions([ + rxn("EX_A", {"A_e": -1.0}, -1000.0, 0.0), + rxn("TR_A", {"A_e": -1.0, "A_c": 1.0}, 0.0, 1000.0), + rxn("R1", {"A_c": -1.0, "prot_E1": -c1, "B_c": 1.0}, 0.0, 1000.0), + rxn("R2", {"B_c": -1.0, "prot_E2": -c2, "bio_met": 1.0}, 0.0, 1000.0), + rxn("biomass", {"bio_met": -1.0}, 0.0, 1000.0), + rxn("prot_pool_exchange", {"pool": 1.0}, 0.0, pool), + rxn("usage_prot_E1", {"pool": -1.0, "prot_E1": 1.0}, 0.0, 1000.0), + rxn("usage_prot_E2", {"pool": -1.0, "prot_E2": 1.0}, 0.0, 1000.0), + ]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R1", "R2"], + kcat=np.array(kcats, dtype=float), + source=["initial", "initial"], + notes=["", ""], + eccodes=["", ""], + genes=["g1", "g2"], + enzymes=["E1", "E2"], + mw=np.array(mws, dtype=float), + sequence=["", ""], + concs=np.array([np.nan, np.nan]), + rxn_enz_mat=sparse.csr_matrix(np.eye(2)), + ) + return model + + +def test_tuned_reactions_are_listed_in_model_order(tmp_path): + """R2 costs more pool, so it is tuned first; the result lists R1 before R2.""" + model = _chain(_toy_adapter(tmp_path), kcats=[1.0, 0.5]) + result = sensitivity_tuning(model, desired_growth_rate=1.0, verbose=False) + assert result.rxns == ["R1", "R2"] + np.testing.assert_allclose(result.old_kcat, [1.0, 0.5]) + np.testing.assert_allclose(result.new_kcat, [10.0, 5.0]) + + +def test_rerun_reports_the_kcat_and_source_of_the_input_model(tmp_path): + model = _chain(_toy_adapter(tmp_path), kcats=[1.0, 0.5]) + sensitivity_tuning(model, desired_growth_rate=1.0, verbose=False) + result = sensitivity_tuning(model, desired_growth_rate=30.0, verbose=False) + assert result.rxns + assert not np.isnan(result.old_kcat).any() + assert all(src == "sensitivityTuning" for src in result.source) + # The values the second run started from are the tuned ones of the first. + assert set(result.old_kcat) <= {10.0, 5.0} + + +@pytest.mark.parametrize("value, text", [ + (10.0, "10"), + (float(np.pi), "3.1416"), + (0.5, "0.5"), + (123456.789, "123456.789"), + (1e7, "10000000"), + (float("nan"), "NaN"), +]) +def test_num2str_prints_a_scalar_as_matlab_does(value, text): + from geckopy.kcat_tuning.sensitivity_tuning import _num2str + assert _num2str(value) == text + + +def test_notes_hold_the_pre_tune_kcat_as_num2str_prints_it(tmp_path): + model = _chain(_toy_adapter(tmp_path), kcats=[10.0, float(np.pi)]) + sensitivity_tuning(model, desired_growth_rate=30.0, verbose=False) + assert "preTuneKcat=10 |" in model.ec.notes[0] + assert "preTuneKcat=3.1416 |" in model.ec.notes[1] + + +def test_a_tuned_kcat_has_no_upper_limit(tmp_path): + """kcat 2e6 and mw 7.2e9 give growth 1; a target of 100 needs 2e8, past the old 1e7 cap.""" + model = _chain(_toy_adapter(tmp_path), kcats=[2e6, 1e12], mws=(7.2e9, 3600.0)) + result = sensitivity_tuning(model, desired_growth_rate=99.0, verbose=False) + assert result.rxns == ["R1"] + assert result.new_kcat[0] == pytest.approx(2e8) + + +def test_prot_to_ignore_matches_without_regard_to_case(tmp_path): + model = _chain(_toy_adapter(tmp_path), kcats=[1.0, 0.5]) + result = sensitivity_tuning( + model, desired_growth_rate=1.0, prot_to_ignore=["e2"], verbose=False, + ) + # E2 is ignored, so only R1 is ever tuned. + assert result.rxns == ["R1"] + + +def test_prot_to_ignore_with_an_unknown_id_raises(tmp_path): + model = _chain(_toy_adapter(tmp_path), kcats=[1.0, 0.5]) + with pytest.raises(ValueError, match="prot_to_ignore"): + sensitivity_tuning( + model, desired_growth_rate=1.0, prot_to_ignore=["NOPE"], verbose=False, + ) + + +# --------------------------------------------------------------------------- # +# sensitivity_tuning, gecko-light model +# --------------------------------------------------------------------------- # + +def _light_model(kcats): + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + conv = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + for rxn_id in ("R2", "R4"): + conv.reactions.get_by_id(rxn_id).bounds = (0.0, 0.0) + model = make_ec_model(conv, adapter, gecko_light=True) + model.adapter = adapter + for i, ec_id in enumerate(model.ec.rxns): + base = str(ec_id).split("_", 1)[1] + if base in kcats: + model.ec.kcat[i] = float(kcats[base]) + model.objective = "R5" + apply_kcat_constraints(model) + set_prot_pool_size(model) + return model + + +def test_gecko_light_model_is_tuned(): + model = _light_model({"R3": 0.05, "R5": 200.0}) + before = float(model.optimize().objective_value) + kcat_before = np.array(model.ec.kcat, dtype=float) + source_before = list(model.ec.source) + result = sensitivity_tuning( + model, desired_growth_rate=30.0, bio_rxn="R5", verbose=False, + ) + assert result.rxns + assert float(model.optimize().objective_value) > before + # ec.source is left as it was on a light model, as in sensitivityTuning.m. + assert list(model.ec.source) == source_before + # Every isozyme row of a tuned reaction was multiplied by the same power of 10. + for rxn_id in result.rxns: + rows = [i for i, ec_id in enumerate(model.ec.rxns) + if str(ec_id).split("_", 1)[1] == rxn_id and kcat_before[i] > 0] + ratios = [model.ec.kcat[i] / kcat_before[i] for i in rows] + assert ratios and np.allclose(ratios, ratios[0]) + assert ratios[0] >= 10.0 - 1e-9 diff --git a/tests/test_parity_round3_u2.py b/tests/test_parity_round3_u2.py new file mode 100644 index 0000000..8de6203 --- /dev/null +++ b/tests/test_parity_round3_u2.py @@ -0,0 +1,140 @@ +"""An empty update list in apply_kcat_constraints, and a quoted reaction list in +apply_custom_kcats, as in GECKO.""" +import copy +from pathlib import Path +from types import SimpleNamespace + +import cobra +import numpy as np +import pytest + +from geckopy import EcModel, ModelAdapter, make_ec_model +from geckopy.ec_model.pipeline import ( + apply_custom_kcats, + apply_kcat_constraints, + set_kcat_for_reactions, +) + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" +HEADER = "proteins\tgenes\tgene_name\tkcat\trxns\tnotes\tstoicho\n" + +_CACHE: dict[bool, EcModel] = {} + + +def _model(light: bool = False) -> EcModel: + if light not in _CACHE: + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + cobra_model = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + _CACHE[light] = make_ec_model(cobra_model, adapter, gecko_light=light) + return copy.deepcopy(_CACHE[light]) + + +def _coef(model: EcModel, rxn_id: str, met_id: str) -> float: + for met, coef in model.reactions.get_by_id(rxn_id).metabolites.items(): + if met.id == met_id: + return coef + return 0.0 + + +@pytest.mark.parametrize("empty", [[], (), np.array([], dtype=str)]) +def test_empty_update_list_updates_all_reactions(empty): + model = _model() + model.ec.kcat[model.ec.rxns.index("R3")] = 10.0 + assert _coef(model, "R3", "prot_P4") == 0.0 + + apply_kcat_constraints(model, update_rxns=empty) + + assert _coef(model, "R3", "prot_P4") < 0.0 + + +def test_empty_update_list_updates_all_reactions_light(): + model = _model(light=True) + model.ec.kcat[model.ec.rxns.index("001_R3")] = 10.0 + assert _coef(model, "R3", "prot_pool") == 0.0 + + apply_kcat_constraints(model, update_rxns=[]) + + assert _coef(model, "R3", "prot_pool") < 0.0 + + +def test_a_list_naming_one_reaction_still_updates_only_that_one(): + model = _model() + model.ec.kcat[model.ec.rxns.index("R3")] = 10.0 + model.ec.kcat[model.ec.rxns.index("R5")] = 10.0 + + apply_kcat_constraints(model, update_rxns=["R3"]) + + assert _coef(model, "R3", "prot_P4") < 0.0 + assert _coef(model, "R5", "prot_P5") == 0.0 + + +def test_set_kcat_for_reactions_with_an_empty_id_list_changes_nothing(): + model = _model() + model.ec.kcat[model.ec.rxns.index("R3")] = 10.0 + + updated = set_kcat_for_reactions(model, [], 5.0, apply=True) + + assert updated == [] + assert _coef(model, "R3", "prot_P4") == 0.0 + + +def test_set_kcat_for_reactions_with_only_unknown_ids_does_not_rewrite_the_model(): + # an unknown id changes no row; apply=True must not hand an empty list to + # apply_kcat_constraints, where an empty list now means every reaction + model = _model() + model.ec.kcat[model.ec.rxns.index("R3")] = 10.0 + + updated = set_kcat_for_reactions(model, ["nonsense"], 5.0, apply=True) + + assert updated == [] + assert _coef(model, "R3", "prot_P4") == 0.0 + + +def test_evotune_scoring_skips_the_constraint_update_for_no_tunable_rows(monkeypatch): + from geckopy.kcat_tuning.evotune import tuning + + calls = [] + + class _Stop(Exception): + pass + + def _record(model, update_rxns=None): + calls.append(update_rxns) + + def _stop(model): + raise _Stop + + monkeypatch.setattr(tuning, "apply_kcat_constraints", _record) + monkeypatch.setattr(tuning, "_reset_solver_basis", _stop) + model = SimpleNamespace(ec=SimpleNamespace(kcat=np.zeros(3))) + + with pytest.raises(_Stop): + tuning._score_kcat_vector( + model, np.array([], dtype=int), [], None, {}, "bio", + np.array([]), make_anaerobic=None, change_protein_biomass=None, + ) + + assert calls == [] + + +def test_quoted_reaction_list_matches_the_listed_reactions(tmp_path): + model = _model() + tsv = tmp_path / "customKcats.tsv" + tsv.write_text(HEADER + '\t\t\t50\t"R3,R5"\tquoted\t1\n') + + apply_custom_kcats(model, path=tsv, apply=False) + + for rxn in ("R3", "R5"): + assert model.ec.kcat[model.ec.rxns.index(rxn)] == 50.0 + assert model.ec.source[model.ec.rxns.index(rxn)] == "custom" + + +def test_unquoted_and_single_reaction_lists_are_unchanged(tmp_path): + model = _model() + tsv = tmp_path / "customKcats.tsv" + tsv.write_text(HEADER + "\t\t\t50\tR3,R5\tplain\t1\n") + + apply_custom_kcats(model, path=tsv, apply=False) + + assert model.ec.kcat[model.ec.rxns.index("R3")] == 50.0 + assert model.ec.kcat[model.ec.rxns.index("R5")] == 50.0 diff --git a/tests/test_review_fixes_r1p.py b/tests/test_review_fixes_r1p.py new file mode 100644 index 0000000..8d52768 --- /dev/null +++ b/tests/test_review_fixes_r1p.py @@ -0,0 +1,49 @@ +"""Edge cases of get_conc_control_coeffs, fit_sigma and the docstring of fill_kcats_from_isozymes.""" +import numpy as np +import pytest + +from geckopy.kcat_tuning import fit_sigma +from geckopy.limit_proteins import get_conc_control_coeffs +from geckopy.limit_proteins.get_conc_control_coeffs import _finite_difference_coeffs + +from test_get_conc_control_coeffs import _build_toy_ec_model +from test_sigma_fitter import _adapter, _build_pool_only_model + +SOLVERS = ["glpk", "glpk_exact", "gurobi"] + + +def _with_solver(model, solver): + try: + model.solver = solver + model.optimize() + except Exception as exc: # solver not installed or not licensed + pytest.skip(f"solver {solver} unavailable: {exc}") + return model + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_small_coefficient_on_a_small_bound_is_kept(solver): + # One unit of enzyme per 2000 units of flux: growth = ub / 2000, so the coefficient is 5e-4 + # while a 0.01 percent rise of a bound of 1e-3 raises growth by only 5e-11. + model = _with_solver(_build_toy_ec_model(enzyme_ub=1e-3, enzyme_coeff_in_R=-2000.0), solver) + assert model.optimize().objective_value == pytest.approx(5e-7, rel=1e-6) + enz, coeffs = get_conc_control_coeffs(model) + _, local = _finite_difference_coeffs(model, ["E"], 2.0, 0.0) + assert enz[0] + assert local[0] == pytest.approx(5e-4, rel=1e-6) + assert coeffs[0] == pytest.approx(5e-4, rel=1e-6) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_large_coefficient_is_unchanged(solver): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) + + +@pytest.mark.parametrize("method", ["grid", "bisect"]) +def test_a_model_that_cannot_grow_gives_the_first_grid_sigma(tmp_path, method): + adapter = _adapter(tmp_path, p_tot=0.0) + model = _build_pool_only_model(adapter) + result = fit_sigma(model, growth_rate=0.5, p_tot=0.0, f=1.0, method=method) + assert result.sigma == pytest.approx(0.01) diff --git a/tests/test_review_fixes_r2p.py b/tests/test_review_fixes_r2p.py new file mode 100644 index 0000000..906d9ea --- /dev/null +++ b/tests/test_review_fixes_r2p.py @@ -0,0 +1,126 @@ +"""get_conc_control_coeffs: the verification rule of the reduced-cost path. + +A reduced cost is kept when raising the usage bound by 0.01 percent raises the optimum by at +least half of the rise the reduced cost predicts and by at least 1e-12 times max(1, |optimum|). +Every nonzero reduced cost is checked, however small. A usage reaction held at its lower bound +is not checked. The check solve is replaced by a stub that reports a chosen rise, so each case +fixes the realised gain exactly. +""" +import sys +from types import SimpleNamespace + +import pytest + +from geckopy.limit_proteins import get_conc_control_coeffs + +from test_get_conc_control_coeffs import _build_toy_ec_model + +SOLVERS = ["glpk", "glpk_exact", "gurobi"] +MODULE = "geckopy.limit_proteins.get_conc_control_coeffs" + + +def _with_solver(model, solver): + try: + model.solver = solver + model.optimize() + except Exception as exc: # solver not installed or not licensed + pytest.skip(f"solver {solver} unavailable: {exc}") + return model + + +def _stub_check_solve(monkeypatch, model, gain): + """The first solve is real; later solves report the first optimum plus ``gain``.""" + real = model.optimize + state = {"calls": 0, "initial": None} + + def optimize(*args, **kwargs): + sol = real(*args, **kwargs) + state["calls"] += 1 + if state["calls"] == 1: + state["initial"] = sol.objective_value + return sol + return SimpleNamespace( + status="optimal", objective_value=state["initial"] + gain, fluxes=sol.fluxes + ) + + monkeypatch.setattr(model, "optimize", optimize, raising=False) + return state + + +def _fixed_reduced_cost(monkeypatch, value): + module = sys.modules[MODULE] + monkeypatch.setattr(module, "_forward_reduced_cost", lambda rxn, flux=None: value) + + +# The toy with one enzyme of bound 5 and one unit of enzyme per unit of flux has reduced cost 1, +# so the predicted gain of the 0.01 percent step is 1 * 5 * 1e-4 = 5e-4. + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_gain_below_half_the_prediction_gives_zero(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _stub_check_solve(monkeypatch, model, 0.3 * 5e-4) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == 0.0 + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_gain_above_half_the_prediction_keeps_the_reduced_cost(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _stub_check_solve(monkeypatch, model, 0.6 * 5e-4) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_gain_below_the_floor_is_noise(solver, monkeypatch): + # Reduced cost 1e-8 on a bound of 5 predicts 5e-12; the floor is 1e-12 * max(1, 5) = 5e-12. + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _fixed_reduced_cost(monkeypatch, 1e-8) + _stub_check_solve(monkeypatch, model, 4e-12) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == 0.0 + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_gain_above_the_floor_is_kept(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _fixed_reduced_cost(monkeypatch, 1e-8) + _stub_check_solve(monkeypatch, model, 6e-12) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == pytest.approx(1e-8) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_reduced_cost_below_1e_10_is_checked_and_kept(solver, monkeypatch): + # Reduced cost 5e-11 on a bound of 500 predicts 2.5e-12, above the floor of 1e-12. + model = _with_solver(_build_toy_ec_model(enzyme_ub=500.0, enzyme_coeff_in_R=-2000.0), solver) + _fixed_reduced_cost(monkeypatch, 5e-11) + _stub_check_solve(monkeypatch, model, 2.5e-12) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == pytest.approx(5e-11) + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_zero_reduced_cost_is_not_checked(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + _fixed_reduced_cost(monkeypatch, 0.0) + state = _stub_check_solve(monkeypatch, model, 1.0) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == 0.0 + assert state["calls"] == 1 + + +@pytest.mark.parametrize("solver", SOLVERS) +def test_a_usage_reaction_held_at_its_lower_bound_is_not_checked(solver, monkeypatch): + model = _with_solver(_build_toy_ec_model(enzyme_ub=5.0), solver) + usage = model.reactions.get_by_id("usage_prot_E") + usage.lower_bound = 2.0 + usage.upper_bound = 2.0 + _fixed_reduced_cost(monkeypatch, 1.0) + state = _stub_check_solve(monkeypatch, model, 1.0) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == 0.0 + assert state["calls"] == 1 diff --git a/tests/test_review_fixes_r3p.py b/tests/test_review_fixes_r3p.py new file mode 100644 index 0000000..08a540e --- /dev/null +++ b/tests/test_review_fixes_r3p.py @@ -0,0 +1,83 @@ +"""Reduced costs without usable duals in flexibilize_enz_concs, and the positional kcat rule of +set_kcat_for_reactions when an input id matches nothing.""" +import copy +import importlib +from pathlib import Path + +import cobra +import numpy as np +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.ec_model.pipeline import set_kcat_for_reactions +from geckopy.limit_proteins.flexibilize_enz_concs import _unit_scores + +from test_conc_control_coeffs_shadow import _with_solver +from test_flexibilize_enz_concs import _adapter, _build_toy + +gcc = importlib.import_module("geckopy.limit_proteins.get_conc_control_coeffs") + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" +_CACHE = {} + + +def _ectestgem(): + if "m" not in _CACHE: + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + cobra_model = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + _CACHE["m"] = make_ec_model(cobra_model, adapter) + return copy.deepcopy(_CACHE["m"]) + + +def _scipy_toy(tmp_path): + model = _with_solver(_build_toy(_adapter(tmp_path, gr_exp=5.0, bio_rxn="biomass"), enzyme_ub=1.0, measured_conc=1.0), "scipy") + reduced = model.optimize().reduced_costs + if not np.isnan(np.asarray(reduced, dtype=float)).all(): + pytest.skip("this cobra returns usable reduced costs with the scipy solver") + return model + + +def _optimal_in_sequence(monkeypatch, answers): + calls = iter(answers) + real = gcc._solution_is_optimal + monkeypatch.setattr(gcc, "_solution_is_optimal", lambda sol: next(calls, real(sol))) + + +def test_unit_scores_are_zero_not_nan_when_the_lifted_pool_solve_is_not_optimal(tmp_path, monkeypatch): + model = _scipy_toy(tmp_path) + _optimal_in_sequence(monkeypatch, [True, False]) + scores = _unit_scores(model, ["E"], [[0]], 1.1, 0.75) + assert not np.isnan(scores).any() + assert (scores == 0.0).all() + + +def test_unit_scores_are_zero_not_nan_when_a_probe_is_not_optimal(tmp_path, monkeypatch): + model = _scipy_toy(tmp_path) + _optimal_in_sequence(monkeypatch, [True, True, False]) + scores = _unit_scores(model, ["E"], [[0]], 1.1, 0.75) + assert not np.isnan(scores).any() + assert (scores == 0.0).all() + + +def test_unit_scores_with_nan_reduced_costs_fall_back_to_the_probes(tmp_path): + model = _scipy_toy(tmp_path) + scores = _unit_scores(model, ["E"], [[0]], 1.1, 0.75) + assert not np.isnan(scores).any() + assert scores[0] > 0.0 + + +def test_a_known_id_next_to_an_unknown_id_takes_the_values_by_position(): + model = _ectestgem() + rows = [i for i, r in enumerate(model.ec.rxns) if r.split("_EXP_")[0] == "R2"] + assert len(rows) == 2 + updated = set_kcat_for_reactions(model, ["R2", "nonsense"], [1.0, 2.0]) + assert sorted(updated) == sorted(model.ec.rxns[i] for i in rows) + assert [float(model.ec.kcat[i]) for i in rows] == [1.0, 2.0] + + +def test_one_value_per_input_id_still_broadcasts_when_every_id_matches(): + model = _ectestgem() + set_kcat_for_reactions(model, ["R2", "R3"], [7.0, 9.0]) + r2 = [float(model.ec.kcat[i]) for i, r in enumerate(model.ec.rxns) if r.split("_EXP_")[0] == "R2"] + assert r2 == [7.0, 7.0] + assert float(model.ec.kcat[model.ec.rxns.index("R3")]) == 9.0 diff --git a/tests/test_scipy_duals_t3.py b/tests/test_scipy_duals_t3.py new file mode 100644 index 0000000..463fe16 --- /dev/null +++ b/tests/test_scipy_duals_t3.py @@ -0,0 +1,141 @@ +"""get_conc_control_coeffs with optlang's scipy backend: the reduced cost comes from HiGHS's bound marginal. + +optlang raises on ``Variable.dual`` for the scipy solver but keeps the ``linprog`` result, whose +``upper.marginals`` hold d(objective)/d(upper bound) of the minimised problem. The function reads +that value, so scipy takes the same one-solve path as the other solvers. A result without usable +marginals, or one that does not belong to the current variables, falls back to the +finite-difference loop. +""" +import importlib +import logging +from types import SimpleNamespace + +import numpy as np +import pytest + +from geckopy.limit_proteins import get_conc_control_coeffs +from geckopy.limit_proteins.get_conc_control_coeffs import ( + _finite_difference_coeffs, + _forward_reduced_cost, + _scipy_upper_marginal, +) + +from test_conc_control_coeffs_shadow import _degenerate_model, _pool_binding_model, _with_solver +from test_get_conc_control_coeffs import _build_toy_ec_model + +mod = importlib.import_module("geckopy.limit_proteins.get_conc_control_coeffs") + + +def _forbid_finite_difference(monkeypatch): + def fail(*args, **kwargs): + raise AssertionError("the finite-difference loop was used") + + monkeypatch.setattr(mod, "_finite_difference_coeffs", fail) + + +def _scipy_model(**kwargs): + return _with_solver(_build_toy_ec_model(**kwargs), "scipy") + + +def _flux(model, rxn): + """The flux cobra reports for ``rxn`` in a fresh solve of the model.""" + return float(model.optimize().fluxes[rxn.id]) + + +def test_scipy_reads_the_bound_marginal_and_not_the_finite_difference(monkeypatch): + model = _scipy_model(enzyme_ub=5.0) + _forbid_finite_difference(monkeypatch) + enz, coeffs = get_conc_control_coeffs(model) + assert enz[0] + assert coeffs[0] == pytest.approx(1.0, rel=1e-6) + + +def test_scipy_marginal_of_a_binding_bound_equals_the_marginal_growth(): + model = _scipy_model(enzyme_ub=5.0, enzyme_coeff_in_R=-2.0) + rxn = model.reactions.get_by_id("usage_prot_E") + assert _forward_reduced_cost(rxn, _flux(model, rxn)) == pytest.approx(0.5, rel=1e-6) + _, local = _finite_difference_coeffs(model, ["E"], 1.001, 0.0) + assert local[0] == pytest.approx(0.5, rel=1e-6) + + +def test_scipy_gives_zero_for_a_bound_that_does_not_bind_while_the_pool_binds(monkeypatch): + model = _pool_binding_model("scipy") + rxn = model.reactions.get_by_id("usage_prot_E") + _forbid_finite_difference(monkeypatch) + assert _forward_reduced_cost(rxn, _flux(model, rxn)) == 0.0 + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == 0.0 + + +def test_scipy_degenerate_vertex_is_resolved_by_the_verification_step(monkeypatch): + model = _degenerate_model("scipy") + rxn = model.reactions.get_by_id("usage_prot_E") + # HiGHS may put the dual of this vertex on R_AB's bound, which leaves nothing to verify: force + # the nonzero dual another vertex gives, so the verification step is what removes it + assert _forward_reduced_cost(rxn, _flux(model, rxn)) >= 0.0 + monkeypatch.setattr(mod, "_forward_reduced_cost", lambda rxn, flux=None: 1.0) + _forbid_finite_difference(monkeypatch) + _, coeffs = get_conc_control_coeffs(model) + assert coeffs[0] == 0.0 + + +@pytest.mark.parametrize("other", ["glpk", "gurobi"]) +@pytest.mark.parametrize( + "build", + [ + lambda s: _with_solver(_build_toy_ec_model(enzyme_ub=5.0), s), + lambda s: _with_solver(_build_toy_ec_model(enzyme_ub=5.0, enzyme_coeff_in_R=-2.0), s), + _pool_binding_model, + _degenerate_model, + ], + ids=["binding", "coefficient_two", "pool_binding", "degenerate"], +) +def test_scipy_agrees_with_the_other_solvers(build, other): + reference = build(other) + scipy_model = build("scipy") + enz_ref, ref = get_conc_control_coeffs(reference) + enz_sp, sp = get_conc_control_coeffs(scipy_model) + np.testing.assert_array_equal(enz_ref, enz_sp) + np.testing.assert_allclose(sp, ref, rtol=1e-6, atol=1e-12) + + +def test_scipy_result_without_bound_marginals_falls_back(caplog, monkeypatch): + model = _scipy_model(enzyme_ub=8.0) + model.reactions.get_by_id("EX_A").lower_bound = -10.0 + problem_class = type(model.solver.problem) + original = problem_class.optimize + + def optimize_without_marginals(self, *args, **kwargs): + original(self, *args, **kwargs) + self._solution = SimpleNamespace(x=self._solution.x, status=self._solution.status) + + monkeypatch.setattr(problem_class, "optimize", optimize_without_marginals) + mod._FALLBACK_REPORTED.clear() + with caplog.at_level(logging.INFO, logger=mod.__name__): + pm, pc = get_conc_control_coeffs(model, fold_change=2.0) + fm, fc = _finite_difference_coeffs(model, ["E"], 2.0, 0.0) + np.testing.assert_array_equal(pm, fm) + np.testing.assert_allclose(pc, fc, rtol=1e-6, atol=1e-12) + assert "returns no usable variable duals" in caplog.text + + +def test_scipy_result_that_does_not_belong_to_the_current_variables_gives_none(): + model = _scipy_model(enzyme_ub=5.0) + rxn = model.reactions.get_by_id("usage_prot_E") + flux = _flux(model, rxn) + problem = model.solver.problem + assert _scipy_upper_marginal(rxn, flux) is not None + assert _scipy_upper_marginal(rxn, flux + 1.0) is None # not the solution cobra returned + assert _scipy_upper_marginal(rxn, None) is None + problem._solution.x[:] = 0.0 # a stored result that is not the current solution + assert _scipy_upper_marginal(rxn, flux) is None + + +def test_nan_marginal_falls_back(monkeypatch): + model = _scipy_model(enzyme_ub=5.0) + monkeypatch.setattr(mod, "_scipy_upper_marginal", lambda rxn, flux: float("nan")) + rxn = model.reactions.get_by_id("usage_prot_E") + assert _forward_reduced_cost(rxn, 5.0) is None + enz, coeffs = get_conc_control_coeffs(model) + _, local = _finite_difference_coeffs(model, ["E"], 2.0, 0.0) + assert coeffs[0] == pytest.approx(local[0], rel=1e-6) diff --git a/tests/test_sensitivity_tuning.py b/tests/test_sensitivity_tuning.py index b28446f..294aa8f 100644 --- a/tests/test_sensitivity_tuning.py +++ b/tests/test_sensitivity_tuning.py @@ -108,14 +108,6 @@ def _build_toy( # Pre-flight # --------------------------------------------------------------------------- # -def test_gecko_light_raises(tmp_path): - adapter = _adapter(tmp_path) - model = _build_toy(adapter) - model.ec.gecko_light = True - with pytest.raises(NotImplementedError, match="gecko-light"): - sensitivity_tuning(model) - - def test_no_adapter_no_growth_rate_raises(tmp_path): adapter = _adapter(tmp_path) model = _build_toy(adapter) diff --git a/tests/test_set_kcat.py b/tests/test_set_kcat.py index 6413c71..7187af6 100644 --- a/tests/test_set_kcat.py +++ b/tests/test_set_kcat.py @@ -181,17 +181,31 @@ def test_kcat_list_wrong_length_raises(): # Unknown ID, empty input # --------------------------------------------------------------------------- # -def test_unknown_rxn_id_raises(): +def test_unknown_rxn_id_is_ignored(): ec_model = _ectestgem_ec_model() - with pytest.raises(ValueError, match="matched no entries"): - set_kcat_for_reactions(ec_model, ["nonsense"], 1.0) + before = list(ec_model.ec.kcat) + assert set_kcat_for_reactions(ec_model, ["nonsense"], 1.0) == [] + assert list(ec_model.ec.kcat) == before -def test_unknown_with_known_still_raises(): - """If any one ID is unknown, the whole call must fail.""" +def test_unknown_with_known_changes_the_known_one(): + """An unknown ID is ignored; the known one is still set.""" ec_model = _ectestgem_ec_model() - with pytest.raises(ValueError, match="matched no entries"): - set_kcat_for_reactions(ec_model, ["R3", "nonsense"], 1.0) + updated = set_kcat_for_reactions(ec_model, ["R3", "nonsense"], 1.0) + assert updated == ["R3"] + assert _kcat_at(ec_model, "R3") == 1.0 + + +def test_ids_match_without_regard_to_letter_case(): + """The ``_EXP_`` suffix itself is recognised in capitals only, as in + MATLAB; the rest of the ID is compared without regard to case.""" + ec_model = _ectestgem_ec_model() + updated = set_kcat_for_reactions(ec_model, ["r3"], 7.0) + assert updated == ["R3"] + assert _kcat_at(ec_model, "R3") == 7.0 + updated = set_kcat_for_reactions(ec_model, ["r2_EXP_1"], 8.0) + assert updated == ["R2_EXP_1"] + assert _kcat_at(ec_model, "R2_EXP_1") == 8.0 def test_empty_rxn_ids_returns_empty(): @@ -221,12 +235,14 @@ def test_apply_false_leaves_s_matrix_unchanged(): assert _get_s_coef(r3, "prot_P4") == 0.0 -def test_apply_default_is_true(): - """The default value of apply should be True.""" +def test_apply_default_is_false(): + """As MATLAB setKcatForReactions, only ec.kcat and ec.source change by + default; the S matrix follows once apply_kcat_constraints runs.""" ec_model = _ectestgem_ec_model() set_kcat_for_reactions(ec_model, ["R3"], 10.0) # no apply arg r3 = ec_model.reactions.get_by_id("R3") - assert _get_s_coef(r3, "prot_P4") != 0.0 + assert _get_s_coef(r3, "prot_P4") == 0.0 + assert _kcat_at(ec_model, "R3") == 10.0 def test_apply_only_updates_passed_reactions(): @@ -240,7 +256,7 @@ def test_apply_only_updates_passed_reactions(): assert r5_pre != 0.0 # Now use set_kcat_for_reactions to update R3. R5 must be untouched. - set_kcat_for_reactions(ec_model, ["R3"], 10.0) + set_kcat_for_reactions(ec_model, ["R3"], 10.0, apply=True) r5_post = _get_s_coef(ec_model.reactions.get_by_id("R5"), "prot_P5") assert r5_post == r5_pre @@ -287,24 +303,29 @@ def test_light_set_single_explicit_prefixed_id(): assert _source_at(ec, "001_R3") == "setKcatForReactions" -def test_light_unsuffixed_id_expands_to_all_isozyme_rows(): - """Bare ``R2`` should match both ``001_R2`` and ``002_R2`` (R2's two - isozymes), but not ``001_R2_REV`` / ``002_R2_REV``.""" +def test_light_base_id_matches_nothing(): + """The ``###_`` prefix is part of a light ec.rxns ID, as in MATLAB: the + bare cobra ID ``R2`` matches none of ``001_R2`` and ``002_R2``.""" ec = _ectestgem_light_ec_model() - updated = set_kcat_for_reactions(ec, ["R2"], 5.0, apply=False) - assert sorted(updated) == ["001_R2", "002_R2"] + updated = set_kcat_for_reactions(ec, ["R2"], 5.0) + assert updated == [] + assert _kcat_at(ec, "001_R2") == 0.0 + assert _kcat_at(ec, "002_R2") == 0.0 + + +def test_light_both_prefixed_ids_set_both_isozyme_rows(): + ec = _ectestgem_light_ec_model() + updated = set_kcat_for_reactions(ec, ["001_R2", "002_R2"], 5.0) + assert updated == ["001_R2", "002_R2"] assert _kcat_at(ec, "001_R2") == 5.0 assert _kcat_at(ec, "002_R2") == 5.0 assert _kcat_at(ec, "001_R2_REV") == 0.0 assert _kcat_at(ec, "002_R2_REV") == 0.0 -def test_light_unsuffixed_with_list_kcat_assigns_positionally(): - """Same as the full-layout case: an un-suffixed ID that expands to - multiple matches now accepts a per-match kcat list, assigned - positionally in ec.rxns order.""" +def test_light_prefixed_ids_with_list_kcat(): ec = _ectestgem_light_ec_model() - updated = set_kcat_for_reactions(ec, ["R2"], [1.0, 2.0], apply=False) + updated = set_kcat_for_reactions(ec, ["001_R2", "002_R2"], [1.0, 2.0]) assert updated == ["001_R2", "002_R2"] assert _kcat_at(ec, "001_R2") == 1.0 assert _kcat_at(ec, "002_R2") == 2.0 @@ -333,10 +354,10 @@ def test_light_apply_true_writes_prot_pool_coefficient(): def test_light_apply_picks_lowest_cost_isozyme(): - """If both isozymes of R2 get a kcat via the base-name shorthand, - apply_kcat must pick the cheapest (lowest MW_sum/kcat) one.""" + """If both isozymes of R2 get a kcat, apply_kcat must pick the + cheapest (lowest MW_sum/kcat) one.""" ec = _ectestgem_light_ec_model() - set_kcat_for_reactions(ec, ["R2"], 5.0, apply=True) + set_kcat_for_reactions(ec, ["001_R2", "002_R2"], 5.0, apply=True) r2 = ec.reactions.get_by_id("R2") coef = _get_s_coef(r2, "prot_pool") @@ -347,7 +368,6 @@ def test_light_apply_picks_lowest_cost_isozyme(): assert coef == pytest.approx(expected) -def test_light_unknown_rxn_id_raises(): +def test_light_unknown_rxn_id_is_ignored(): ec = _ectestgem_light_ec_model() - with pytest.raises(ValueError, match="matched no entries"): - set_kcat_for_reactions(ec, ["nonsense"], 1.0) + assert set_kcat_for_reactions(ec, ["nonsense"], 1.0) == [] diff --git a/tests/test_sigma_fitter.py b/tests/test_sigma_fitter.py index ef9dd8c..9d82325 100644 --- a/tests/test_sigma_fitter.py +++ b/tests/test_sigma_fitter.py @@ -100,16 +100,15 @@ def test_no_adapter_no_defaults_raises(tmp_path): def test_best_sigma_close_to_target_growth(tmp_path): """With p_tot=1, f=1, growth = 1000*sigma, so for growth_rate=10 the - exactly-optimal sigma is 0.01; the search should land on it with - near-zero error.""" + exactly-optimal sigma is 0.01; the grid scan should land on it with + zero error.""" adapter = _adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=10.0) model = _build_pool_only_model(adapter) result = fit_sigma(model) assert isinstance(result, SigmaFitterResult) # growth = 1000*sigma; for growth=10, sigma=0.01. assert result.sigma == pytest.approx(0.01, abs=1e-4) - # Error should be near zero for the best sigma (bisect resolves the - # crossing to ~1e-6 in sigma => ~0.01% relative growth error). + # The 0.01 grid contains the exactly-optimal sigma, so the error is zero. best_idx = int(np.argmin(result.error_grid)) assert result.error_grid[best_idx] < 0.1 @@ -160,8 +159,7 @@ def test_optimal_sigma_applied_to_model_at_end(tmp_path): def test_model_at_optimal_sigma_reaches_growth_rate(tmp_path): """After fitting, the model's LP optimum matches growth_rate (within - the bisect resolution: sigma is converged to ~1e-6, which here means - growth within ~1e-3 of the target on this linear model).""" + the grid resolution: here the grid contains the optimal sigma).""" adapter = _adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=50.0) model = _build_pool_only_model(adapter) fit_sigma(model) @@ -228,12 +226,11 @@ def test_sigma_grid_starts_above_zero_ends_at_one(tmp_path): assert result.sigma_grid[-1] == pytest.approx(1.0) -def test_default_method_is_bisect_with_fewer_solves(tmp_path): - """Default method is bisect: same answer as the grid search, with far - fewer evaluated sigmas.""" +def test_default_method_is_the_full_grid(tmp_path): + """Default method is the grid of sigmaFitter.m: every one of the 100 + sigmas is evaluated.""" adapter = _adapter(tmp_path, p_tot=1.0, f=1.0, gr_exp=10.0) model = _build_pool_only_model(adapter) - result = fit_sigma(model) # default == bisect - assert result.sigma == pytest.approx(0.01, abs=1e-4) - # Far fewer evaluations than the 100-point grid. - assert len(result.sigma_grid) < 30 + result = fit_sigma(model) # default == grid + assert result.sigma == pytest.approx(0.01, abs=1e-12) + assert len(result.sigma_grid) == 100 diff --git a/tutorials/full_ecModel/protocol.py b/tutorials/full_ecModel/protocol.py index 1bc50e3..f673bc0 100644 --- a/tutorials/full_ecModel/protocol.py +++ b/tutorials/full_ecModel/protocol.py @@ -302,7 +302,7 @@ sa = 2.15 # umol/min/mg protein kcat_per_sec = sa / 1000 / 60 * enz_mw print(f"Computed kcat for r_0079: {kcat_per_sec:.4f} /sec") -set_kcat_for_reactions(ec_model, ["r_0079"], kcat_per_sec) +set_kcat_for_reactions(ec_model, ["r_0079"], kcat_per_sec, apply=True) # %% save_ec_model(ec_model, "ecYeastGEM_stage3.yml", adapter=adapter) diff --git a/tutorials/light_ecModel/protocol.py b/tutorials/light_ecModel/protocol.py index cbafeff..d07d589 100644 --- a/tutorials/light_ecModel/protocol.py +++ b/tutorials/light_ecModel/protocol.py @@ -103,17 +103,17 @@ # pipeline). Here we set a couple of kcats by hand to demonstrate the # light dispatch. # -# ``set_kcat_for_reactions`` understands two ID shapes for light models: -# -# - **Base name** (``"R2"``): broadcasts the kcat to every isozyme row -# for that cobra reaction. Both ``001_R2`` and ``002_R2`` are updated. -# - **Explicit prefixed name** (``"001_R2"``): targets exactly one -# isozyme row, letting you set different kcats per isozyme. +# In a light model the ``ec.rxns`` IDs carry a ``###_`` prefix, which is +# part of the ID for ``set_kcat_for_reactions``, as in GECKO's +# ``setKcatForReactions``: ``"001_R2"`` targets exactly that isozyme row, +# so you can set different kcats per isozyme, and the bare cobra ID +# ``"R2"`` matches nothing. The call changes ``ec.kcat`` and ``ec.source`` +# only; ``apply_kcat_constraints`` writes the S matrix afterwards. # %% -set_kcat_for_reactions(ec_model, ["R2"], 5.0, apply=False) -set_kcat_for_reactions(ec_model, ["001_R3"], 10.0, apply=False) -set_kcat_for_reactions(ec_model, ["001_R5"], 20.0, apply=False) +set_kcat_for_reactions(ec_model, ["001_R2", "002_R2"], 5.0) +set_kcat_for_reactions(ec_model, ["001_R3"], 10.0) +set_kcat_for_reactions(ec_model, ["001_R5"], 20.0) print("ec.kcat after manual edits:") for rid in ec_model.ec.rxns: idx = ec_model.ec.rxns.index(rid)