Edge cases of geckopy agree with GECKO

CHANGELOG.md

@@ -2,12 +2,35 @@
 
 ## Unreleased
 
+* `ec_fseof` ranks the targets by the score `|slope| * |correlation|` rounded to 9 significant digits, with scores that are equal after rounding in model order, as GECKO's `ecFSEOF` does; scores that differ by solver noise used to decide the order.
+* `geckopy.kcat_tuning.evotune` exports `best_parsimonious`, `n_changed`, `parsimony_frontier` and `isozyme_tie_map`, the helpers behind the parsimony sweep and the tying of isozyme copies, which GECKO has as `bestParsimonious`, `nChanged`, `parsimonyFrontier` and `isozymeTieMap`.
+* `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""<TAB>""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number.
+* `read_dlkcat_output` drops a row whose kcat is infinite (`inf`), as it drops `NA` and `None`, so the kcat list holds finite numbers only.
+* `write_dlkcat_input` and `extract_enzyme_substrate_pairs` write the rows in the order of `model.ec.rxns` whatever the order of `ec_rxns`, and a reaction named twice once, as GECKO's `writeDLKcatInput` does with its logical vector.
+* `ECTree` raises `ValueError` for `k = NaN`, as GECKO's `ecTreeEstimate` does; it returned NaN estimates.
+* `kcats_from_ec_tree` with an empty `ec_rxns` returns no kcats instead of raising, as an empty selection means no reaction, and counts a reaction named twice once in the `ec_rxns` parameter of the source file.
+* `assign_standard_kcat` leaves out every reaction without a substrate or without a product (RAVEN's `getExchangeRxns`), such as `A + B =>`; it left out single-metabolite boundary reactions only.
+* `review_assignment` recognises the source labels `custom` and `standard` without regard to letter case, as MATLAB's `reviewKcatAssignment` does: a source `Custom` is checked for a duplicate of a database value and `STANDARD` is the fallback.
+* `find_met_smiles` skips a metabolite without a name with a logged warning and gives it no SMILES; it used to look the metabolite up by its id. When the cache file holds a name more than once, the last row is used (as before, now stated in the docstring).
+* `review_assignment`, `EcStats` and `ec_stats_from_brenda` are exported from `geckopy.gather_kcats` and `geckopy`, and listed in `__all__`.
+* `apply_kcat_constraints` takes cobra reaction ids as well as ec reaction ids (with the `###_` prefix) in `update_rxns` on a gecko-light model, as GECKO's `applyKcatConstraints` does for `model.rxns` ids; naming the cobra reaction updates it as a whole. An id that is both kinds or neither raises a `ValueError`.
+* `constrain_enz_concs` resets each processed usage reaction to a lower bound of 0 and a `prot_pool` coefficient of -1 before it sets the upper bound, as GECKO's `constrainEnzConcs` does; a prior lower bound above the new upper bound used to raise.
+* `fill_enz_concs` raises a `ValueError` for a negative abundance of a protein that is an enzyme of the model, as GECKO's `fillEnzConcs` does with the error `fillEnzConcs:negativeAbundance`; it used to skip the value with a warning.
+* `constrain_enz_concs` raises a `ValueError` for a negative concentration in `model.ec.concs` before it changes any bound, as GECKO's `constrainEnzConcs` does with the error `constrainEnzConcs:negativeConcentration`; it used to skip the enzyme with a warning. `remove_constraints=True` does not read the concentrations.
+* `load_flux_data` raises a `ValueError` for a row with a different number of cells than the header, as GECKO's `loadFluxData` does with the error `loadFluxData:raggedRow`; it used to pad the missing cells with NaN. A header `X (EX_a(e))` gives the metabolite `X` and the reaction id `EX_a(e)`: the reaction id is the text between the parenthesis that matches the final one and that one; a header with nested parentheses used to give the whole header for both.
+* `enzyme_usage` logs a warning that lists the enzymes without a `usage_prot_<enzyme>` reaction, which it leaves out of the result, as GECKO's `getEnzymeUsage` does with the warning `getEnzymeUsage:missingUsageReaction`.
+* `report_enzyme_usage` raises a `ValueError` for a negative `top_abs_usage`; it used to list all enzymes except the last ones, counted from the end of the ranking.
+* `save_ec_model` keeps the description of a model that has one (`metaData` of the file it was loaded from) and writes `Enzyme-constrained model of <id>` only when there is none; it used to overwrite it. GECKO's `saveEcModel` does the same.
+* `best_parsimonious` skips points with a NaN `rmse`, as MATLAB's `bestParsimonious` does with `min(rmses)`: a NaN at the first point of the frontier made every comparison false and raised `ValueError`. It still raises when every `rmse` is NaN.
+* `screen_kcat_leverage` keeps rows of equal rank in the order of the tie groups (a stable sort), as MATLAB's `screenKcatLeverage` does, and `compute_excarbon` warns when no metabolite of the model has a formula (every exchange then gets the carbon weight 1).
+* `cmaes_kcat_tuning` runs `max_generations` generations on a flat objective instead of stopping after the first (the flat-fitness criteria of `cma` do not end the search), as MATLAB's `cmaesKcatTuning` does.
 * `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_<id>` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. When no enzyme is at its usage bound, as on ecYeastGEM where the protein pool is the active constraint, a warning says so and the rows are in model order.
