diff --git a/CHANGELOG.md b/CHANGELOG.md index dd38ef1..228aa58 100644 --- a/CHANGELOG.md +++ b/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 (`""""`), 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_` 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 ` 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_` 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. diff --git a/docs/yaml_format.md b/docs/yaml_format.md index 283c58a..516ada8 100644 --- a/docs/yaml_format.md +++ b/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: diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index c5b5711..a23a0ff 100644 --- a/src/geckopy/__init__.py +++ b/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", diff --git a/src/geckopy/databases/flux_data.py b/src/geckopy/databases/flux_data.py index 3665544..4971512 100644 --- a/src/geckopy/databases/flux_data.py +++ b/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 @@ # ` ()` -> 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 `` ()``, 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 ``' ()'`` 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 `' ()'` 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()) diff --git a/src/geckopy/databases/pubchem.py b/src/geckopy/databases/pubchem.py index edcee5a..fc9ef0d 100644 --- a/src/geckopy/databases/pubchem.py +++ b/src/geckopy/databases/pubchem.py @@ -73,9 +73,11 @@ def find_met_smiles( For every metabolite in the model: 1. Skip if the name starts with ``prot_`` (enzyme pseudometabolite). - 2. Skip if ``annotation['smiles']`` is already set. - 3. Look up the metabolite name in the local TSV cache. - 4. If still not found, query PubChem; cache the result. + 2. Skip, with a warning, if the metabolite has no name. + 3. Skip if ``annotation['smiles']`` is already set. + 4. Look up the metabolite name in the local TSV cache. When the cache + holds a name more than once, the last row is used. + 5. If still not found, query PubChem; cache the result. The same metabolite name may appear on many ``cobra.Metabolite`` instances (different compartments). Lookup is done once per unique @@ -105,14 +107,24 @@ def find_met_smiles( # Build map: metabolite name -> list of metabolite objects sharing it. name_to_mets: dict[str, list] = {} + unnamed = [] for met in model.metabolites: - name = met.name or met.id - if not name or _PROT_PATTERN.match(name): + name = met.name + if not name: + unnamed.append(met.id) + continue + if _PROT_PATTERN.match(name): continue if met.annotation.get("smiles"): continue name_to_mets.setdefault(name, []).append(met) + if unnamed: + logger.warning( + "find_met_smiles: %d metabolite(s) have no name and get no SMILES: %s", + len(unnamed), ", ".join(unnamed), + ) + if not name_to_mets: logger.info("find_met_smiles: nothing to look up.") return diff --git a/src/geckopy/ec_model/pipeline/apply_kcat.py b/src/geckopy/ec_model/pipeline/apply_kcat.py index 989251c..221a8d1 100644 --- a/src/geckopy/ec_model/pipeline/apply_kcat.py +++ b/src/geckopy/ec_model/pipeline/apply_kcat.py @@ -65,7 +65,9 @@ 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. ``None`` or an + Reaction IDs (from ``model.ec.rxns``) to update. A gecko-light + model also takes the IDs of its cobra reactions (see + ``_apply_kcat_constraints_light``). ``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. @@ -216,10 +218,13 @@ 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`` 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. + IDs to update: ec.rxns IDs (with the ``###_`` prefix) or cobra + reaction IDs. ``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, or the + cobra reaction, recomputes the cheapest cost over all its rows. + An ID that is both an ec.rxns ID and a cobra reaction ID, or + neither, raises ``ValueError``. """ ec_rxn_ids = model.ec.rxns @@ -230,14 +235,34 @@ def _apply_kcat_constraints_light( selected_idx = np.arange(len(ec_rxn_ids)) else: id_to_idx = {rxn_id: i for i, rxn_id in enumerate(ec_rxn_ids)} - unknown = [r for r in update_rxns if r not in id_to_idx] + rows_of_cobra: dict[str, list[int]] = defaultdict(list) + for i, rxn_id in enumerate(ec_rxn_ids): + rows_of_cobra[split_light_rxn_id(rxn_id)[1]].append(i) + cobra_ids = {r.id for r in model.reactions} + ambiguous = [r for r in update_rxns if r in id_to_idx and r in cobra_ids] + if ambiguous: + raise ValueError( + f"update_rxns contains IDs that are both ec.rxns IDs and " + f"reaction IDs: {ambiguous[:5]}" + ) + unknown = [ + r for r in update_rxns + if r not in id_to_idx and r not in rows_of_cobra + ] if unknown: raise ValueError( - f"update_rxns contains IDs not present in ec.rxns: " - f"{unknown[:5]}" + f"update_rxns contains IDs present neither in ec.rxns nor " + f"as reactions of the model: {unknown[:5]}" ) selected_idx = np.array( - [id_to_idx[r] for r in update_rxns], dtype=int, + [ + i + for r in update_rxns + for i in ( + [id_to_idx[r]] if r in id_to_idx else rows_of_cobra[r] + ) + ], + dtype=int, ) if selected_idx.size == 0: diff --git a/src/geckopy/gather_kcats/__init__.py b/src/geckopy/gather_kcats/__init__.py index 795816b..f5688ab 100644 --- a/src/geckopy/gather_kcats/__init__.py +++ b/src/geckopy/gather_kcats/__init__.py @@ -18,6 +18,7 @@ ) from .read_dlkcat_output import read_dlkcat_output from .remove_standard_kcat import remove_standard_kcat +from .review import EcStats, ec_stats_from_brenda, review_assignment from .run_dlkcat import run_dlkcat from .select_kcat_value import ( apply_kcat_list, @@ -27,6 +28,7 @@ from .write_dlkcat_input import extract_enzyme_substrate_pairs, write_dlkcat_input __all__ = [ + "EcStats", "OKPClient", "OKPError", "apply_kcat_list", @@ -34,6 +36,7 @@ "assign_standard_kcat", "collect_kcats", "ec_rxn_keys", + "ec_stats_from_brenda", "extract_enzyme_substrate_pairs", "fetch_open_kinetics_predictor", "format_kcat_source", @@ -51,6 +54,7 @@ "read_kcat_source", "read_dlkcat_output", "remove_standard_kcat", + "review_assignment", "run_dlkcat", "save_kcat_source", "select_kcat_value", # deprecated; alias of apply_kcat_list diff --git a/src/geckopy/gather_kcats/ec_tree.py b/src/geckopy/gather_kcats/ec_tree.py index a5fce75..a58dd6f 100644 --- a/src/geckopy/gather_kcats/ec_tree.py +++ b/src/geckopy/gather_kcats/ec_tree.py @@ -21,8 +21,8 @@ class ECTree: """EC-tree estimates from a BRENDA measurement table (``ec``, ``log10_kcat``).""" def __init__(self, measurements: pd.DataFrame, k: float = 5.0): - if k <= 0: - raise ValueError("k must be positive") + if not k > 0: + raise ValueError("k must be a positive number") self.k = float(k) ec = measurements["ec"].astype(str) logk = measurements["log10_kcat"].to_numpy(dtype=float) diff --git a/src/geckopy/gather_kcats/get_standard_kcat.py b/src/geckopy/gather_kcats/get_standard_kcat.py index b57cf13..d8760f0 100644 --- a/src/geckopy/gather_kcats/get_standard_kcat.py +++ b/src/geckopy/gather_kcats/get_standard_kcat.py @@ -57,9 +57,10 @@ def assign_standard_kcat( (identified by ``model.ec.source == "standard"``) are stripped before re-applying. - Filtering: exchange, transport (same metabolite name across - compartments), pseudo (name contains "pseudoreaction"), SLIME - (name contains "SLIME rxn"), spontaneous (via + Filtering: exchange (a reaction without a substrate or without a + product), transport (same metabolite name across compartments), pseudo + (name contains "pseudoreaction", in any letter case), SLIME (name + contains "SLIME rxn", in any letter case), spontaneous (via ``adapter.get_spontaneous_reactions(model)``), and reactions listed in ``adapter.params.path/data/pseudoRxns.tsv`` are excluded from receiving the standard pseudoenzyme. @@ -237,6 +238,14 @@ def _find_reactions_missing_enzyme(model: "EcModel") -> list[str]: return no_gpr + has_gpr_no_ec +def _is_exchange(rxn: cobra.Reaction) -> bool: + """Whether a reaction only consumes or only produces metabolites, as + RAVEN's ``getExchangeRxns`` counts it: ``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 _detect_transport_reactions(model: "EcModel") -> set[str]: """A transport reaction has the same metabolite name appearing in different compartments.""" @@ -257,7 +266,7 @@ def _classify_reactions_to_ignore( ignore: set[str] = set(custom_rxns) for rxn in model.reactions: - if rxn.boundary: + if _is_exchange(rxn): ignore.add(rxn.id) name_lower = (rxn.name or "").lower() if "pseudoreaction" in name_lower or "slime rxn" in name_lower: diff --git a/src/geckopy/gather_kcats/kcat_sources.py b/src/geckopy/gather_kcats/kcat_sources.py index 41ef74c..f3437b1 100644 --- a/src/geckopy/gather_kcats/kcat_sources.py +++ b/src/geckopy/gather_kcats/kcat_sources.py @@ -344,8 +344,10 @@ def read_kcat_source(path: str | Path) -> tuple[pd.DataFrame, dict[str, str]]: The file is read as UTF-8 (a leading byte-order mark is ignored). The ``#`` lines before the column header are the header; ``#`` lines and - blank lines further down are skipped. Fields wrapped in double quotes, - as a spreadsheet may save them, are unwrapped. ``kcat`` and + blank lines further down are skipped, as are rows without a filled field + (a row of empty quoted fields included). Fields wrapped in double quotes, + as a spreadsheet may save them, are unwrapped; a field that opens a + double quote and does not close it raises ``ValueError``. ``kcat`` and ``log10_sd`` must be plain numbers (``12``, ``0.5``, ``1e-05``); other text, such as ``5,5`` with a decimal comma, is read as NaN, which :func:`load_kcat_source` reports. A row with more filled fields than @@ -389,9 +391,16 @@ def read_kcat_source(path: str | Path) -> tuple[pd.DataFrame, dict[str, str]]: raise ValueError(f"{path}, line {j + 1} has more fields than the {len(columns)} columns.") fields = fields[: len(columns)] fields += [""] * (len(columns) - len(fields)) - if not any(f.strip() for f in fields): - continue # a row of tabs only - rows.append([_unquote(f) for f in fields]) + for f in fields: + t = f.strip() + if t.startswith('"') and not (len(t) >= 2 and t.endswith('"')): + raise ValueError( + f"{path}, line {j + 1} has a field with an opening double quote and no closing one." + ) + values = [_unquote(f) for f in fields] + if not any(values): + continue # a row without a filled field + rows.append(values) line_numbers.append(j + 1) keep = [i for i, c in enumerate(columns) if c] records = pd.DataFrame([[r[i] for i in keep] for r in rows], columns=[columns[i] for i in keep], dtype=object) diff --git a/src/geckopy/gather_kcats/kcats_from.py b/src/geckopy/gather_kcats/kcats_from.py index f338343..0e90894 100644 --- a/src/geckopy/gather_kcats/kcats_from.py +++ b/src/geckopy/gather_kcats/kcats_from.py @@ -191,8 +191,9 @@ def kcats_from_ec_tree( Shrinkage constant: the number of measurements at which an EC's own median and its class count equally. ec_rxns - ``model.ec.rxns`` to estimate kcats for (default: all). With the - default ``save_as`` the file then only holds those reactions. + ``model.ec.rxns`` to estimate kcats for (default: all). An empty + list selects no reaction. With the default ``save_as`` the file then + only holds those reactions. measurements BRENDA measurement table (default: ``measurements.tsv`` of the adapter's BRENDA folder, see @@ -209,7 +210,7 @@ def kcats_from_ec_tree( (``n``) and the deepest EC level with measurements (``depth``). """ adapter = _adapter(model, adapter, "kcats_from_ec_tree needs the BRENDA folder and project folder") - ec_rxns = _rxn_subset(ec_rxns) + ec_rxns = None if ec_rxns is None else list(ec_rxns) if save: path = kcat_source_file(save_as, model=model, adapter=adapter) # fails before estimating _warn_subset(ec_rxns, model, path) @@ -232,7 +233,7 @@ def kcats_from_ec_tree( if save: parameters = {"k": k} if ec_rxns is not None: - parameters["ec_rxns"] = len(ec_rxns) + parameters["ec_rxns"] = len(set(ec_rxns)) path = save_kcat_source( model, kcat_list, save_as, adapter=adapter, generator="kcats_from_ec_tree", parameters=parameters, data=measurements.attrs.get("release", ""), diff --git a/src/geckopy/gather_kcats/read_dlkcat_output.py b/src/geckopy/gather_kcats/read_dlkcat_output.py index f91e46d..2795fc8 100644 --- a/src/geckopy/gather_kcats/read_dlkcat_output.py +++ b/src/geckopy/gather_kcats/read_dlkcat_output.py @@ -9,6 +9,7 @@ from pathlib import Path from typing import TYPE_CHECKING +import numpy as np import pandas as pd if TYPE_CHECKING: @@ -43,8 +44,9 @@ def read_dlkcat_output( The input file is the 6-column TSV that DLKcat writes (see ``run_dlkcat``), with one header line and columns ``rxn_id, gene, substrate, smiles, sequence, kcat``. Rows whose - kcat column is non-numeric (``"NA"``, empty, etc.) are silently - dropped before validation against the model. + kcat column is not a finite number (``"NA"``, ``"None"``, empty, + ``"inf"``, ``"1.2.3"``) are dropped before the reactions are validated + against the model. The output DataFrame matches the schema produced by ``fuzzy_kcat_matching`` so downstream functions @@ -87,8 +89,8 @@ def read_dlkcat_output( ValueError If the file contains no numeric kcat values; if ``strict=True`` and any substrate name from the file is not - in ``model.metabolites.`` (case-insensitive); or if any - ``rxn_id`` is not in ``model.ec.rxns``. + in ``model.metabolites.`` (case-insensitive); or if the + ``rxn_id`` of a row that is kept is not in ``model.ec.rxns``. """ file_path = Path(file_path) if not file_path.is_file(): @@ -104,6 +106,7 @@ def read_dlkcat_output( ) numeric_kcat = pd.to_numeric(df["kcat"], errors="coerce") + numeric_kcat = numeric_kcat.where(np.isfinite(numeric_kcat)) valid_mask = numeric_kcat.notna() if not valid_mask.any(): diff --git a/src/geckopy/gather_kcats/review.py b/src/geckopy/gather_kcats/review.py index 177ef90..1ebbf24 100644 --- a/src/geckopy/gather_kcats/review.py +++ b/src/geckopy/gather_kcats/review.py @@ -113,6 +113,9 @@ def review_assignment( flagged rows are a minority of the model's leverage. include Names from :data:`OPTIONAL_CHECKS` to report as well. + sources + The source labels ``custom`` and ``standard`` are recognised + without regard to letter case. Returns ------- @@ -176,7 +179,7 @@ def review_assignment( if kcats[i] < slow or kcats[i] > fast: checks.append("magnitude") detail.append(f"{kcats[i]:g} 1/s is outside {slow:g}-{fast:g}") - if (str(sources[i]) == "custom" and st is not None + if (str(sources[i]).lower() == "custom" and st is not None and round(float(kcats[i]), 6) in st.values): checks.append("custom-duplicate") detail.append("custom value equals a database value for its own EC") @@ -184,7 +187,7 @@ def review_assignment( if "tied-isozymes" in include: checks.append("tied-isozymes") detail.append(f"{len(idx)} isozyme copies share this value") - if str(sources[i]) == "standard": + if str(sources[i]).lower() == "standard": if "standard-fallback" not in include: continue checks.append("standard-fallback") diff --git a/src/geckopy/gather_kcats/write_dlkcat_input.py b/src/geckopy/gather_kcats/write_dlkcat_input.py index fc17cc9..2cd7737 100644 --- a/src/geckopy/gather_kcats/write_dlkcat_input.py +++ b/src/geckopy/gather_kcats/write_dlkcat_input.py @@ -82,7 +82,9 @@ def write_dlkcat_input( ``load_dlkcat_ignore_lists``). ec_rxns Optional iterable of reaction IDs (must be in - ``model.ec.rxns``). ``None`` means all. + ``model.ec.rxns``). ``None`` means all; an empty iterable selects + no reaction and writes an empty file. The rows follow the order of + ``model.ec.rxns``, whatever the order of ``ec_rxns``. only_with_smiles When True (default), drop rows whose substrate has no SMILES. When False, write ``"None"`` in the SMILES column @@ -150,7 +152,7 @@ def extract_enzyme_substrate_pairs( Loaded ``DLKcatIgnoreLists``. ec_rxns Optional iterable of reaction IDs (must be in ``model.ec.rxns``). - ``None`` means all. + ``None`` means all; the rows follow the order of ``model.ec.rxns``. only_with_smiles Drop rows with no SMILES when True; otherwise write ``"None"``. @@ -172,7 +174,7 @@ def extract_enzyme_substrate_pairs( f"{len(unknown)} reaction ID(s) in ec_rxns are not present " f"in model.ec.rxns (examples: {preview})" ) - ec_rxn_indices = [index_by_id[rid] for rid in ec_rxns_list] + ec_rxn_indices = sorted({index_by_id[rid] for rid in ec_rxns_list}) if not ec_rxn_indices: return pd.DataFrame(columns=_PAIR_COLUMNS) diff --git a/src/geckopy/kcat_tuning/evotune/__init__.py b/src/geckopy/kcat_tuning/evotune/__init__.py index 2a6e00f..7243d87 100644 --- a/src/geckopy/kcat_tuning/evotune/__init__.py +++ b/src/geckopy/kcat_tuning/evotune/__init__.py @@ -25,6 +25,7 @@ recommend_target_impact_share, sweep_tunable_mask_size, ) +from .parsimony import best_parsimonious, n_changed, parsimony_frontier from .priors import ( UNLABELLED_GROUP, build_sigma0_log, @@ -32,6 +33,7 @@ classify_kcat_sources, ) from .simulate import ConditionSimResult, simulate_evotune_dataset +from .tying import isozyme_tie_map from .tuning import ( EvotuneResult, cmaes_kcat_tuning, @@ -49,13 +51,17 @@ "ConditionSimResult", "MaskSizePoint", "evotune_distance", + "best_parsimonious", "build_sigma0_log", "classify_kcat_source", "classify_kcat_sources", "cmaes_kcat_tuning", "compute_excarbon", "dataset_rmse", + "isozyme_tie_map", "load_evotune_data", + "n_changed", + "parsimony_frontier", "recommend_target_impact_share", "screen_kcat_leverage", "select_tunable_mask", diff --git a/src/geckopy/kcat_tuning/evotune/distance.py b/src/geckopy/kcat_tuning/evotune/distance.py index f6cc757..ee84195 100644 --- a/src/geckopy/kcat_tuning/evotune/distance.py +++ b/src/geckopy/kcat_tuning/evotune/distance.py @@ -16,6 +16,7 @@ """ from __future__ import annotations +import warnings from typing import TYPE_CHECKING, Iterable import numpy as np @@ -49,7 +50,9 @@ def compute_excarbon( hand-rolled formula parser). A missing/unparseable formula, or a computed carbon count of exactly 0 (e.g. O2, water), both fall back to 1 -- matching MATLAB's ``Ematrix(isnan(Ematrix)) = 1`` and - ``excarbon(excarbon == 0) = 1`` respectively. + ``excarbon(excarbon == 0) = 1`` respectively. A model in which no + metabolite has a formula gives every exchange the weight 1 and + emits a warning. Parameters ---------- @@ -67,6 +70,12 @@ def compute_excarbon( dict[str, float] """ excarbon: dict[str, float] = {} + if not any(m.formula for m in model.metabolites): + warnings.warn( + "The model has no metabolite formulas; every exchange " + "reaction gets the carbon weight 1.", + UserWarning, stacklevel=2, + ) for rxn_id in rxn_ids: if rxn_id == bio_rxn_id: excarbon[rxn_id] = BIOMASS_CARBON_EQUIV diff --git a/src/geckopy/kcat_tuning/evotune/parsimony.py b/src/geckopy/kcat_tuning/evotune/parsimony.py index 7f3ffd3..3da3b74 100644 --- a/src/geckopy/kcat_tuning/evotune/parsimony.py +++ b/src/geckopy/kcat_tuning/evotune/parsimony.py @@ -239,10 +239,15 @@ def best_parsimonious( Defaults to 2%, which is the scale of the alternate-optimum noise measured on this model -- below that, two vectors are not - meaningfully different in fit, so the sparser one wins. + meaningfully different in fit, so the sparser one wins. A point with + a NaN RMSE is never chosen and does not set the best RMSE; when every + RMSE is NaN no point qualifies and ``ValueError`` is raised. """ if not points: raise ValueError("points is empty.") - best_rmse = min(p.rmse for p in points) + rmses = [p.rmse for p in points if not np.isnan(p.rmse)] + if not rmses: + raise ValueError("every point has a NaN rmse.") + best_rmse = min(rmses) within = [p for p in points if p.rmse <= best_rmse * (1.0 + tolerance)] return min(within, key=lambda p: p.n_changed) diff --git a/src/geckopy/kcat_tuning/evotune/tuning.py b/src/geckopy/kcat_tuning/evotune/tuning.py index 6c448f2..9585e68 100644 --- a/src/geckopy/kcat_tuning/evotune/tuning.py +++ b/src/geckopy/kcat_tuning/evotune/tuning.py @@ -185,6 +185,10 @@ def _resolve_context(model, adapter, params, evotune_data, bio_rxn, okp_method): return params, evotune_data, bio_rxn, okp_method +#: Termination criteria of ``cma`` that a flat objective triggers. +_FLAT_STOPS = frozenset({"tolfun", "tolfunhist", "tolflatfitness"}) + + def _tunable_context( model, params, evotune_data, bio_rxn, okp_method, tunable_mask, ): @@ -351,7 +355,8 @@ def screen_kcat_leverage( df = pd.DataFrame(rows) df["rank"] = df["leverage"] * df["sigma0_log"] - df = df.sort_values("rank", ascending=False, ignore_index=True) + df = df.sort_values( + "rank", ascending=False, ignore_index=True, kind="stable") total = float(df["rank"].sum()) df["cum_leverage_share"] = ( (df["rank"].cumsum() / total) if total > 0 else 0.0 @@ -412,7 +417,9 @@ def cmaes_kcat_tuning( Configuration comes entirely from ``EvotuneParams``: trust tiers, ``max_growth_weight``, ``prior_penalty_weight`` and ``tie_isozymes`` control the search, and ``max_generations``/``rmse_threshold`` are - its stopping conditions. + its stopping conditions. A generation whose candidates all score + the same does not end the search early; ``cma``'s other convergence + criteria (such as ``tolx``) still can. Parameters ---------- @@ -556,12 +563,20 @@ def _score_batch(kcat_matrix, pool): rmse_trace: list[float] = [] objective_trace: list[float] = [] generation = 0 - while not es.stop(): + flat = False + while True: + stop = es.stop() + # cma ends a search whose candidates all score the same + # (tolfun, tolfunhist, tolflatfitness); that is not a stopping + # condition here, so a flat generation does not end the run. + if stop and not (flat and set(stop) <= _FLAT_STOPS): + break asks = es.ask() xm = np.array(asks).T scores = _score_batch(expand(xm), pool) es.tell(asks, list(scores[:, 0])) generation += 1 + flat = bool(np.ptp(scores[:, 0]) == 0.0) j = int(np.argmin(scores[:, 0])) if scores[j, 0] < best_obj: diff --git a/src/geckopy/limit_proteins/constrain_enz_concs.py b/src/geckopy/limit_proteins/constrain_enz_concs.py index f7ce43c..a1e1442 100644 --- a/src/geckopy/limit_proteins/constrain_enz_concs.py +++ b/src/geckopy/limit_proteins/constrain_enz_concs.py @@ -35,9 +35,10 @@ def constrain_enz_concs( ``usage_prot_`` reaction's upper bound is set to that concentration. Enzymes with NaN concentrations have their usage reaction reset to the default ``upper_bound = 1000``, drawing - freely from ``prot_pool``. Only the reaction's upper bound is - changed; its stoichiometry (as set by - ``add_protein_usage_reactions``) is left untouched. + freely from ``prot_pool``. Each processed usage reaction is first + reset to ``lower_bound = 0``, ``upper_bound = 1000`` and a + ``prot_pool`` coefficient of -1, so a prior lower bound above the + new upper bound does not conflict with it. Ported from GECKO MATLAB: src/geckomat/limit_proteins/constrainEnzConcs.m. @@ -50,7 +51,8 @@ def constrain_enz_concs( machinery already in place. Mutated in place. remove_constraints When True, all usage reactions are reset to their default - ``upper_bound = 1000``, regardless of ``ec.concs``. + ``upper_bound = 1000`` (with ``lower_bound = 0`` and a + ``prot_pool`` coefficient of -1), regardless of ``ec.concs``. ``ec.concs`` itself is left unchanged. restrict_to If supplied, only the listed enzyme uniprot IDs are @@ -62,7 +64,9 @@ def constrain_enz_concs( Raises ------ ValueError - If the ``prot_pool`` metabolite is missing, or if any + If ``model.ec.concs`` holds a negative concentration for a + processed enzyme (unless ``remove_constraints`` is set), if the + ``prot_pool`` metabolite is missing, or if any ``usage_prot_`` reaction is missing for an enzyme in ``model.ec.enzymes`` (e.g. because the model is gecko-light, which has no usage reactions). @@ -94,20 +98,26 @@ def constrain_enz_concs( if enz in restrict_set ) - for i, enz in index_iter: + index_list = list(index_iter) + if not remove_constraints: + negative = [ + enz for i, enz in index_list if model.ec.concs[i] < 0 + ] + if negative: + raise ValueError( + "model.ec.concs has a negative concentration for " + + ", ".join(negative) + "." + ) + + pool_met = model.metabolites.get_by_id(POOL_ID) + for i, enz in index_list: rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enz}") - rxn.upper_bound = 1000.0 + rxn.bounds = (0.0, 1000.0) + if rxn.metabolites.get(pool_met) != -1.0: + rxn.add_metabolites({pool_met: -1.0}, combine=False) if remove_constraints: continue conc = model.ec.concs[i] if np.isnan(conc): continue - if conc < 0: - # A negative bound would make the usage reaction infeasible - # (ub < lb == 0). Skip it and leave the enzyme unconstrained. - logger.warning( - "constrain_enz_concs: negative concentration %g for %r; " - "leaving it unconstrained.", float(conc), enz, - ) - continue rxn.upper_bound = float(conc) diff --git a/src/geckopy/limit_proteins/fill_enz_concs.py b/src/geckopy/limit_proteins/fill_enz_concs.py index 4829b39..ad86e17 100644 --- a/src/geckopy/limit_proteins/fill_enz_concs.py +++ b/src/geckopy/limit_proteins/fill_enz_concs.py @@ -59,12 +59,13 @@ def fill_enz_concs( Raises ------ + ValueError + If a protein that is an enzyme of the model has a negative + abundance. IndexError If ``data_col`` is out of range for the supplied abundances. """ n_enz = model.ec.n_enzymes - model.ec.concs = np.full(n_enz, np.nan, dtype=float) - if prot_data.abundances.ndim == 2: if not 0 <= data_col < prot_data.abundances.shape[1]: raise IndexError( @@ -80,7 +81,9 @@ def fill_enz_concs( ) column = prot_data.abundances + concs = np.full(n_enz, np.nan, dtype=float) if n_enz == 0 or len(prot_data.uniprot_ids) == 0: + model.ec.concs = concs return enzyme_to_idx = {enz: i for i, enz in enumerate(model.ec.enzymes)} @@ -90,11 +93,8 @@ def fill_enz_concs( continue value = float(conc) if value < 0: - # A measured concentration can't be negative; skip it so it - # doesn't become an infeasible bound in constrain_enz_concs. - logger.warning( - "fill_enz_concs: negative concentration %g for %r; skipped.", - value, uid, + raise ValueError( + f"prot_data has a negative abundance ({value:g}) for {uid!r}." ) - continue - model.ec.concs[idx] = value + concs[idx] = value + model.ec.concs = concs diff --git a/src/geckopy/utilities/ec_fseof.py b/src/geckopy/utilities/ec_fseof.py index 9bebcee..2ba1283 100644 --- a/src/geckopy/utilities/ec_fseof.py +++ b/src/geckopy/utilities/ec_fseof.py @@ -116,7 +116,7 @@ def ec_fseof( rows filtered out of both ``scan`` and ``targets``. ``targets`` has an added boolean column ``transport`` and lists the reaction targets first, then the transport targets, each group in the - order of the score ``|slope| * |correlation|`` (equal scores in + order of the score ``|slope| * |correlation|`` rounded to 9 significant digits (equal scores in model order); a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting the ``prot_*`` ones. A caller that takes the top N rows gets the @@ -219,14 +219,19 @@ def _group_transport_targets(result: FSEOFResult, model: "EcModel") -> FSEOFResu """Add the ``transport`` column and list reaction targets before transport targets. As in ``ecFSEOF.m``, a reaction is a transport reaction when a metabolite name repeats - among its metabolites after the ``prot_*`` metabolites are dropped. The order within a - group is the order of ``result.targets``. + among its metabolites after the ``prot_*`` metabolites are dropped. Within a group the + reactions are ordered by their score rounded to 9 significant digits, equal scores in + model order. """ targets = result.targets.copy() targets["transport"] = [ _is_transport(model.reactions.get_by_id(rid)) for rid in targets["reaction"] ] - targets = targets.sort_values("transport", kind="stable", ignore_index=True) + position = {rxn.id: i for i, rxn in enumerate(model.reactions)} + targets["_model_order"] = [position[rid] for rid in targets["reaction"]] + targets["_score"] = [float(f"{x:.9g}") for x in targets["score"]] + targets = targets.sort_values(["transport", "_score", "_model_order"], ascending=[True, False, True], ignore_index=True) + targets = targets.drop(columns=["_model_order", "_score"]) return FSEOFResult(scan=result.scan, enforced=result.enforced, targets=targets) diff --git a/src/geckopy/utilities/enzyme_usage.py b/src/geckopy/utilities/enzyme_usage.py index 328d462..ac780dc 100644 --- a/src/geckopy/utilities/enzyme_usage.py +++ b/src/geckopy/utilities/enzyme_usage.py @@ -5,6 +5,7 @@ """ from __future__ import annotations +import logging from dataclasses import dataclass, field from typing import TYPE_CHECKING, Mapping, Union @@ -17,6 +18,8 @@ from ..ec_model.constants import USAGE_PREFIX +logger = logging.getLogger(__name__) + @dataclass class EnzymeUsageResult: @@ -83,6 +86,9 @@ def enzyme_usage( When False, enzymes with zero absolute usage are dropped from the result. + An enzyme without a ``usage_prot_`` reaction in the model is + left out of the result, with a logged warning. + Raises ------ NotImplementedError @@ -98,9 +104,11 @@ def enzyme_usage( abs_usage_list: list[float] = [] ub_list: list[float] = [] + left_out: list[str] = [] for enzyme in model.ec.enzymes: rxn_id = f"{USAGE_PREFIX}{enzyme}" if rxn_id not in cobra_rxn_ids: + left_out.append(enzyme) continue flux = float(_lookup_flux(fluxes, rxn_id)) rxn_ub = float(model.reactions.get_by_id(rxn_id).upper_bound) @@ -108,6 +116,12 @@ def enzyme_usage( abs_usage_list.append(abs(flux)) ub_list.append(rxn_ub) + if left_out: + logger.warning( + "enzyme_usage: the model has no usage reaction for the enzymes " + "%s; they are left out.", ", ".join(left_out), + ) + abs_usage_arr = np.array(abs_usage_list, dtype=float) ub_arr = np.array(ub_list, dtype=float) with np.errstate(divide="ignore", invalid="ignore"): diff --git a/src/geckopy/utilities/load_ec_model.py b/src/geckopy/utilities/load_ec_model.py index 18cc53a..7d5a918 100644 --- a/src/geckopy/utilities/load_ec_model.py +++ b/src/geckopy/utilities/load_ec_model.py @@ -23,7 +23,9 @@ 4. **Tolerated input**: a file that the YAML reader rejects is read again after the two forms that GECKO's ``loadEcModel`` (through RAVEN's ``readYAMLmodel``) accepts are rewritten as valid YAML: ``enzymes: !!omap`` - with no entries, as GECKO writes it for a reaction without enzymes, and a + with no entries, as earlier GECKO and RAVEN versions wrote it for a + reaction without enzymes (current writers omit the ``enzymes`` entry, which + is read as no enzymes), and a value in double quotes that holds unescaped double quotes, such as ``""text""``, which RAVEN 2 wrote for a value that has quotes of its own. The value loaded is the one MATLAB loads: the outer pair of quotes is @@ -145,9 +147,9 @@ def count_turned(message, category, filename, lineno, file=None, line=None): ec_model.ec = ec if on_fast_solver: switch_to_scipy(ec_model) - # raven-toolbox turns the reactions with a negative lower bound; GECKO - # turns those whose stoichiometry has the GECKO 3 sign, which also - # covers a usage reaction fixed to zero. + # raven-toolbox turns the reactions whose stoichiometry has the GECKO 3 + # sign; the same rule is applied here to a model that was not read from + # YAML. turned += len(upgrade_ec_model(ec_model)) if turned: warnings.warn( diff --git a/src/geckopy/utilities/pfba_enzymes.py b/src/geckopy/utilities/pfba_enzymes.py index 39581f7..66be1c9 100644 --- a/src/geckopy/utilities/pfba_enzymes.py +++ b/src/geckopy/utilities/pfba_enzymes.py @@ -33,22 +33,23 @@ _PFBA_OBJ_NAME = "_pfba_enzymes_objective" -# RAVEN's solveLP fixes the objective at this fraction of its optimum in the flux-minimising -# stage (getPfbaEnzymes calls it with minFlux 1). A solver reports the optimum to within its -# feasibility tolerance, so a second solve that demands the full value can be infeasible. -_PFBA_FRACTION = 0.999999 +# RAVEN's solveLP multiplies the fraction of the optimum by this margin in the +# flux-minimising stage (getPfbaEnzymes calls it with minFlux 1). A solver reports the +# optimum to within its feasibility tolerance, so a second solve that demands the full +# value can be infeasible. +_PFBA_MARGIN = 0.999999 def pfba_enzymes( model: "EcModel", objective: Union[cobra.Reaction, str, None] = None, *, - fraction_of_optimum: float = _PFBA_FRACTION, + fraction_of_optimum: float = 1.0, ) -> cobra.Solution: """Parsimonious FBA minimising the L1 norm of ``usage_prot_*`` fluxes. Sets the objective (if given), fixes it at - ``fraction_of_optimum * optimal_value``, then minimises + ``fraction_of_optimum * 0.999999 * optimal_value``, then minimises ``sum(usage_prot_*.forward_variable)``. Usage reactions are forward-only in the GECKO 3 layout (``lb=0``); if a user has flipped one to allow reverse flux, its ``reverse_variable`` is @@ -69,10 +70,12 @@ def pfba_enzymes( used as-is. fraction_of_optimum Fraction of the objective's optimal value to fix as a - constraint before minimising enzyme usage. The default, - 0.999999, is the fraction MATLAB's ``getPfbaEnzymes`` uses - (RAVEN's ``solveLP``); a lower value trades some objective - value for potentially lower enzyme usage. + constraint before minimising enzyme usage (default 1.0). The + fraction is multiplied by 0.999999, the margin of RAVEN's + ``solveLP`` that MATLAB's ``getPfbaEnzymes`` applies, so the + default fixes the objective at 0.999999 of its optimum and + 0.9 fixes it at 0.8999991. A lower value trades some + objective value for potentially lower enzyme usage. Returns ------- @@ -105,7 +108,9 @@ def pfba_enzymes( f"Objective is already {_PFBA_OBJ_NAME}; " "pfba_enzymes appears to be nested." ) - cobra.util.fix_objective_as_constraint(m, fraction=fraction_of_optimum) + cobra.util.fix_objective_as_constraint( + m, fraction=fraction_of_optimum * _PFBA_MARGIN, + ) usage_rxns = [r for r in m.reactions if r.id.startswith(USAGE_PREFIX)] variables = [] diff --git a/src/geckopy/utilities/report_enzyme_usage.py b/src/geckopy/utilities/report_enzyme_usage.py index 9526e5c..65af600 100644 --- a/src/geckopy/utilities/report_enzyme_usage.py +++ b/src/geckopy/utilities/report_enzyme_usage.py @@ -108,7 +108,16 @@ def report_enzyme_usage( Returns ------- EnzymeUsageReport + + Raises + ------ + ValueError + If ``top_abs_usage`` is negative or NaN. """ + if not top_abs_usage >= 0: + raise ValueError( + f"top_abs_usage must be a non-negative number, got {top_abs_usage!r}." + ) pool_rxn = model.reactions.get_by_id(POOL_EXCHANGE_ID) total_pool = float(pool_rxn.upper_bound) diff --git a/src/geckopy/utilities/save_ec_model.py b/src/geckopy/utilities/save_ec_model.py index f49c97a..4561153 100644 --- a/src/geckopy/utilities/save_ec_model.py +++ b/src/geckopy/utilities/save_ec_model.py @@ -73,7 +73,9 @@ def save_ec_model( Returns ------- Path - The absolute path of the file that was written. + The absolute path of the file that was written. An ec reaction + without enzymes (a zero row in ``model.ec.rxn_enz_mat``) is written + without an ``enzymes`` entry. Raises ------ @@ -139,8 +141,8 @@ def _build_metadata(model: EcModel) -> dict: format the file follows (as ``project_format`` in ``model_adapter.toml``, written as text as GECKO's ``saveEcModel`` writes it); ``generator``, the release that wrote the file; and a - description string of the form ``"Enzyme-constrained model of "``. - GECKO's ``saveEcModel`` writes the same ``gecko_format`` and + description string of the form ``"Enzyme-constrained model of "``, + written only when the model has no description of its own. GECKO's ``saveEcModel`` writes the same ``gecko_format`` and ``generator``. """ try: @@ -152,5 +154,6 @@ def _build_metadata(model: EcModel) -> dict: meta["date"] = meta.get("date") or date.today().isoformat() meta["gecko_format"] = "4" meta["generator"] = f"geckopy {geckopy_version}" - meta["description"] = f"Enzyme-constrained model of {model.id}" + if not meta.get("description"): + meta["description"] = f"Enzyme-constrained model of {model.id}" return meta diff --git a/tests/test_add_new_rxns_to_ec.py b/tests/test_add_new_rxns_to_ec.py index 7818750..8ade448 100644 --- a/tests/test_add_new_rxns_to_ec.py +++ b/tests/test_add_new_rxns_to_ec.py @@ -318,3 +318,25 @@ def test_multiple_new_rxns_added(tmp_path): assert sorted(result.enz_added) == ["Y", "Z"] assert "R_AY" in model.ec.rxns assert "R_BZ" in model.ec.rxns + + +# --------------------------------------------------------------------------- # +# Genes that already have enzymes +# --------------------------------------------------------------------------- # + +def test_reaction_with_genes_that_already_have_enzymes_is_added(tmp_path): + adapter = _adapter(tmp_path) + model = _build_base_model(adapter) + A_c = model.metabolites.get_by_id("A_c") + B_c = cobra.Metabolite("B_c", compartment="c") + model.add_metabolites([B_c]) + n_enzymes = len(model.ec.enzymes) + rxn = _make_rxn( + "R_AX", metabolites={A_c: -1.0, B_c: 1.0}, gpr="g_X", + ) + + result = add_new_rxns_to_ec(model, [rxn], []) + + assert result.enz_added == [] + assert len(model.ec.enzymes) == n_enzymes + assert "R_AX" in model.ec.rxns diff --git a/tests/test_apply_complex_data.py b/tests/test_apply_complex_data.py index 29f5cb7..8f26a34 100644 --- a/tests/test_apply_complex_data.py +++ b/tests/test_apply_complex_data.py @@ -283,3 +283,20 @@ def test_missing_path_and_no_complex_data_raises(tmp_path): with pytest.raises(FileNotFoundError): apply_complex_data(ec_model) + + +def test_partial_zero_stoichiometry_counts_the_subunit_as_one(): + """A subunit with stoichiometry 0 has no stated stoichiometry: it is a + subunit with stoichiometry 1, so the complex still matches.""" + ec_model = _ectestgem_ec_model() + data = [ComplexPortalEntry( + complex_id="R2Compl", + name="x", + species="t", + gene_names=["G1", "G2"], + protein_ids=["P1", "P2"], + stoichiometry=[2, 0], + )] + apply_complex_data(ec_model, complex_data=data, apply=False) + assert _stoich(ec_model, "R2_EXP_1", "P1") == 2.0 + assert _stoich(ec_model, "R2_EXP_1", "P2") == 1.0 diff --git a/tests/test_apply_kcat.py b/tests/test_apply_kcat.py index 73e87f5..c18e1f1 100644 --- a/tests/test_apply_kcat.py +++ b/tests/test_apply_kcat.py @@ -324,10 +324,31 @@ def test_light_setting_kcat_to_zero_and_reapplying_clears_old_constraint(): def test_light_update_rxns_unknown_id_raises(): ec = _ectestgem_light_model() - with pytest.raises(ValueError, match="not present in ec.rxns"): + with pytest.raises(ValueError, match="neither in ec.rxns nor"): apply_kcat_constraints(ec, update_rxns=["999_nonsense"]) +def test_light_update_rxns_takes_cobra_ids_and_ec_ids(): + ec = _ectestgem_light_model() + ec.ec.kcat[ec.ec.rxns.index("001_R3")] = 10.0 + r3 = ec.reactions.get_by_id("R3") + apply_kcat_constraints(ec, update_rxns=["R3"]) + coef = _get_s_coef(r3, "prot_pool") + assert coef < 0 + apply_kcat_constraints(ec, update_rxns=["001_R5"]) + assert _get_s_coef(r3, "prot_pool") == coef + ec.ec.kcat[ec.ec.rxns.index("001_R3")] = 0.0 + apply_kcat_constraints(ec, update_rxns=["001_R3"]) + assert _get_s_coef(r3, "prot_pool") == 0.0 + + +def test_light_update_rxns_id_of_both_kinds_raises(): + ec = _ectestgem_light_model() + ec.reactions.get_by_id("R5").id = "001_R3" + with pytest.raises(ValueError, match="both"): + apply_kcat_constraints(ec, update_rxns=["001_R3"]) + + def test_light_warns_when_no_valid_kcats(): """When every selected reaction has kcat=0 across all isozymes, emit the same warning the full branch does and leave the model unchanged.""" diff --git a/tests/test_calculate_f_factor.py b/tests/test_calculate_f_factor.py index cae1d40..cc23daf 100644 --- a/tests/test_calculate_f_factor.py +++ b/tests/test_calculate_f_factor.py @@ -178,3 +178,14 @@ def test_integration_with_pax_db_loader(tmp_path): f = calculate_f_factor(model, pd_data) # P1+P2 = 80, total = 110. assert f == pytest.approx(80.0 / 110.0) + + +# --------------------------------------------------------------------------- # +# Explicit empty enzyme list +# --------------------------------------------------------------------------- # + +def test_explicit_empty_enzyme_list_means_no_enzymes(): + model = _ec_model(["P1", "P2"]) + pd = _prot_data([("P1", 100.0), ("P2", 100.0)]) + assert calculate_f_factor(model, pd, enzymes=[]) == 0.0 + assert calculate_f_factor(model, pd) == pytest.approx(1.0) diff --git a/tests/test_constrain_enz_concs.py b/tests/test_constrain_enz_concs.py index 9089435..25045d1 100644 --- a/tests/test_constrain_enz_concs.py +++ b/tests/test_constrain_enz_concs.py @@ -215,3 +215,36 @@ def test_integration_with_fill_enz_concs(): assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 4.0 assert model.reactions.get_by_id("usage_prot_P2").upper_bound == 1000.0 assert model.reactions.get_by_id("usage_prot_P3").upper_bound == 9.0 + + +# --------------------------------------------------------------------------- # +# Negative concentration +# --------------------------------------------------------------------------- # + +def test_negative_conc_raises_before_changing_the_model(): + model = _ec_model_with_pool(["P1", "P2"], concs=[3.0, -0.5]) + with pytest.raises(ValueError, match="negative concentration"): + constrain_enz_concs(model) + assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 1000.0 + assert model.reactions.get_by_id("usage_prot_P2").upper_bound == 1000.0 + + +def test_negative_conc_is_not_read_when_removing_constraints(): + model = _ec_model_with_pool(["P1"], concs=[-0.5]) + constrain_enz_concs(model, remove_constraints=True) + assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 1000.0 + + +# --------------------------------------------------------------------------- # +# Prior bounds and coefficient +# --------------------------------------------------------------------------- # + +def test_prior_lower_bound_above_the_new_upper_bound_is_reset(): + model = _ec_model_with_pool(["P1"], concs=[0.01]) + rxn = model.reactions.get_by_id("usage_prot_P1") + rxn.lower_bound = 0.05 + pool = model.metabolites.get_by_id("prot_pool") + rxn.add_metabolites({pool: -2.0}, combine=False) + constrain_enz_concs(model) + assert rxn.bounds == (0.0, 0.01) + assert rxn.metabolites[pool] == -1.0 diff --git a/tests/test_ec_fseof.py b/tests/test_ec_fseof.py index 79048fa..6c988fc 100644 --- a/tests/test_ec_fseof.py +++ b/tests/test_ec_fseof.py @@ -415,3 +415,30 @@ def spy(*args, **kwargs): result = ec_fseof(_build_fseof_model(_adapter(tmp_path)), "EX_PROD", n_steps=4) assert seen["equality"] is True assert "transport" in result.targets.columns + + +def test_scores_equal_to_nine_digits_keep_model_order(): + """Scores that differ only below 9 significant digits count as equal and are listed in + model order.""" + import importlib + + import pandas as pd + + module = importlib.import_module("geckopy.utilities.ec_fseof") + model = cobra.Model("m") + met = cobra.Metabolite("a_c", name="a", compartment="c") + reactions = [] + for rid in ("R1", "R2", "R3"): + rxn = cobra.Reaction(rid) + rxn.add_metabolites({met: -1.0}) + reactions.append(rxn) + model.add_reactions(reactions) + result = FSEOFResult( + scan=pd.DataFrame(), + enforced=[1.0], + targets=pd.DataFrame({"reaction": ["R3", "R2", "R1"], "score": [2.0 + 1e-12, 2.0, 2.5]}), + ) + + grouped = module._group_transport_targets(result, model).targets + + assert list(grouped["reaction"]) == ["R1", "R2", "R3"] diff --git a/tests/test_enzyme_usage.py b/tests/test_enzyme_usage.py index b6b2888..4b3869d 100644 --- a/tests/test_enzyme_usage.py +++ b/tests/test_enzyme_usage.py @@ -189,3 +189,16 @@ def test_fluxes_carried_through_to_result(): result = enzyme_usage(model, sol.fluxes) # Same Series object (or equivalent values). assert result.fluxes is sol.fluxes + + +# --------------------------------------------------------------------------- # +# Enzyme without a usage reaction +# --------------------------------------------------------------------------- # + +def test_enzyme_without_usage_reaction_is_left_out_with_a_warning(caplog): + model = _build_toy(enzyme_ub=5.0) + model.ec.enzymes.append("F") + with caplog.at_level("WARNING"): + result = enzyme_usage(model, {"usage_prot_E": 3.0}) + assert result.prot_id == ["E"] + assert "no usage reaction for the enzymes F" in caplog.text diff --git a/tests/test_evotune.py b/tests/test_evotune.py index b926671..fc73710 100644 --- a/tests/test_evotune.py +++ b/tests/test_evotune.py @@ -443,3 +443,41 @@ def test_result_reports_its_own_corrections(): assert rows[0].source == "brenda" assert rows[0].fold_change == pytest.approx(10.0) assert rows[0].cumulative_share == pytest.approx(3.0 / 4.01) + + +def test_screen_orders_rows_of_equal_rank_as_the_tie_groups(tmp_path, monkeypatch): + from geckopy.kcat_tuning.evotune import tuning + + monkeypatch.setattr(tuning, "_score_kcat_vector", lambda *a, **k: (1.0, 1.0)) + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + screen = screen_kcat_leverage( + model, adapter=adapter, evotune_data=_evotune_data(), n_proc=1, + ) + assert list(screen["leverage"]) == [0.0, 0.0] + assert list(screen["rxn_id"]) == list(model.ec.rxns) + + +def test_cmaes_tuning_runs_all_generations_on_a_flat_objective(tmp_path, monkeypatch): + from geckopy.kcat_tuning.evotune import tuning + + monkeypatch.setattr(tuning, "_score_kcat_vector", lambda *a, **k: (1.0, 1.0)) + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + params = EvotuneParams(max_generations=6, rmse_threshold=-1.0) + result = cmaes_kcat_tuning( + model, adapter=adapter, params=params, evotune_data=_evotune_data(), + n_proc=1, seed=0, verbose=False, + ) + assert result.n_generations == 6 + + +def test_screen_warns_when_the_model_has_no_formulas(tmp_path): + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + for met in model.metabolites: + met.formula = None + with pytest.warns(UserWarning, match="no metabolite formulas"): + screen_kcat_leverage( + model, adapter=adapter, evotune_data=_evotune_data(), n_proc=1, + ) diff --git a/tests/test_evotune_parsimony.py b/tests/test_evotune_parsimony.py index ddb3e5a..3d360da 100644 --- a/tests/test_evotune_parsimony.py +++ b/tests/test_evotune_parsimony.py @@ -72,6 +72,20 @@ def test_best_parsimonious_prefers_fewer_changes_within_tolerance(): assert best_parsimonious(pts, tolerance=0.0).n_changed == 4000 +def test_best_parsimonious_skips_nan_rmse(): + from geckopy.kcat_tuning.evotune.parsimony import FrontierPoint + + pts = [ + FrontierPoint(0.0, np.array([1.0]), 4000, 7.0, float("nan")), + FrontierPoint(3.0, np.array([1.0]), 3000, 6.0, 0.905), + FrontierPoint(6.0, np.array([1.0]), 1000, 3.0, 1.400), + ] + assert best_parsimonious(pts).n_changed == 3000 + nan_only = [FrontierPoint(0.0, np.array([1.0]), 1, 1.0, float("nan"))] + with pytest.raises(ValueError): + best_parsimonious(nan_only) + + def test_fold_change_is_symmetric_and_at_least_one(): from geckopy.kcat_tuning.evotune.parsimony import fold_change diff --git a/tests/test_fill_enz_concs.py b/tests/test_fill_enz_concs.py index 02c5b7b..069c6be 100644 --- a/tests/test_fill_enz_concs.py +++ b/tests/test_fill_enz_concs.py @@ -191,3 +191,22 @@ def test_integration_with_load_pax_db(tmp_path): fill_enz_concs(model, prot) # P2 -> index 0 (60), P1 -> index 1 (20). np.testing.assert_array_equal(model.ec.concs, [60.0, 20.0]) + + +# --------------------------------------------------------------------------- # +# Negative abundance +# --------------------------------------------------------------------------- # + +def test_negative_abundance_of_an_enzyme_raises(): + model = _ec_model(["P1", "P2"], initial_concs=[1.0, 2.0]) + pd = _prot_data([("P1", 5.0), ("P2", -0.5)]) + with pytest.raises(ValueError, match="negative abundance"): + fill_enz_concs(model, pd) + np.testing.assert_array_equal(model.ec.concs, [1.0, 2.0]) + + +def test_negative_abundance_of_a_protein_outside_the_model_is_not_read(): + model = _ec_model(["P1"]) + pd = _prot_data([("P1", 5.0), ("X1", -1.0)]) + fill_enz_concs(model, pd) + np.testing.assert_array_equal(model.ec.concs, [5.0]) diff --git a/tests/test_kcat_contract.py b/tests/test_kcat_contract.py new file mode 100644 index 0000000..f2e260b --- /dev/null +++ b/tests/test_kcat_contract.py @@ -0,0 +1,143 @@ +"""Tests for the behaviour GECKO's kcatContractTests check as well.""" +import shutil +from pathlib import Path + +import cobra +import numpy as np +import pandas as pd +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.databases import UniprotDB, load_dlkcat_ignore_lists +from geckopy.gather_kcats import ( + assign_standard_kcat, + kcats_from_ec_tree, + read_dlkcat_output, + read_kcat_source, + write_dlkcat_input, +) +from geckopy.gather_kcats.ec_tree import ECTree +from geckopy.get_enzyme_data.ec_from_gem import fill_eccodes_from_gem + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + + +def _project(tmp_path: Path, extra_reactions: bool = False, light: bool = False): + project = tmp_path / "ecTestGEM" + shutil.copytree(EXAMPLE_DIR, project) + brenda = project / "data" / "brenda" + brenda.mkdir() + for name in ("kcat.tsv", "sa.tsv", "mw.tsv"): + shutil.copy(project / "data" / name, brenda / name) + adapter = ModelAdapter.from_folder(project) + conv = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + if extra_reactions: + m1 = conv.metabolites.get_by_id("m1c") + m2 = conv.metabolites.get_by_id("m2c") + for rxn_id, name, stoich in (("R6", "R6", {m1: -1.0, m2: 1.0}), + ("R7", "R7", {m1: -1.0, m2: -1.0}), + ("R8", "Biomass Pseudoreaction", {m1: -1.0, m2: 1.0})): + rxn = cobra.Reaction(rxn_id, name=name) + rxn.lower_bound = 0.0 + rxn.upper_bound = 1000.0 + rxn.add_metabolites(stoich) + conv.add_reactions([rxn]) + model = make_ec_model(conv, adapter, gecko_light=light) + fill_eccodes_from_gem(model) + return model, adapter, project + + +def _uniprot() -> UniprotDB: + mws = [100.0, 200.0, 300.0] + return UniprotDB(ids=["P000", "P001", "P002"], genes=["g0", "g1", "g2"], eccodes=["1.1.1.1"] * 3, + mw=np.array(mws), sequences=[""] * 3) + + +def test_read_kcat_source_quoted_rows(tmp_path): + """A row of empty quoted fields is blank and a quoted comment is a comment; a + field that opens a quote and does not close it is an error.""" + file = tmp_path / "odd.tsv" + head = ["# source: odd", "rxn\tenzymes\tkcat\tsource", "R3\tP4\t25\todd"] + file.write_text("\n".join(head + ['"# a quoted comment"', '""\t""\t""\t""', "R5\tP5\t30\todd"]) + "\n", + encoding="utf-8") + records, _ = read_kcat_source(file) + assert list(records["kcat"]) == [25.0, 30.0] + assert list(records["line"]) == [3, 6] + file.write_text("\n".join(head + ['"# a comment that was quoted']) + "\n", encoding="utf-8") + with pytest.raises(ValueError, match="line 4"): + read_kcat_source(file) + file.write_text("\n".join(head[:2] + ['"R3\tP4\t25\todd']) + "\n", encoding="utf-8") + with pytest.raises(ValueError, match="no closing one"): + read_kcat_source(file) + + +def test_read_dlkcat_output_kcat_fields(tmp_path): + """A kcat field that is not a finite number drops the row before the reactions + are checked; a file without a number raises.""" + model, _, _ = _project(tmp_path) + file = tmp_path / "DLKcat_out.tsv" + head = "rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat" + + def write(*rows): + file.write_text("\n".join([head] + [f"{r}\tG1\t{s}\tsmiles\tseq\t{k}" for r, s, k in rows]) + "\n", + encoding="utf-8") + + write(("R2_EXP_1", "m1", "12.5"), ("R2_EXP_2", "m1", "inf"), ("R2_EXP_2", "m1", "1.2.3")) + out = read_dlkcat_output(model, file) + assert list(out["rxn_id"]) == ["R2_EXP_1"] and list(out["kcat"]) == [12.5] + write(("R2_EXP_1", "m1", "None"), ("R2_EXP_2", "m1", "None")) + with pytest.raises(ValueError, match="no numeric kcat"): + read_dlkcat_output(model, file) + write(("R2_EXP_1", "m1", "12.5"), ("NOTAREACTION_EXP_1", "m1", "NA")) + assert list(read_dlkcat_output(model, file)["rxn_id"]) == ["R2_EXP_1"] + write(("R2_EXP_1", "m1", "12.5"), ("NOTAREACTION_EXP_1", "m1", "7")) + with pytest.raises(ValueError, match="not present in model.ec.rxns"): + read_dlkcat_output(model, file) + + +def test_write_dlkcat_input_selection(tmp_path): + """An empty selection writes an empty file; the rows follow model.ec.rxns.""" + model, _, _ = _project(tmp_path) + ignore = load_dlkcat_ignore_lists() + target = tmp_path / "DLKcat_in.tsv" + df = write_dlkcat_input(model, target, ignore, ec_rxns=[], overwrite=True, only_with_smiles=False) + assert df.empty and target.read_text(encoding="utf-8") == "" + df = write_dlkcat_input(model, target, ignore, ec_rxns=["R5", "R3", "R5"], overwrite=True, + only_with_smiles=False) + wanted = [r for r in model.ec.rxns if r in {"R3", "R5"}] + assert list(dict.fromkeys(df["rxn_id"])) == wanted + + +def test_ec_tree_rejects_nan_k(): + table = pd.DataFrame({"ec": ["1.1.1.1", "1.1.1.1"], "log10_kcat": [1.0, 2.0]}) + with pytest.raises(ValueError, match="positive"): + ECTree(table, k=float("nan")) + + +def test_kcats_from_ec_tree_subset(tmp_path): + """An empty ec_rxns selects no reaction; a reaction named twice gives one row.""" + model, adapter, _ = _project(tmp_path) + ecs = sorted({e.strip() for c in model.ec.eccodes for e in str(c or "").split(";") if e.strip()}) + table = pd.DataFrame([(ec, float(i + 1)) for i, ec in enumerate(ecs) for _ in range(2)], + columns=["ec", "log10_kcat"]) + table["kcat"] = 10.0 ** table["log10_kcat"] + table.attrs["release"] = "BRENDA release test" + with_ec = [r for r, c in zip(model.ec.rxns, model.ec.eccodes) if str(c or "").strip()] + assert len(with_ec) >= 2 + none = kcats_from_ec_tree(model, measurements=table, adapter=adapter, ec_rxns=[], save=False) + assert none.empty + twice = kcats_from_ec_tree(model, measurements=table, adapter=adapter, save=False, + ec_rxns=[with_ec[0], with_ec[1], with_ec[1]]) + assert list(twice["rxn_id"]) == with_ec[:2] + + +def test_standard_kcat_exchange_and_pseudoreactions(tmp_path): + """R7 consumes two metabolites and produces none: an exchange reaction. R8 is + named like a pseudoreaction in another letter case. R6 is an ordinary reaction + without a GPR.""" + model, _, _ = _project(tmp_path, extra_reactions=True) + model.ec.kcat[:] = 10.0 + assign_standard_kcat(model, _uniprot()) + assert "R6" in model.ec.rxns + assert "R7" not in model.ec.rxns + assert "R8" not in model.ec.rxns diff --git a/tests/test_load_ec_model.py b/tests/test_load_ec_model.py index 01febe0..cad9903 100644 --- a/tests/test_load_ec_model.py +++ b/tests/test_load_ec_model.py @@ -554,7 +554,7 @@ def read_with_warning(path): # Input that MATLAB's readYAMLmodel reads and a YAML parser rejects # --------------------------------------------------------------------------- # -# The ec reactions as GECKO writes them: R1_EXP_1 has no enzymes (`enzymes: !!omap` +# The ec reactions as earlier GECKO versions wrote them: R1_EXP_1 has no enzymes (`enzymes: !!omap` # with no entries, as for a `standard` row), R3 is the last row before the next section. # The notes of R1_EXP_1 and R3 hold quotes that RAVEN 2 did not escape. _GECKO_FORM = """\ diff --git a/tests/test_load_flux_data.py b/tests/test_load_flux_data.py index baa5188..7c32541 100644 --- a/tests/test_load_flux_data.py +++ b/tests/test_load_flux_data.py @@ -289,3 +289,29 @@ def test_loaded_data_can_feed_apply_flux_data_constraints(tmp_path): apply_flux_data_constraints(model, fd) assert model.reactions.get_by_id("EX_glc").lower_bound == -5.0 assert model.reactions.get_by_id("biomass").upper_bound == pytest.approx(0.4) + + +# --------------------------------------------------------------------------- # +# Malformed rows and headers +# --------------------------------------------------------------------------- # + +def test_ragged_row_raises(tmp_path): + p = tmp_path / "ragged.tsv" + _write_tsv(p, [ + ["condition", "Ptot", "grRate", "Glc (EX_a)", "O2 (EX_b)"], + ["c1", "0.5", "0.4", "-1", "2"], + ["c2", "0.5", "0.4", "-1"], + ]) + with pytest.raises(ValueError, match="row 3 has 4 cell"): + load_flux_data(p) + + +def test_nested_parentheses_in_reaction_id(tmp_path): + p = tmp_path / "nested.tsv" + _write_tsv(p, [ + ["condition", "Ptot", "grRate", "X (EX_a(e))", "O2 (EX_b)"], + ["c1", "0.5", "0.4", "-1", "2"], + ]) + fd = load_flux_data(p) + assert fd.exch_mets == ["X", "O2"] + assert fd.exch_rxn_ids == ["EX_a(e)", "EX_b"] diff --git a/tests/test_pfba_enzymes.py b/tests/test_pfba_enzymes.py index fd5aa2f..2259839 100644 --- a/tests/test_pfba_enzymes.py +++ b/tests/test_pfba_enzymes.py @@ -45,11 +45,18 @@ def test_returns_solution(): assert sol.status == "optimal" -def test_objective_preserved_at_full_fraction(): +def test_full_fraction_keeps_0_999999_of_the_optimum(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() sol = pfba_enzymes(model, fraction_of_optimum=1.0) - assert sol.fluxes["R3"] == pytest.approx(fba_optimum, abs=1e-6) + assert sol.fluxes["R3"] / fba_optimum == pytest.approx(0.999999, abs=2e-7) + + +def test_given_fraction_is_multiplied_by_the_margin(): + model = _ectestgem_ec_model_with_kcats() + fba_optimum = model.slim_optimize() + sol = pfba_enzymes(model, fraction_of_optimum=0.9) + assert sol.fluxes["R3"] / fba_optimum == pytest.approx(0.9 * 0.999999, abs=2e-7) def test_default_fixes_the_objective_at_0_999999_of_the_optimum(): @@ -65,7 +72,7 @@ def test_objective_preserved_at_partial_fraction(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() sol = pfba_enzymes(model, fraction_of_optimum=0.5) - assert sol.fluxes["R3"] >= 0.5 * fba_optimum - 1e-6 + assert sol.fluxes["R3"] >= 0.5 * 0.999999 * fba_optimum - 1e-6 def test_total_usage_le_plain_fba(): diff --git a/tests/test_pubchem.py b/tests/test_pubchem.py index 4ba33d2..3a4e41d 100644 --- a/tests/test_pubchem.py +++ b/tests/test_pubchem.py @@ -368,3 +368,22 @@ def test_met_smiles_accepts_the_shape_a_loaded_model_carries(): del met.annotation["smiles"] assert met_smiles(met) == "" + + +def test_duplicate_cache_name_uses_the_last_row(monkeypatch, tmp_path): + cache_path = tmp_path / "smilesDB.tsv" + cache_path.write_text("glucose\tFIRST\nglucose\tLAST\n") + model = _build_model([("glc_c", "glucose")]) + find_met_smiles(model, cache_path=cache_path) + assert _smiles(model, "glc_c") == "LAST" + + +def test_metabolite_without_name_is_skipped_with_a_warning(monkeypatch, tmp_path, caplog): + cache_path = tmp_path / "smilesDB.tsv" + cache_path.write_text("m1_c\tFROM_ID\nglucose\tCCO\n") + model = _build_model([("m1_c", ""), ("glc_c", "glucose")]) + with caplog.at_level(logging.WARNING): + find_met_smiles(model, cache_path=cache_path) + assert _smiles(model, "m1_c") == "" + assert _smiles(model, "glc_c") == "CCO" + assert any("m1_c" in r.getMessage() and "no name" in r.getMessage() for r in caplog.records) diff --git a/tests/test_relax_proteomics_greedy.py b/tests/test_relax_proteomics_greedy.py index 80a62aa..a7fdc1a 100644 --- a/tests/test_relax_proteomics_greedy.py +++ b/tests/test_relax_proteomics_greedy.py @@ -99,3 +99,12 @@ def test_max_iterations_raises(): relax_proteomics_greedy( model, minimal_growth=1e9, max_iterations=1, ) + + +def test_empty_enzyme_set_makes_no_enzyme_eligible(): + model, _ = _ec_with_tight_concs() + result = relax_proteomics_greedy( + model, minimal_growth=0.001, enzyme_set=[], max_iterations=20, + ) + assert result.relaxed == {} + assert result.trace == [] diff --git a/tests/test_report_enzyme_usage.py b/tests/test_report_enzyme_usage.py index 2c41082..3a6eda9 100644 --- a/tests/test_report_enzyme_usage.py +++ b/tests/test_report_enzyme_usage.py @@ -314,3 +314,14 @@ def test_top_abs_usage_inf_returns_all(): ud = enzyme_usage(model, sol.fluxes) rep = report_enzyme_usage(model, ud, top_abs_usage=float("inf")) assert len(rep.top_abs_usage) == 1 + + +# --------------------------------------------------------------------------- # +# Invalid top_abs_usage +# --------------------------------------------------------------------------- # + +def test_negative_top_abs_usage_raises(): + model = _build_single_enzyme_model() + ud = enzyme_usage(model, model.optimize().fluxes) + with pytest.raises(ValueError, match="top_abs_usage"): + report_enzyme_usage(model, ud, top_abs_usage=-2) diff --git a/tests/test_review_public_api.py b/tests/test_review_public_api.py new file mode 100644 index 0000000..90cab47 --- /dev/null +++ b/tests/test_review_public_api.py @@ -0,0 +1,53 @@ +"""The kcat review functions are public API of geckopy and geckopy.gather_kcats.""" +import geckopy +from geckopy import gather_kcats +from geckopy.databases import load_brenda_data + +_NAMES = ("EcStats", "ec_stats_from_brenda", "review_assignment") + + +def test_review_functions_are_exported_from_both_packages(): + for name in _NAMES: + assert name in geckopy.__all__ + assert name in gather_kcats.__all__ + assert getattr(geckopy, name) is getattr(gather_kcats, name) + + +def test_ec_stats_keep_the_case_of_the_ec_number(tmp_path): + (tmp_path / "kcat.tsv").write_text( + "# BRENDA release test\n" + "ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n" + "3.5.1.n3\tm1\torg\t40\t4\t2\t*\n" + "1.1.1.1\tm1\torg\t10\t2\t2\t*\n" + "1.1.1.1\tm2\torg\t6\t3\t2\t*\n", + encoding="utf-8", + ) + (tmp_path / "sa.tsv").write_text( + "# BRENDA release test\n" + "ec_code\tsubstrate\torganism\tsa_max\tsa_median\tn\treferences\n", + encoding="utf-8", + ) + (tmp_path / "mw.tsv").write_text( + "# BRENDA release test\nec_code\tsubstrate\torganism\tmw\tn\treferences\n", + encoding="utf-8", + ) + stats = geckopy.ec_stats_from_brenda(load_brenda_data(tmp_path, filter_outlier_kcats=False)) + assert sorted(stats) == ["1.1.1.1", "3.5.1.n3"] + assert stats["1.1.1.1"].kcat_max == 10.0 + assert stats["1.1.1.1"].kcat_median == 2.5 + assert stats["3.5.1.n3"].values == frozenset({4.0, 40.0}) + findings = geckopy.review_assignment( + ["R1"], [40.0], ["brenda"], ["3.5.1.n3"], ec_stats=stats, + ) + assert [f.checks for f in findings] == [("ec-maximum",)] + + +def test_review_recognises_custom_and_standard_without_regard_to_case(): + stats = {"1.1.1.1": geckopy.EcStats(n_kcat=1, kcat_max=2.0, kcat_median=2.0, values=frozenset({2.0}))} + found = geckopy.review_assignment( + ["R1", "R2"], [2.0, 20000.0], ["Custom", "STANDARD"], ["1.1.1.1", ""], ec_stats=stats, + include=("standard-fallback",), + ) + by_rxn = {f.rxn_id: f.checks for f in found} + assert "custom-duplicate" in by_rxn["R1"] + assert "standard-fallback" in by_rxn["R2"] diff --git a/tests/test_save_ec_model.py b/tests/test_save_ec_model.py index 9c89601..3c36817 100644 --- a/tests/test_save_ec_model.py +++ b/tests/test_save_ec_model.py @@ -1,4 +1,6 @@ """Tests for save_ec_model.""" +import os +import re from datetime import date from pathlib import Path @@ -201,6 +203,14 @@ def test_doc_top_level_keys(tmp_path): }.issubset(doc.keys()) +def test_existing_description_is_kept(tmp_path): + adapter = _adapter(tmp_path) + model = _build_full_ec_model(adapter) + model.notes["metaData"] = {"description": "my own description"} + path = save_ec_model(model) + assert _read_yaml(path)["metaData"]["description"] == "my own description" + + def test_metadata_includes_date_version_description(tmp_path): adapter = _adapter(tmp_path) model = _build_full_ec_model(adapter) @@ -413,3 +423,73 @@ def test_round_trip_with_gecko_light_true(tmp_path): path = save_ec_model(model) loaded = load_ec_model(path) assert loaded.ec.gecko_light is True + + +# --------------------------------------------------------------------------- # +# ec reaction without enzymes +# --------------------------------------------------------------------------- # + +def _build_model_with_enzymeless_rxn(adapter: ModelAdapter) -> EcModel: + """`_build_full_ec_model` with R2 as an ec reaction without enzymes.""" + model = _build_full_ec_model(adapter) + mat = model.ec.rxn_enz_mat.tolil() + mat[1, 1] = 0.0 + model.ec.rxn_enz_mat = mat.tocsr() + model.ec.rxn_enz_mat.eliminate_zeros() + model.ec.source[1] = "standard" + return model + + +def test_rxn_without_enzymes_has_no_enzymes_entry(tmp_path): + adapter = _adapter(tmp_path) + path = save_ec_model(_build_model_with_enzymeless_rxn(adapter)) + text = path.read_text(encoding="utf-8") + ec_rxns = text.split("- ec-rxns:")[1].split("- ec-enzymes:")[0] + assert ec_rxns.count("enzymes") == 1 + r2 = re.split(r'- id: "?R2"?\n', ec_rxns)[1] + assert "enzymes" not in r2 + assert "standard" in r2 + assert "!!omap []" not in text + doc = _read_yaml(path) + assert "enzymes" not in doc["ec-rxns"][1] + + +def test_round_trip_rxn_without_enzymes(tmp_path): + adapter = _adapter(tmp_path) + model = _build_model_with_enzymeless_rxn(adapter) + loaded = load_ec_model(save_ec_model(model)) + assert loaded.ec.rxns == ["R1", "R2"] + assert loaded.ec.source == ["brenda", "standard"] + np.testing.assert_array_equal(loaded.ec.kcat, model.ec.kcat) + np.testing.assert_array_equal( + loaded.ec.rxn_enz_mat.toarray(), [[1.0, 0.0], [0.0, 0.0]], + ) + assert loaded.ec.enzymes == ["E1", "E2"] + + +_LIGHT_ECMODEL = Path(os.environ.get("LIGHT_ECMODEL_YML", "/nonexistent")) + + +@pytest.mark.skipif( + not _LIGHT_ECMODEL.is_file(), + reason="set LIGHT_ECMODEL_YML to GECKO's tutorials/light_ecModel/models/ecModel.yml", +) +def test_light_tutorial_model_loses_only_the_empty_enzymes_entries(tmp_path): + """The light tutorial ecModel, written by GECKO 3, has `enzymes: !!omap` with + no entries for its reactions without enzymes. Loaded and saved again, the + file has no `enzymes` entry for them and the ec data is unchanged.""" + loaded = load_ec_model(_LIGHT_ECMODEL) + n_zero = int((loaded.ec.rxn_enz_mat.getnnz(axis=1) == 0).sum()) + assert n_zero == 52 + path = save_ec_model(loaded, tmp_path / "out.yml", adapter=None) + saved = path.read_text(encoding="utf-8") + saved_rxns = saved.split("- ec-rxns:")[1].split("- ec-enzymes:")[0] + assert "!!omap []" not in saved_rxns + assert saved_rxns.count("- id:") == len(loaded.ec.rxns) + assert saved_rxns.count("enzymes:") == len(loaded.ec.rxns) - n_zero + again = load_ec_model(path) + assert again.ec.rxns == loaded.ec.rxns + np.testing.assert_array_equal(again.ec.kcat, loaded.ec.kcat) + assert again.ec.source == loaded.ec.source + assert again.ec.enzymes == loaded.ec.enzymes + assert (again.ec.rxn_enz_mat != loaded.ec.rxn_enz_mat).nnz == 0