+* `save_ec_model` writes an ec reaction without enzymes with no `enzymes` entry, as GECKO's `saveEcModel` does (it wrote `enzymes: !!omap []`). `load_ec_model` reads a missing `enzymes` entry as no enzymes, and still reads `enzymes: !!omap` with no entries and `enzymes: !!omap []` from files written by earlier versions of GECKO, RAVEN and geckopy.
 * `get_conc_control_coeffs` takes the control coefficient by finite difference by default, as GECKO's `getConcControlCoeffs` does: it raises one usage bound by `fold_change` (2), re-solves and divides the rise of the optimum by the rise of the bound. The reduced-cost path, which reads the dual of the usage reaction and is unreliable at a degenerate vertex, is `method="shadow_price"`. The finite-difference path needs one solve per protein, so each iteration of `flexibilize_enz_concs`, which calls it, takes minutes on a genome-scale ecModel.
 * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model.
 * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason.
 * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets.
-* `pfba_enzymes` fixes the objective at 0.999999 of its optimum by default, as MATLAB's `getPfbaEnzymes` does through RAVEN's `solveLP`; the default was 1.0, so `objective / optimum` of the solution was 1.0 against MATLAB's 0.999999. Pass `fraction_of_optimum=1.0` for the old behaviour.
+* `pfba_enzymes` fixes the objective at 0.999999 of its optimum by default, as MATLAB's `getPfbaEnzymes` does through RAVEN's `solveLP`; a given `fraction_of_optimum` is multiplied by 0.999999 as well, so `fraction_of_optimum=0.9` fixes the objective at 0.8999991 of its optimum and the default 1.0 at 0.999999.
 * `model.tolerance` reaches HiGHS on the scipy solver. optlang's scipy interface ignored every tolerance, so a model set to `1e-9` (as the Gurobi and GLPK runs of the parity scenarios are) was solved at HiGHS's `1e-7`; the flux ranges of `ec_fva` on ecYeastGEM then differed from Gurobi's by up to 0.2 percent. The feasibility tolerance becomes `primal_feasibility_tolerance` and the optimality tolerance `dual_feasibility_tolerance`; both stay at HiGHS's default until set, and survive a copy or a pickle of the model. HiGHS ends some LPs with numerical difficulties (status 4) at `1e-9`; those are solved again at its default, and a solve that still ends so has the status `numeric` instead of raising `KeyError`.
 * A scipy `EcModel` is copied and pickled into the sparse problem. optlang pickles a model as JSON and restores it by adding the variables and constraints to a dense problem, so every worker process of `ec_fva` (and `cobra.Model.copy`, `copy.deepcopy`) spent tens of minutes rebuilding the matrix of a genome-scale ecModel; ecYeastGEM now restores in 1.6 s and copies in 2.2 s. The solver of a repaired or converted scipy model is a `geckopy.solver_support._SparseModel`.
 * `load_ec_model` reads the file on glpk and converts the ecModel to scipy when scipy is cobra's default solver, as `load_conventional_gem` does, and setting `model.solver = "scipy"` on an `EcModel` that is not on scipy builds the sparse problem of `geckopy.solver_support.switch_to_scipy`. optlang's own clone grows the dense coefficient matrix row by row: loading ecYeastGEM (8001 reactions) on scipy did not finish in 25 minutes and is now a matter of seconds. The scenarios that load a genome-scale ecModel on scipy spent most of their hours there. `geckopy.solver_support.is_scipy_interface` tells whether a solver name or interface is scipy.

docs/yaml_format.md

@@ -122,7 +122,7 @@ entry here. Fields:
 | `source` | string | optional | Where the kcat came from (`"brenda"`, `"dlkcat"`, `"manual"`, ...). Default `""` |
 | `notes` | string | optional | Free-form note. Default `""` |
 | `eccodes` | string \| list | optional | One EC code as a string (`"1.1.2.4"`), or a list when several apply. Default `""` |
-| `enzymes` | mapping `enzyme_id -> stoich` | yes | Which enzymes catalyse this reaction, and the subunit count of each. Every key must appear in `ec-enzymes` |
+| `enzymes` | mapping `enzyme_id -> stoich` | optional | Which enzymes catalyse this reaction, and the subunit count of each. Every key must appear in `ec-enzymes`. A reaction without enzymes (such as a `standard` kcat row) has no `enzymes` entry; the loader also reads the `enzymes: !!omap` and `enzymes: !!omap []` forms that earlier versions of GECKO, RAVEN and raven-toolbox wrote for it |
 
 Example:
 

src/geckopy/__init__.py

@@ -73,10 +73,12 @@
 
 # Gather-kcats
 from geckopy.gather_kcats import (
+    EcStats,
     apply_kcat_list,
     assign_kcats,
     assign_standard_kcat,
     collect_kcats,
+    ec_stats_from_brenda,
     fetch_open_kinetics_predictor,
     fuzzy_kcat_matching,
     get_standard_kcat,  # deprecated; alias of assign_standard_kcat
@@ -91,6 +93,7 @@
     read_dlkcat_output,
     read_kcat_source,
     remove_standard_kcat,
+    review_assignment,
     run_dlkcat,
     save_kcat_source,
     select_kcat_value,  # deprecated; alias of apply_kcat_list
@@ -175,6 +178,7 @@
     # Core data + builder + per-enzyme accessor
     "EcData",
     "EcModel",
+    "EcStats",
     "Enzyme",
     "make_ec_model",
     # Databases (loaders + dataclasses)
@@ -209,6 +213,7 @@
     "assign_kcats",
     "assign_standard_kcat",
     "collect_kcats",
+    "ec_stats_from_brenda",
     "fetch_open_kinetics_predictor",
     "fuzzy_kcat_matching",
     "kcats_from_brenda",
@@ -222,6 +227,7 @@
     "read_dlkcat_output",
     "read_kcat_source",
     "remove_standard_kcat",
+    "review_assignment",
     "run_dlkcat",
     "save_kcat_source",
     "submit_open_kinetics_predictor",

src/geckopy/databases/flux_data.py

@@ -4,7 +4,6 @@
 from __future__ import annotations
 
 import logging
-import re
 from dataclasses import dataclass, field
 from pathlib import Path
 from typing import Optional
@@ -15,7 +14,6 @@
 
 
 # `<met_name> (<rxn_id>)` -> capture met_name and rxn_id.
-_EXCH_HEADER_RE = re.compile(r"^(.*)\s*\(([^()]+)\)\s*$")
 
 
 @dataclass
@@ -81,6 +79,11 @@ def load_flux_data(path: str | Path) -> FluxData:
     * ``evotuneRMSEweight`` (float): weights for evotune kcat tuning.
     * ``source`` (str): free-text describing the data origin.
 
+    An exchange column header has the form ``<met> (<rxn>)``, split at
+    the parenthesis that matches the final closing one, so a reaction id
+    may contain parentheses (``X (EX_a(e))`` gives ``X`` and
+    ``EX_a(e)``).
+
     Both columns are dropped from the matrix and surfaced on the
     returned ``FluxData`` as ``evotune_rmse_weight`` and ``source``.
 
@@ -101,7 +104,8 @@ def load_flux_data(path: str | Path) -> FluxData:
     FileNotFoundError
         If ``path`` does not exist.
     ValueError
-        If the file has no header row or no exchange columns.
+        If the file has no header row, no exchange columns, or a row
+        with a different number of cells than the header.
     """
     path = Path(path)
     if not path.is_file():
@@ -123,6 +127,12 @@ def load_flux_data(path: str | Path) -> FluxData:
 
     header = rows[0]
     data = rows[1:]
+    for line_no, r in enumerate(data, start=2):
+        if len(r) != len(header):
+            raise ValueError(
+                f"{path}: row {line_no} has {len(r)} cell(s), the header "
+                f"has {len(header)}."
+            )
 
     # Pull out optional columns first; collect indices to drop.
     optional_indices: list[int] = []
@@ -134,15 +144,12 @@ def load_flux_data(path: str | Path) -> FluxData:
 
     if tuning_idx is not None:
         tuning_weights = np.array(
-            [_parse_float(r[tuning_idx]) if tuning_idx < len(r) else float("nan")
-             for r in data],
+            [_parse_float(r[tuning_idx]) for r in data],
             dtype=float,
         )
         optional_indices.append(tuning_idx)
     if source_idx is not None:
-        source_values = [
-            r[source_idx] if source_idx < len(r) else "" for r in data
-        ]
+        source_values = [r[source_idx] for r in data]
         optional_indices.append(source_idx)
 
     # Build the kept-column index list (preserving order).
@@ -158,22 +165,13 @@ def load_flux_data(path: str | Path) -> FluxData:
     exch_headers = kept_header[3:]
     exch_mets, exch_rxn_ids = _split_exch_headers(exch_headers)
 
-    conds = [r[kept[0]] if kept[0] < len(r) else "" for r in data]
-    p_tot = np.array(
-        [_parse_float(r[kept[1]]) if kept[1] < len(r) else float("nan")
-         for r in data],
-        dtype=float,
-    )
-    gr_rate = np.array(
-        [_parse_float(r[kept[2]]) if kept[2] < len(r) else float("nan")
-         for r in data],
-        dtype=float,
-    )
+    conds = [r[kept[0]] for r in data]
+    p_tot = np.array([_parse_float(r[kept[1]]) for r in data], dtype=float)
+    gr_rate = np.array([_parse_float(r[kept[2]]) for r in data], dtype=float)
     exch_fluxes = np.full((len(data), len(exch_headers)), np.nan, dtype=float)
     for i, r in enumerate(data):
         for ci, j in enumerate(kept[3:]):
-            if j < len(r):
-                exch_fluxes[i, ci] = _parse_float(r[j])
+            exch_fluxes[i, ci] = _parse_float(r[j])
 
     return FluxData(
         conds=conds,
@@ -206,6 +204,25 @@ def _parse_float(s: str) -> float:
         return float("nan")
 
 
+def _split_header(header: str) -> Optional[tuple[str, str]]:
+    """Split ``'<met> (<rxn>)'`` at the parenthesis that matches the final
+    closing one; return None for a header of another form."""
+    if not header.endswith(")"):
+        return None
+    depth = 0
+    for k in range(len(header) - 1, -1, -1):
+        if header[k] == ")":
+            depth += 1
+        elif header[k] == "(":
+            depth -= 1
+            if depth == 0:
+                rxn = header[k + 1:-1].strip()
+                if not rxn:
+                    return None
+                return header[:k].strip(), rxn
+    return None
+
+
 def _split_exch_headers(headers: list[str]) -> tuple[list[str], list[str]]:
     """Split each header `'<met> (<rxn>)'` into (met, rxn) pairs.
 
@@ -217,10 +234,10 @@ def _split_exch_headers(headers: list[str]) -> tuple[list[str], list[str]]:
     rxns: list[str] = []
     no_paren: list[str] = []
     for h in headers:
-        m = _EXCH_HEADER_RE.match(h.strip())
-        if m:
-            mets.append(m.group(1).strip())
-            rxns.append(m.group(2).strip())
+        parts = _split_header(h.strip())
+        if parts is not None:
+            mets.append(parts[0])
+            rxns.append(parts[1])
         else:
             mets.append(h.strip())
             rxns.append(h.strip())