diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index ef078fc..747f59e 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -38,5 +38,11 @@ jobs: - name: Install run: pip install -e ".[dev,evotune]" + # ec_fseof needs the `equality` option of raven_toolbox.analysis.fseof, + # which is not in a release yet. Remove this step once geckopy requires + # a raven-toolbox release that has it. + - name: Install raven-toolbox from develop + run: pip install --upgrade --force-reinstall --no-deps "raven-toolbox @ git+https://github.com/CompMetEng/raven-toolbox.git@develop" + - name: Run pytest run: pytest -q diff --git a/CHANGELOG.md b/CHANGELOG.md index d2a9f5b..f96b82e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,22 @@ ## Unreleased +* `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. +* `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. +* `get_subset_ec_model` removes the metabolites and genes that no reaction of the subset uses, as the `removeGenes` and `removeReactions` calls of GECKO's `getSubsetEcModel` do, also those that were unused in the big ecModel. The HT29-GEM subset of the light Human-GEM ecModel had the pseudo-gene `standard` and six metabolites (`MAM00579m`, `MAM01009c`, `MAM01099c`, `MAM01420m`, `MAM02549c`, `MAM02579m`) more than MATLAB's. A context model equal to the big model gives a copy without its unused metabolites and genes. +* `load_ec_model` loads the two forms of ecModel file that GECKO's `loadEcModel` reads and the YAML parser of raven-toolbox rejects, with the values MATLAB loads. `- enzymes: !!omap` with no entries, as GECKO writes it for an ec reaction without enzymes (the `standard` rows), raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`; the row has no enzymes. A value in double quotes that holds unescaped double quotes, such as the notes `""text""` of the published eciYali5-GEM (RAVEN 2's `writeYAMLmodel` wraps a value in quotes without escaping the quotes in it, and GECKO 3 read the notes column of `customKcats.tsv` with `%s`, which keeps the quotes of the file), raised a parser error; the note is `"text"`, the outer pair of quotes removed and the inner ones kept. A file that the parser rejects is read again after these two forms are rewritten as YAML. +* `load_uniprot_tsv` and `ModelAdapter.get_uniprot_ids_from_table` read the cells as MATLAB's `textscan` with `%q` does, which the docstring of `load_uniprot_tsv` already claimed: a cell in double quotes is given without the quotes (`"1.1.1.25; 2.5.1.19"` gives `1.1.1.25; 2.5.1.19`), `""` inside a quoted cell is one quote, tabs and line breaks inside quotes belong to the cell, and leading spaces of a cell are dropped. A line with fewer cells is padded with empty cells and the cells of a line beyond the last column start a new row, as `textscan` does. The conversion table maps a gene that is listed twice to its first row, as `ismember` does, and a gene on a line without a second cell to an empty UniProt ID. `geckopy.databases.uniprot_loader.read_quoted_tsv` does the splitting. +* `assign_standard_kcat` compares the subsystem names of the reactions without regard to letter case, as MATLAB's `strcmpi` does: `Pentose phosphate pathway` and `Pentose Phosphate Pathway` are one subsystem, so their reactions count together toward the `threshold` and give one median. A reaction with an empty subsystem has none and gets the model-wide standard kcat. +* `get_subset_ec_model` removes the genes of the big ecModel that the context model lacks with `cobra.manipulation.remove_genes`, so a reaction that survives the removal (an isozyme remains, as with the rule `G1 and G2 or G3` without `G3`) loses the gene from its rule. Before, the gene left `model.genes` while the rule still named it, and removing that reaction later because the context model lacks it raised `ValueError: G3 not found`, as it did for the HT29-GEM subset of the light Human-GEM ecModel. In a gecko-light model the `001_` prefix of the ec reaction ids is no longer part of the id the ec rows are matched against: the subset used to drop every row of `ec.rxns`, `kcat`, `source`, `notes`, `eccodes` and `rxn_enz_mat`. A context model equal to the big model gives a copy, as before. +* The scipy solver (optlang's HiGHS interface) gives the results of glpk and gurobi. optlang renumbers the variable columns wrongly after a variable is removed (it stops at `len(constraints)` instead of `len(variables)`), so cobra read fluxes from the wrong columns after `make_ec_model`, `Reaction.delete` or any other removal: `ec_fva` returned 0 for every reaction of ecTestGEM although the LP optimum was 3272.7. An `EcModel` now repairs the problem whenever `model.solver` is read; `geckopy.solver_support.repair_scipy_problem(model)` does the same for a plain `cobra.Model`. The repaired problem stores the coefficients sparsely, because optlang copies a dense matrix on every added or removed variable or constraint, and it keeps the solution when `cobra.Model.optimize` sets the objective direction again, so `Reaction.flux` and `Enzyme.flux` work after `model.optimize()`. After an infeasible solve it returns NaN values and objective, so `model.optimize(raise_error=False)` returns a solution with status `infeasible` instead of raising a `TypeError`. +* With scipy as cobra's default solver, `load_conventional_gem` reads the model on glpk and `make_ec_model` builds the ecModel on glpk, and both convert the result to scipy (`geckopy.solver_support.switch_to_scipy`). Building the yeast-GEM ecModel on the scipy interface took more than 40 minutes in optlang's dense matrix operations without finishing. +* `Enzyme.shadow_price`, `get_enzyme_bottlenecks` and `relax_proteomics_greedy` work on the scipy solver, which has no constraint duals in optlang: `geckopy.solver_support.constraint_shadow_prices` reads the HiGHS row marginals of the last solve and uses them only if they are an optimum of the problem as it is now. `relax_proteomics_greedy` computes the shadow prices of all candidates once per iteration and takes the lowest ID among enzymes with equal absolute shadow price. +* `flexibilize_enz_concs` and `sensitivity_tuning` decide that the starting LP is infeasible from the solver status, not from the objective value, so they behave the same on every solver. Before, an infeasible start raised on Gurobi (no objective value) but not on GLPK (a finite objective value). `flexibilize_enz_concs` now returns an empty result for it on all solvers, as MATLAB's `flexibilizeEnzConcs` does, instead of raising on Gurobi. `sensitivity_tuning` raises on all solvers and stops with a warning when a later solve is not optimal. * `apply_kcat_constraints` treats an empty `update_rxns` as all reactions, as MATLAB's `applyKcatConstraints` does with an empty `updateRxns`; it used to do nothing. The callers that build the list from a filter (`set_kcat_for_reactions`, the evotune scoring and reset) diff --git a/src/geckopy/adapter/adapter.py b/src/geckopy/adapter/adapter.py index 30d8329..c3317bb 100644 --- a/src/geckopy/adapter/adapter.py +++ b/src/geckopy/adapter/adapter.py @@ -141,19 +141,22 @@ def get_uniprot_ids_from_table(self, model_genes: list[str]) -> list[str]: If the conversion table is absent, returns the input unchanged. The table is a two-column TSV with header row: model gene ID, - UniProt ID. + UniProt ID. The cells are read as MATLAB's ``textscan`` with ``%q`` + reads them (see :func:`geckopy.databases.uniprot_loader.read_quoted_tsv`), + so a cell in double quotes is given without the quotes. """ table_path = self.params.path / "data" / "uniprotConversion.tsv" if not table_path.is_file(): return list(model_genes) + from ..databases.uniprot_loader import read_quoted_tsv + + with open(table_path, "r", encoding="utf-8", newline="") as f: + text = f.read() + # A gene listed twice takes its first row, as MATLAB's ismember does. mapping: dict[str, str] = {} - with open(table_path, "r", encoding="utf-8") as f: - next(f) # skip header - for line in f: - parts = line.rstrip("\n").split("\t") - if len(parts) >= 2: - mapping[parts[0]] = parts[1] + for gene, uniprot_id in read_quoted_tsv(text, 2, header_lines=1): + mapping.setdefault(gene, uniprot_id) return [mapping.get(g, g) for g in model_genes] diff --git a/src/geckopy/databases/uniprot_loader.py b/src/geckopy/databases/uniprot_loader.py index c7a63a0..a9fd474 100644 --- a/src/geckopy/databases/uniprot_loader.py +++ b/src/geckopy/databases/uniprot_loader.py @@ -115,7 +115,9 @@ def load_uniprot_tsv( Notes ----- Empty cells (common for EC number and Sequence) become empty - strings, matching MATLAB ``textscan`` with ``%q`` behavior. + strings. Cells are read as :func:`read_quoted_tsv` reads them, as + MATLAB's ``textscan`` with ``%q`` does: the quotes around a cell + are dropped, so ``"1.1.1.25; 2.5.1.19"`` gives ``1.1.1.25; 2.5.1.19``. Duplicate entries raise, matching MATLAB's dispEM failure mode. """ path = Path(path) @@ -144,44 +146,36 @@ def load_uniprot_tsv( # expands to. source_row_ids: list[str] = [] - with open(path, "r", encoding="utf-8") as f: - header_line = f.readline() - if not header_line: - raise ValueError(f"{path} is empty.") - header_fields = header_line.rstrip("\n").split("\t") - if len(header_fields) < 5: - raise ValueError( - f"{path} header has {len(header_fields)} fields, expected 5." - ) + with open(path, "r", encoding="utf-8", newline="") as f: + text = f.read() + if not text: + raise ValueError(f"{path} is empty.") + header_end = _first_line_end(text) + header_fields = text[:header_end].rstrip("\r\n").split("\t") + if len(header_fields) < 5: + raise ValueError( + f"{path} header has {len(header_fields)} fields, expected 5." + ) + + for entry, gene_cell, ec, mass_str, seq in read_quoted_tsv(text, 5, header_lines=1): + # MATLAB str2double returns NaN for empty or non-numeric; mirror. + try: + mass_da = float(mass_str) if mass_str.strip() else float("nan") + except ValueError: + mass_da = float("nan") - for line_no, line in enumerate(f, start=2): - line = line.rstrip("\n") - if not line: - continue - parts = line.split("\t") - # Pad short lines so all five fields are always addressable. - while len(parts) < 5: - parts.append("") - entry, gene_cell, ec, mass_str, seq = parts[:5] - - # MATLAB str2double returns NaN for empty or non-numeric; mirror. - try: - mass_da = float(mass_str) if mass_str.strip() else float("nan") - except ValueError: - mass_da = float("nan") - - if split_gene_cells and gene_cell.strip(): - gene_names = gene_cell.split() - else: - gene_names = [gene_cell] - - source_row_ids.append(entry) - for gene_name in gene_names: - ids.append(entry) - genes.append(gene_name) - eccodes.append(ec) - mw_da.append(mass_da) - sequences.append(seq) + if split_gene_cells and gene_cell.strip(): + gene_names = gene_cell.split() + else: + gene_names = [gene_cell] + + source_row_ids.append(entry) + for gene_name in gene_names: + ids.append(entry) + genes.append(gene_name) + eccodes.append(ec) + mw_da.append(mass_da) + sequences.append(seq) _check_duplicate_ids(source_row_ids, path) @@ -214,3 +208,82 @@ def _check_duplicate_ids(source_row_ids: list[str], path: Path) -> None: f"Duplicate Entry values in {path}: {preview}. " f"Manually curate uniprot.tsv or adjust adapter parameters." ) + + +def _first_line_end(text: str, start: int = 0) -> int: + """Index just after the end of the line of ``text`` that begins at + ``start`` (the length of ``text`` for the last line).""" + for i in range(start, len(text)): + c = text[i] + if c == "\n": + return i + 1 + if c == "\r": + return i + 2 if text[i + 1:i + 2] == "\n" else i + 1 + return len(text) + + +def read_quoted_tsv(text: str, n_columns: int, header_lines: int = 0) -> list[list[str]]: + """Split tab-delimited ``text`` into rows of ``n_columns`` cells the way + MATLAB's ``textscan(fid, '%q ...', 'Delimiter', '\\t')`` does. + + A cell that starts with a double quote (after spaces) ends at the next + double quote that is not doubled; a doubled quote inside stands for one + quote, and tabs and line breaks inside are part of the cell. Text between + the closing quote and the next tab or line break is appended to the cell. + A quote that is never closed takes the rest of the text. In any other + cell the quotes are kept, and leading spaces are dropped. A line break + ends the row, so a short line is padded with empty cells; the cells of a + line beyond ``n_columns`` start a new row. A line without any text is + skipped. The first ``header_lines`` lines are skipped as they are, as + ``textscan`` does with ``'HeaderLines'``. + """ + rows: list[list[str]] = [] + row: list[str] = [] + n = len(text) + i = 0 + for _ in range(header_lines): + i = _first_line_end(text, i) + + def finish(cell: str, at_line_end: bool) -> None: + row.append(cell) + if len(row) == n_columns or at_line_end: + row.extend([""] * (n_columns - len(row))) + rows.append(row.copy()) + row.clear() + + while i < n: + while i < n and text[i] in " \b": + i += 1 + quoted = i < n and text[i] == '"' + cell: list[str] = [] + if quoted: + i += 1 + while i < n: + if text[i] == '"': + if text[i + 1:i + 2] == '"': + cell.append('"') + i += 2 + continue + i += 1 + break + cell.append(text[i]) + i += 1 + while i < n and text[i] not in "\t\r\n": + cell.append(text[i]) + i += 1 + value = "".join(cell) + if i >= n: + if value or quoted or row: + finish(value, True) + break + c = text[i] + if c == "\t": + i += 1 + finish(value, False) + if i >= n and row: + finish("", True) + continue + i += 2 if c == "\r" and text[i + 1:i + 2] == "\n" else 1 + if value or quoted or row: + finish(value, True) + return rows diff --git a/src/geckopy/ec_model/ec_model.py b/src/geckopy/ec_model/ec_model.py index 6aaf415..806b170 100644 --- a/src/geckopy/ec_model/ec_model.py +++ b/src/geckopy/ec_model/ec_model.py @@ -25,6 +25,7 @@ import cobra +from ..solver_support import is_scipy_interface, repair_scipy_problem, switch_to_scipy from .ec_data import EcData if TYPE_CHECKING: @@ -71,6 +72,24 @@ def __init__( from .enzyme import EnzymeView self.enzymes = EnzymeView(self) + @property + def solver(self): + """The optlang model. With the scipy interface, reading it also repairs the + problem's variable index map (see :mod:`geckopy.solver_support`), so values read + after any removal of variables come from the right columns.""" + repair_scipy_problem(self) + return cobra.Model.solver.fget(self) + + @solver.setter + def solver(self, value) -> None: + """Set the solver. Moving a model to the scipy interface goes through + :func:`geckopy.solver_support.switch_to_scipy`: optlang's clone grows the dense + coefficient matrix row by row, which takes hours for a genome-scale model.""" + if not self._contexts and is_scipy_interface(value) and not is_scipy_interface(self.problem): + switch_to_scipy(self) + else: + cobra.Model.solver.fset(self, value) + @classmethod def from_cobra( cls, diff --git a/src/geckopy/ec_model/enzyme.py b/src/geckopy/ec_model/enzyme.py index c44e322..1e6cddb 100644 --- a/src/geckopy/ec_model/enzyme.py +++ b/src/geckopy/ec_model/enzyme.py @@ -36,6 +36,7 @@ import cobra +from ..solver_support import constraint_shadow_prices from .constants import PROT_PREFIX, USAGE_PREFIX from .pipeline.populate_ec import split_light_rxn_id @@ -186,9 +187,9 @@ def upper_bound(self, value: float) -> None: @property def shadow_price(self) -> float: - """Dual of the prot_ mass-balance constraint.""" + """Dual of the prot_ mass-balance constraint, on every solver (scipy included).""" self._require_full_model("shadow_price") - return float(self._model.constraints[self.prot_metabolite_id].dual) + return float(constraint_shadow_prices(self._model, [self.prot_metabolite_id])[0]) @property def reactions(self) -> frozenset: diff --git a/src/geckopy/ec_model/make_ec_model.py b/src/geckopy/ec_model/make_ec_model.py index 5734201..8a4f3db 100644 --- a/src/geckopy/ec_model/make_ec_model.py +++ b/src/geckopy/ec_model/make_ec_model.py @@ -30,6 +30,7 @@ from raven_toolbox.utils.sort import sort_identifiers from ..databases import UniprotDB, load_uniprot_tsv +from ..solver_support import repair_scipy_problem, switch_to_scipy from .ec_model import EcModel from .pipeline import ( add_protein_pool_exchange_reaction, @@ -203,6 +204,15 @@ def make_ec_model( # (splitting reactions, adding pseudometabolites, ...), but the caller's # input GEM must be left untouched. model = model.copy() + # On the scipy solver the pipeline runs on glpk and the result is converted at the end: + # optlang grows the scipy matrix densely, so the thousands of edits of a genome-scale + # build take hours there (see geckopy.solver_support). + on_scipy = repair_scipy_problem(model) + if on_scipy: + try: + model.solver = "glpk" + except Exception: # glpk unavailable: build on scipy + on_scipy = False # Stages 1-4: preprocess on a plain cobra.Model. Stage 5 (isozyme # expansion) is full-only; light keeps isozymes singular and tracks @@ -279,6 +289,9 @@ def make_ec_model( more, ) + if on_scipy: + switch_to_scipy(ec_model) + # Catch any internal length/shape drift in the ec arrays before the # model is handed back (cheap, and surfaces build bugs early). ec_model.ec.validate() diff --git a/src/geckopy/ec_model/pipeline/apply_kcat.py b/src/geckopy/ec_model/pipeline/apply_kcat.py index 12b0557..989251c 100644 --- a/src/geckopy/ec_model/pipeline/apply_kcat.py +++ b/src/geckopy/ec_model/pipeline/apply_kcat.py @@ -9,6 +9,7 @@ """ from __future__ import annotations +import contextlib import warnings from collections import defaultdict from typing import TYPE_CHECKING @@ -19,10 +20,15 @@ from ..ec_model import EcModel +from ...solver_support import edit_on_fast_solver from ..constants import POOL_ID, PROT_PREFIX from .populate_ec import split_light_rxn_id +# Number of updated ec reactions from which a scipy model is edited on glpk. +_BULK_UPDATE = 100 + + def apply_kcat_constraints( model: "EcModel", update_rxns: list[str] | None = None, @@ -75,6 +81,17 @@ def apply_kcat_constraints( If every kcat to be applied is 0 or NaN. Old coefficients are still cleared in this case; nothing is written afterward. """ + # A scipy model is edited on glpk when many reactions change: the scipy interface + # rewrites the rows of a constraint for every changed coefficient. + many = update_rxns is None or len(update_rxns) == 0 or len(update_rxns) >= _BULK_UPDATE + with edit_on_fast_solver(model) if many else contextlib.nullcontext(): + _apply_kcat_constraints(model, update_rxns) + + +def _apply_kcat_constraints( + model: "EcModel", + update_rxns: list[str] | None, +) -> None: if model.ec.gecko_light: _apply_kcat_constraints_light(model, update_rxns) return diff --git a/src/geckopy/gather_kcats/get_standard_kcat.py b/src/geckopy/gather_kcats/get_standard_kcat.py index e65f7c2..b57cf13 100644 --- a/src/geckopy/gather_kcats/get_standard_kcat.py +++ b/src/geckopy/gather_kcats/get_standard_kcat.py @@ -49,7 +49,9 @@ def assign_standard_kcat( ``threshold`` reactions with a real kcat) are computed. Reactions that lack a GPR rule, or have a GPR but no entry in ``model.ec.rxns``, are added to ``model.ec`` with this standard - pseudoenzyme. + pseudoenzyme. Only the first subsystem of a reaction counts, and + subsystem names that differ only in letter case are one subsystem. + A reaction with an empty subsystem has none. The function is idempotent: previous "standard" entries (identified by ``model.ec.source == "standard"``) are stripped @@ -172,6 +174,13 @@ def _first_subsystem(rxn: cobra.Reaction) -> str: return sub.split(";")[0].strip() +def _subsystem_key(rxn: cobra.Reaction) -> str: + """The first subsystem of ``rxn`` without regard to letter case, as + MATLAB's ``strcmpi`` compares subsystem names; empty for a reaction + without a subsystem.""" + return _first_subsystem(rxn).lower() + + def _ec_rxn_to_cobra_id(model: "EcModel", ec_rxn_id: str) -> str: return ec_rxn_id[4:] if model.ec.gecko_light else ec_rxn_id @@ -197,7 +206,7 @@ def _compute_subsystem_kcats( rxn = model.reactions.get_by_id(cobra_rxn_id) except KeyError: continue - sub = _first_subsystem(rxn) + sub = _subsystem_key(rxn) if not sub: continue sub_kcats[sub].append(float(kcat)) @@ -379,7 +388,7 @@ def _assign_standard_kcat_to_missing( kcats_to_add: list[float] = [] for rxn_id in rxn_ids: rxn = model.reactions.get_by_id(rxn_id) - sub = _first_subsystem(rxn) + sub = _subsystem_key(rxn) kcat = ( subsystem_kcats.get(sub, standard_kcat) if sub diff --git a/src/geckopy/kcat_tuning/sensitivity_tuning.py b/src/geckopy/kcat_tuning/sensitivity_tuning.py index da903a5..a0a67f0 100644 --- a/src/geckopy/kcat_tuning/sensitivity_tuning.py +++ b/src/geckopy/kcat_tuning/sensitivity_tuning.py @@ -228,12 +228,10 @@ def sensitivity_tuning( model.objective = bio_rxn_id initial_sol = model.optimize() - if ( - initial_sol.objective_value is None - or np.isnan(initial_sol.objective_value) - ): + if initial_sol.status != "optimal": raise ValueError( - "FBA of the input model is infeasible. Reduce the protein " + "FBA of the input model is not optimal (status " + f"{initial_sol.status!r}). Reduce the protein " "pool constraint with set_prot_pool_size and/or check the " "exchange constraints." ) @@ -261,7 +259,13 @@ def sensitivity_tuning( while iteration <= max_iterations: sol = model.optimize() - growth = float(sol.objective_value or 0.0) + if sol.status != "optimal": + logger.warning( + "sensitivity_tuning: the model has no optimal solution with the " + "increased kcats (status %s); tuning stops.", sol.status, + ) + break + growth = float(sol.objective_value) if not np.isnan(last_growth) and growth <= last_growth + _GROWTH_STALL_TOL: logger.warning( diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 14a7398..c2f9f81 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -168,6 +168,8 @@ def flexibilize_enz_concs( Returns ------- FlexEnzResult + Empty when the LP of the input model is not optimal (for example + infeasible), as MATLAB's ``flexibilizeEnzConcs`` returns it. Raises ------ @@ -233,15 +235,14 @@ def flexibilize_enz_concs( unit_frequence = np.zeros(len(units), dtype=int) initial_sol = model.optimize() - if ( - initial_sol.objective_value is None - or np.isnan(initial_sol.objective_value) - ): - raise ValueError( - "FBA of the input model is infeasible, so there is no growth to " - "flexibilize towards. Check the enzyme-concentration and exchange " - "constraints before calling flexibilize_enz_concs." + if initial_sol.status != "optimal": + logger.warning( + "flexibilize_enz_concs: FBA of the input model is %s, so there is " + "no growth to flexibilize towards; no enzyme is flexibilized. " + "Check the enzyme-concentration and exchange constraints.", + initial_sol.status, ) + return FlexEnzResult() pred_growth = float(initial_sol.objective_value) flex_break = False @@ -316,7 +317,7 @@ def flexibilize_enz_concs( ) new_sol = model.optimize() - pred_growth = float(new_sol.objective_value or 0.0) + pred_growth = _optimal_growth(new_sol) if verbose: logger.info( "flexibilize_enz_concs: protein %r relaxed; growth = %g", @@ -334,6 +335,12 @@ def flexibilize_enz_concs( # Helpers # --------------------------------------------------------------------------- # +def _optimal_growth(sol) -> float: + """Objective value of an optimal solution, 0 for any other status; solvers + differ in the objective value they report for an infeasible LP.""" + return float(sol.objective_value) if sol.status == "optimal" else 0.0 + + def _try_relax_pool( model: "EcModel", bio_rxn, @@ -350,7 +357,7 @@ def _try_relax_pool( pool_rxn.upper_bound = 1000.0 bio_rxn.upper_bound = exp_growth sol = model.optimize() - candidate_growth = float(sol.objective_value or 0.0) + candidate_growth = _optimal_growth(sol) if candidate_growth - pred_growth <= _GROWTH_DELTA_THRESHOLD: if verbose: @@ -366,7 +373,7 @@ def _try_relax_pool( model.objective = POOL_EXCHANGE_ID model.objective_direction = "min" sol = model.optimize() - new_pool = float(sol.objective_value or 0.0) + new_pool = _optimal_growth(sol) pool_rxn.upper_bound = new_pool diff --git a/src/geckopy/limit_proteins/relax_proteomics_greedy.py b/src/geckopy/limit_proteins/relax_proteomics_greedy.py index 3c95829..424cf08 100644 --- a/src/geckopy/limit_proteins/relax_proteomics_greedy.py +++ b/src/geckopy/limit_proteins/relax_proteomics_greedy.py @@ -28,6 +28,8 @@ import numpy as np +from ..solver_support import constraint_shadow_prices + if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -113,15 +115,15 @@ def relax_proteomics_greedy( final_growth=float(growth), converged=False, ) - ranked = sorted( - candidates, - key=lambda u: abs(model.enzymes.get_by_id(u).shadow_price), - reverse=True, - ) - target = ranked[0] + ordered = sorted(candidates) + prices = dict(zip( + ordered, + constraint_shadow_prices(model, [model.enzymes.get_by_id(u).prot_metabolite_id for u in ordered]), + )) + target = max(ordered, key=lambda u: abs(prices[u])) enz = model.enzymes.get_by_id(target) original_conc = enz.concentration - sp = enz.shadow_price + sp = prices[target] enz.upper_bound = default_upper_bound relaxed[target] = original_conc diff --git a/src/geckopy/solver_support.py b/src/geckopy/solver_support.py new file mode 100644 index 0000000..01e4c7e --- /dev/null +++ b/src/geckopy/solver_support.py @@ -0,0 +1,649 @@ +"""Support for optlang's scipy (HiGHS) solver interface. + +The scipy interface of optlang differs from the other interfaces in ways that matter for an +ecModel, and this module handles them. + +Index map after a removal + ``scipy_interface.Problem`` keeps a map from variable name to column. After a variable is + removed it renumbers the columns with ``islice(variables, start, len(constraints))``, so + the last ``len(variables) - len(constraints)`` names keep their old column whenever the + problem has more variables than constraints. Primal values, and so every flux that cobra + reads, then come from the wrong columns or raise an ``IndexError``, and a later removal + of such a variable deletes the wrong column of the coefficient matrix. Variables are + removed lazily, by the next ``solver.update()``, which any ``add`` or ``optimize`` runs, so + the damage shows up in code that never removed anything itself. + + :func:`repair_scipy_problem` moves the problem to a subclass that renumbers correctly, so + every removal after it is right, and rebuilds the coefficient matrix from the constraint + expressions if the map was already stale. :class:`~geckopy.ec_model.ec_model.EcModel` + calls it whenever ``model.solver`` is read, so an EcModel is covered whatever removes the + variables: geckopy, cobra, or user code that calls ``model.remove_reactions``, + ``reaction.delete`` or ``model.remove_cons_vars``. geckopy builds ecModels on a different + interface and converts them with :func:`switch_to_scipy`, so no removal happens on an + unrepaired scipy problem inside geckopy. For a plain ``cobra.Model`` on the scipy solver + that you edit yourself, call ``repair_scipy_problem(model)`` before the first removal; a + call after the removals still restores the map and the matrix, but anything read in + between was wrong. + +Coefficient matrix + optlang stores the matrix densely and copies it on every added or removed variable or + constraint, so building or editing a genome-scale model (thousands of reactions) takes + hours. The subclass installed by :func:`repair_scipy_problem` keeps the coefficients in + dictionaries and hands HiGHS a sparse matrix, built once per structural change. + ``problem.A`` is then a ``scipy.sparse.csr_matrix`` (a dense array when there are no rows). + +Solution after a solve + ``cobra.Model.optimize`` sets the objective direction again after the solve, and the scipy + interface discards the solution on every direction assignment, also one that changes + nothing. ``Reaction.flux`` and ``Solver.primal_values`` then find no values after + ``model.optimize()``. The subclass discards the solution only when the direction changes. + + After a solve that finds no optimum (infeasible, unbounded) the scipy interface returns + ``None`` for every primal value and raises for the objective, so ``cobra.Model.optimize`` + fails with a ``TypeError`` instead of returning a solution with that status. The subclass + returns NaN for both, so the solution carries the status and NaN values. + +Duals + The scipy interface raises for ``Constraint.dual`` and ``Variable.dual`` but keeps the + whole ``linprog`` result in ``problem._solution``. :func:`constraint_shadow_prices` reads + the row marginals from it, so the shadow prices of a scipy model equal those of the glpk + and gurobi interfaces. The stored result is used only after it is shown to be an optimal + primal and dual solution of the problem as it is now, so a result left over from before + a bound or an objective changed is never read. +""" +from __future__ import annotations + +import contextlib +from itertools import islice +from typing import Iterable + +import numpy as np +from scipy import sparse + +try: + from optlang import scipy_interface as _scipy_interface + from optlang.util import SolverTolerances + from scipy.optimize import linprog +except ImportError: # pragma: no cover - optlang always ships the interface + _scipy_interface = None + +__all__ = [ + "repair_scipy_problem", + "switch_to_scipy", + "fast_solver_while_building", + "edit_on_fast_solver", + "is_scipy_interface", + "constraint_shadow_prices", +] + +# Relative tolerance of the optimality certificate, scaled by the size of the quantity checked. +_CERTIFICATE_TOL = 1e-6 + + +if _scipy_interface is not None: + + class _SparseProblem(_scipy_interface.Problem): + """``scipy_interface.Problem`` with a sparse coefficient matrix. + + Renumbers every variable column after a removal and keeps its solution when the + direction is set to the value it has. + """ + + def __init__(self): + super().__init__() + self._init_sparse() + + # Feasibility and optimality tolerance of HiGHS; None keeps its default (1e-7). + _feasibility = None + _optimality = None + + def optimize(self, method="highs", verbosity=False, tolerance=1e-9, **kwargs): + """``linprog`` on the problem, with the tolerances that were set on the model.""" + c = np.array([self.objective.get(name, 0) for name in self._variables], dtype=float) + if self.direction == "max": + c *= -1 + options = {"maxiter": 10000, "disp": verbosity} + bounds = list(self.bounds.values()) + + def solve(options): + return linprog(c, self.A, self.upper_bounds, bounds=bounds, method=method, options=options, **kwargs) + + tight = dict(options) + if self._feasibility is not None: + tight["primal_feasibility_tolerance"] = self._feasibility + if self._optimality is not None: + tight["dual_feasibility_tolerance"] = self._optimality + solution = solve(tight) + # HiGHS reports numerical difficulties (status 4) on some LPs at a tolerance + # below its default; those are solved again at the default. + if solution.status == 4 and tight != options: + solution = solve(options) + self._solution = solution + self._status = solution.status + if self._status == 0: + self._var_primals = solution.x + self._slacks = solution.slack + else: + self._var_primals = None + self._slacks = None + self._f = solution.fun + + @property + def status(self): + """The solver status; a solve that ended with numerical difficulties is ``numeric``.""" + if self._status is None: + raise RuntimeError("Problem has not been optimized yet") + return _scipy_interface.SCIPY_STATUS.get(self._status, _scipy_interface.interface.NUMERIC) + + def _init_sparse(self): + # _rows[row name] = {variable name: coefficient}; _cols[variable name] = {row names} + self._rows = {} + self._cols = {name: set() for name in self._variables} + self._matrix = None + self._A = None + + @property + def direction(self): + return self._direction + + @direction.setter + def direction(self, value): + if value != self._direction: + self._direction = value + self._reset_solution() + + @property + def A(self): + """The linear coefficient matrix: sparse, or a dense array without rows.""" + n, m = len(self._variables), len(self._constraints) + if m == 0: + return np.zeros((0, n)) + if self._matrix is None: + data, rows, cols = [], [], [] + columns = self._variables + for row_name, row_index in self._constraints.items(): + for name, coefficient in self._rows[row_name].items(): + data.append(coefficient) + rows.append(row_index) + cols.append(columns[name]) + self._matrix = sparse.csr_matrix((data, (rows, cols)), shape=(m, n), dtype=float) + return self._matrix + + def get_var_primal(self, name): + """The primal value; NaN after a solve that found no optimum, None before any solve.""" + if self._var_primals is None: + return None if self._status is None else float("nan") + return self._var_primals[self._get_var_index(name)] + + @property + def objective_value(self): + """The optimum; NaN after a solve that found none.""" + if self._f is None and self._status is not None: + return float("nan") + return super().objective_value + + def add_variable(self, name): + if name in self._variables: + raise ValueError("A variable named " + name + " already exists.") + self._variables[name] = len(self._variables) + self.bounds[name] = (0, None) + self._cols[name] = set() + self._matrix = None + self._reset_solution() + + def add_constraint(self, name, coefficients={}, ub=0): + if name in self._constraints: + raise ValueError("A constraint named " + name + " already exists.") + self._constraints[name] = len(self._constraints) + self.upper_bounds = np.append(self.upper_bounds, ub) + row = {v: c for v, c in coefficients.items() if c != 0 and v in self._variables} + self._rows[name] = row + for v in row: + self._cols[v].add(name) + self._matrix = None + self._reset_solution() + + def remove_variable(self, name): + index = self._get_var_index(name) + for row_name in self._cols.pop(name): + del self._rows[row_name][name] + del self.bounds[name] + del self._variables[name] + self._update_variable_indices(index) + self._matrix = None + self._reset_solution() + + def remove_constraint(self, name): + index = self._get_constraint_index(name) + self.upper_bounds = np.delete(self.upper_bounds, index) + for v in self._rows.pop(name): + self._cols[v].discard(name) + del self._constraints[name] + self._update_constraint_indices(index) + self._matrix = None + self._reset_solution() + + def change_var_name(self, name, value): + super().change_var_name(name, value) + if name != value: + for row_name in self._cols[name]: + self._rows[row_name][value] = self._rows[row_name].pop(name) + self._cols[value] = self._cols.pop(name) + + def change_constraint_name(self, name, value): + super().change_constraint_name(name, value) + if name != value: + self._rows[value] = self._rows.pop(name) + for v in self._rows[value]: + self._cols[v].discard(name) + self._cols[v].add(value) + + def _update_variable_indices(self, start=0): + i = start + for name in islice(self._variables, start, None): + self._variables[name] = i + i += 1 + + +if _scipy_interface is not None: + + class _SparseModel(_scipy_interface.Model): + """``scipy_interface.Model`` on a `_SparseProblem` that copies and pickles sparsely. + + optlang pickles a model as JSON and restores it by adding the variables and + constraints to a fresh dense problem, which takes tens of minutes for a genome-scale + model. This class restores into a `_SparseProblem`, so ``copy.deepcopy``, + ``cobra.Model.copy`` and the worker processes of a pool get their copy in seconds. + """ + + # optlang finds its interface module (the scipy one) through this attribute. + __module__ = _scipy_interface.__name__ + + def _initialize_problem(self): + self.problem = _SparseProblem() + + def _remove_constraints(self, constraints): + # optlang's scipy interface calls the generic removal inside its loop, once per + # constraint with the whole list, so removing two or more constraints at once + # (up to 350) raises ``LookupError: Constraint ... not in solver``. + for constraint in constraints: + if constraint.lb is not None: + self.problem.remove_constraint(constraint.lower_constraint_name) + if constraint.ub is not None: + self.problem.remove_constraint(constraint.upper_constraint_name) + _scipy_interface.interface.Model._remove_constraints(self, constraints) + + def _initialize_configuration(self): + super()._initialize_configuration() + _install_tolerances(self) + + def __reduce__(self): + problem = self.problem + return (_restore_sparse_model, (self.to_json(), problem._feasibility, problem._optimality)) + + def _restore_sparse_model(state, feasibility, optimality): + model = _SparseModel() + model._init_from_json(state) + model.problem._feasibility = feasibility + model.problem._optimality = optimality + return model + + def _install_tolerances(solver): + """Let ``solver.configuration.tolerances`` (and so ``cobra.Model.tolerance``) set the + feasibility and optimality tolerance of HiGHS, which optlang's scipy interface ignores.""" + + def setter(attribute): + return lambda value: setattr(solver.problem, attribute, float(value)) + + solver.configuration.tolerances = SolverTolerances( + { + "feasibility": (lambda: solver.problem._feasibility or 1e-7, setter("_feasibility")), + "optimality": (lambda: solver.problem._optimality or 1e-7, setter("_optimality")), + } + ) + + +def is_scipy_interface(interface) -> bool: + """True if ``interface`` (a solver name, an optlang interface module or a model's + ``problem`` attribute) is optlang's scipy interface.""" + if _scipy_interface is None: + return False + if isinstance(interface, str): + return interface.lower() == "scipy" + return interface is _scipy_interface + + +def _is_scipy_problem(problem) -> bool: + return _scipy_interface is not None and isinstance(problem, _scipy_interface.Problem) + + +def _map_is_stale(problem) -> bool: + return list(problem._variables.values()) != list(range(len(problem._variables))) + + +def _adopt_dense_matrix(problem) -> None: + """Fill the dictionaries of a `_SparseProblem` from the dense matrix it was swapped from.""" + dense = problem._A + names = list(problem._variables) + row_names = list(problem._constraints) + problem._init_sparse() + problem._rows = {name: {} for name in row_names} + for r, c in zip(*np.nonzero(dense)): + problem._rows[row_names[r]][names[c]] = float(dense[r, c]) + problem._cols[names[c]].add(row_names[r]) + + +def _rebuild_from_expressions(problem, solver) -> None: + """Fill the dictionaries of a `_SparseProblem` from the constraint expressions of ``solver``. + + The expressions are the source the matrix is derived from, so this is correct whatever + happened to the matrix itself. + """ + problem._init_sparse() + problem._rows = {name: {} for name in problem._constraints} + for i, name in enumerate(problem._variables): + problem._variables[name] = i + for constraint in solver.constraints: + coefficients = None + for row_name, sign in ( + (constraint.upper_constraint_name, 1.0), + (constraint.lower_constraint_name, -1.0), + ): + if row_name not in problem._rows: + continue + if coefficients is None: + coefficients = constraint.coefficient_dict() + row = {v: sign * c for v, c in coefficients.items() if c != 0 and v in problem._variables} + problem._rows[row_name] = row + for v in row: + problem._cols[v].add(row_name) + + +def repair_scipy_problem(model) -> bool: + """Make the scipy problem of ``model`` correct and fast to edit. + + Does nothing and returns False when the solver of ``model`` is not optlang's scipy + interface. Otherwise it moves the problem to a class that renumbers the variable columns + correctly after a removal and stores the coefficients sparsely (see the module docstring), + applies the pending changes, restores the column map and the coefficient matrix if a + removal had already left the map stale, and returns True. Once the problem is switched, a + call is a type check. + + Parameters + ---------- + model + A ``cobra.Model`` (or EcModel). Reading ``model.solver`` of an EcModel calls this + function already; call it yourself for a plain ``cobra.Model`` before you remove + variables (``remove_reactions``, ``Reaction.delete``, ``remove_cons_vars``). + """ + solver = getattr(model, "_solver", None) + problem = getattr(solver, "problem", None) + if not _is_scipy_problem(problem): + return False + if isinstance(problem, _SparseProblem): + return True + stale = _map_is_stale(problem) + problem.__class__ = _SparseProblem + solver.__class__ = _SparseModel + _install_tolerances(solver) + if stale: + problem._init_sparse() + problem._rows = {name: {} for name in problem._constraints} + else: + _adopt_dense_matrix(problem) + solver.update() + if stale: + _rebuild_from_expressions(problem, solver) + return True + + +def switch_to_scipy(model) -> None: + """Move ``model`` to the scipy solver, fast, with the sparse problem of this module. + + ``model.solver = "scipy"`` clones the problem through optlang, which grows a dense matrix + row by row. This builds the problem the way ``cobra.Model`` does when it is created, on + the sparse class, so a genome-scale model switches in seconds. The objective and its + direction are kept. A model whose solver holds anything besides one forward and one + reverse variable per reaction and one constraint per metabolite (an extra constraint or + variable of your own) is switched with optlang's clone instead and then repaired. + """ + if repair_scipy_problem(model): + return + solver = model.solver + plain = len(solver.constraints) == len(model.metabolites) and len(solver.variables) == 2 * len( + model.reactions + ) + if _scipy_interface is None or not plain: + from cobra import Model + + Model.solver.fset(model, "scipy") + repair_scipy_problem(model) + return + from cobra.util.solver import linear_reaction_coefficients + from optlang.symbolics import Zero + + coefficients = linear_reaction_coefficients(model) + direction = model.objective.direction + new = _SparseModel() + new.objective = _scipy_interface.Objective(Zero) + model._solver = new + _populate_sparse(model, new) + model.objective = coefficients + model.objective.direction = direction + if getattr(model, "_tolerance", None) is not None: + try: + model.tolerance = model._tolerance + except Exception: # the scipy interface has no tolerance settings to receive it + pass + + +def _populate_sparse(model, solver) -> None: + """Give ``solver`` one forward and one reverse variable per reaction of ``model`` and one + constraint per metabolite, as ``cobra.Model._populate_solver`` does. + + ``_populate_solver`` sets the coefficients of a constraint one at a time, and the scipy + interface removes and adds the rows of a constraint for each call, which takes hours for a + genome-scale model. Here every constraint is created with its whole expression and added + once. + """ + from optlang.symbolics import Zero, add + + problem = solver.interface + variables = {} + for reaction in model.reactions: + variables[reaction.id] = (problem.Variable(reaction.id), problem.Variable(reaction.reverse_id)) + solver.add([v for pair in variables.values() for v in pair]) + terms = {metabolite.id: [] for metabolite in model.metabolites} + for reaction in model.reactions: + forward, reverse = variables[reaction.id] + for metabolite, coefficient in reaction.metabolites.items(): + terms[metabolite.id].append(coefficient * forward) + terms[metabolite.id].append(-coefficient * reverse) + solver.add( + [ + problem.Constraint(add(*expressions) if expressions else Zero, name=name, lb=0, ub=0) + for name, expressions in terms.items() + ] + ) + for reaction in model.reactions: + reaction.update_variable_bounds() + + +@contextlib.contextmanager +def edit_on_fast_solver(model): + """Run structural edits of ``model`` on glpk when it is on the scipy solver. + + Removing reactions and metabolites, or rewriting the coefficients of thousands of + reactions, is quick on glpk and slow on the scipy interface, where optlang rewrites the + constraint expressions (and the rows) for each edit. Inside the block a scipy model is on + glpk; on exit it is converted back with :func:`switch_to_scipy`. The block changes nothing + for a model on another solver, for a model inside a ``with model:`` context (its + rollbacks hold the solver objects) and when glpk is unavailable. The optlang solver of + the model is a different object inside the block and after it. + """ + from cobra import Model + + solver = getattr(model, "_solver", None) + if not _is_scipy_problem(getattr(solver, "problem", None)) or getattr(model, "_contexts", None): + yield False + return + try: + Model.solver.fset(model, "glpk") + except Exception: + yield False + return + try: + yield True + finally: + switch_to_scipy(model) + + +@contextlib.contextmanager +def fast_solver_while_building(): + """Make models created inside the block on a solver that builds large models quickly. + + When cobra's default solver is scipy the block runs with glpk as the default, which + builds a genome-scale model in seconds; the default is restored on exit. Models created + inside are on glpk; convert them with :func:`switch_to_scipy`. With any other default + solver, or when glpk is unavailable, the block changes nothing. + """ + import cobra + + config = cobra.Configuration() + previous = config.solver + switched = False + if getattr(previous, "__name__", str(previous)).endswith("scipy_interface"): + try: + config.solver = "glpk" + switched = True + except Exception: + switched = False + try: + yield switched + finally: + if switched: + config.solver = previous + + +def _objective_vector(problem) -> np.ndarray: + """Cost vector ``linprog`` minimises: the objective, negated for a maximisation.""" + cost = np.array([problem.objective.get(name, 0.0) for name in problem._variables], dtype=float) + if problem.direction == "max": + cost = -cost + return cost + + +def _certified_result(problem): + """The stored ``linprog`` result if it is an optimum of ``problem`` as it is now, else None. + + Checks primal feasibility against the current rows and bounds, stationarity of the cost + vector against the stored marginals, and equality of the primal and dual objective. A + result that passes is an optimal primal-dual pair of the current problem, whatever solve + produced it. + """ + result = getattr(problem, "_solution", None) + try: + x = np.asarray(result.x, dtype=float) + y = np.asarray(result.ineqlin.marginals, dtype=float) + z_low = np.asarray(result.lower.marginals, dtype=float) + z_up = np.asarray(result.upper.marginals, dtype=float) + fun = float(result.fun) + matrix = sparse.csr_matrix(problem.A, dtype=float) + rhs = np.asarray(problem.upper_bounds, dtype=float) + bounds = list(problem.bounds.values()) + except (AttributeError, TypeError, ValueError): + return None + n, m = len(problem._variables), len(problem._constraints) + if matrix.shape != (m, n) or not (x.shape == z_low.shape == z_up.shape == (n,)) or y.shape != (m,): + return None + if len(bounds) != n or not np.all(np.isfinite(x)): + return None + low = np.array([-np.inf if b[0] is None else b[0] for b in bounds], dtype=float) + up = np.array([np.inf if b[1] is None else b[1] for b in bounds], dtype=float) + cost = _objective_vector(problem) + tol = _CERTIFICATE_TOL + + def close(a, b, scale): + return abs(a - b) <= tol * (1.0 + scale) + + if m and np.any(matrix @ x - rhs > tol * (1.0 + np.abs(rhs))): + return None + if np.any(x < low - tol * (1.0 + np.abs(low))) or np.any(x > up + tol * (1.0 + np.abs(up))): + return None + reduced = cost - (matrix.T @ y if m else 0.0) + if np.max(np.abs(reduced - (z_low + z_up)), initial=0.0) > tol * (1.0 + np.max(np.abs(cost), initial=0.0)): + return None + primal = float(cost @ x) + dual = float(rhs @ y) if m else 0.0 + dual += float(np.sum(low[np.isfinite(low)] * z_low[np.isfinite(low)])) + dual += float(np.sum(up[np.isfinite(up)] * z_up[np.isfinite(up)])) + if not (close(primal, fun, abs(fun)) and close(primal, dual, abs(primal))): + return None + return result + + +def constraint_shadow_prices(model, constraint_ids: Iterable[str]) -> np.ndarray: + """Shadow prices of the named constraints, on every solver, as one array. + + The shadow price of a constraint is the change of the optimum, in the direction of the + model objective, per unit rise of its right-hand side: for a metabolite balance + ``S v = 0`` it is the optimum gained per unit of the metabolite supplied. It is the value + of ``cobra.Metabolite.shadow_price`` for the glpk, gurobi and cplex interfaces, and for the + scipy interface (which raises) it is computed from the HiGHS row marginals of the last + solve. + + The model must have been solved. On the scipy interface the stored result is used only if + it is an optimum of the problem as it is now; otherwise the model is solved once more. If + that solve is not optimal, the function raises. + + Parameters + ---------- + model + An EcModel or ``cobra.Model`` that has been optimized. + constraint_ids + IDs of constraints of the model, for example ``prot_``. + + Returns + ------- + numpy.ndarray + One shadow price per ID, in the order given. + + Raises + ------ + RuntimeError + If the scipy problem has no optimal solution to read the marginals from. + KeyError + If an ID is not a constraint of the model. + """ + ids = list(constraint_ids) + constraints = model.constraints + repair_scipy_problem(model) + problem = getattr(model.solver, "problem", None) + if _scipy_interface is None or not isinstance(problem, _scipy_interface.Problem): + return np.array([float(constraints[cid].dual) for cid in ids], dtype=float) + + result = _certified_result(problem) + if result is None: + model.slim_optimize() + result = _certified_result(problem) + if result is None: + raise RuntimeError( + "shadow prices need an optimal solution of the scipy problem; the model has none " + f"(solver status {model.solver.status!r})" + ) + marginals = np.asarray(result.ineqlin.marginals, dtype=float) + rows = problem._constraints + # The rows hold `A x <= ub` and `-A x <= -lb`. The marginal of a row is d(minimised value)/d(rhs + # of that row), so d/d(ub) is the upper row's marginal and d/d(lb) is minus the lower row's. + # A maximisation minimises the negated objective. + sign = -1.0 if problem.direction == "max" else 1.0 + prices = np.zeros(len(ids), dtype=float) + for k, cid in enumerate(ids): + constraint = constraints[cid] + value = 0.0 + if constraint.ub is not None: + value += marginals[rows[constraint.upper_constraint_name]] + if constraint.lb is not None: + value -= marginals[rows[constraint.lower_constraint_name]] + prices[k] = sign * value + return prices diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 2ae3d92..39afee3 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -22,6 +22,8 @@ import pandas as pd +from ..solver_support import constraint_shadow_prices + if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -63,12 +65,14 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: f"Solver status {sol.status!r}; cannot rank bottlenecks" ) + enzymes = list(model.enzymes) + prices = constraint_shadow_prices(model, [enz.prot_metabolite_id for enz in enzymes]) rows = [] - for enz in model.enzymes: + for enz, price in zip(enzymes, prices): rows.append({ "uniprot": enz.id, "gene": enz.gene, - "shadow_price": enz.shadow_price, + "shadow_price": float(price), "flux": enz.flux, "cap_usage": enz.cap_usage, "upper_bound": enz.upper_bound, diff --git a/src/geckopy/utilities/ec_fseof.py b/src/geckopy/utilities/ec_fseof.py index 7e51811..c41b79e 100644 --- a/src/geckopy/utilities/ec_fseof.py +++ b/src/geckopy/utilities/ec_fseof.py @@ -13,6 +13,11 @@ - restricts targets (and the reported scan) to gene-associated reactions whose GPR doesn't reference the "standard" placeholder pseudogene, matching ``ecFSEOF.m``'s candidate-set restriction; +- fixes the target flux at each level (lower and upper bound), as + ``ecFSEOF.m`` does; +- orders ``targets`` as ``ecFSEOF.m`` returns them: reaction targets first, + then transport targets, each group sorted by score, with a ``transport`` + boolean column that marks the group; - filters ``usage_prot_*`` rows out of the returned ``scan`` / ``targets`` DataFrames — those reactions are protein-pool plumbing, never engineering targets. @@ -32,9 +37,9 @@ ``result.targets.to_csv(...)`` themselves. Per-gene essentiality is not computed here; use :func:`cobra.flux_analysis.single_gene_deletion`, which works directly on ecModels since GPRs already gate the catalysed -reactions. Targets come back as one combined table; split it into -reaction vs. transport subsets yourself via the ``subsystem`` column -or by inspecting each reaction's metabolites. +reactions. Targets come back as one table, reaction targets before transport +targets (``ecFSEOF.m``'s ``rxnTargets`` and ``transportTargets``); the +``transport`` column tells them apart. """ from __future__ import annotations @@ -107,7 +112,14 @@ def ec_fseof( ------- FSEOFResult raven-toolbox's FSEOF output dataclass, with ``usage_prot_*`` - rows filtered out of both ``scan`` and ``targets``. See + 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 + 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 + top reaction targets, then transport targets. See :class:`raven_toolbox.analysis.fseof.FSEOFResult` for the full field list (``scan``, ``enforced``, ``targets``, plus the ``amplification`` / ``knockout`` / ``gene_targets`` properties). @@ -131,9 +143,11 @@ def ec_fseof( correlation_threshold=correlation_threshold, flux_eps=flux_eps, min_target=min_target, + equality=True, ) result = _drop_usage_prot(result) - return _drop_ungated_reactions(result, model) + result = _drop_ungated_reactions(result, model) + return _group_transport_targets(result, model) def _check_cs_consistency( @@ -200,6 +214,26 @@ def _drop_ungated_reactions(result: FSEOFResult, model: "EcModel") -> FSEOFResul return FSEOFResult(scan=scan, enforced=result.enforced, targets=targets) +def _group_transport_targets(result: FSEOFResult, model: "EcModel") -> FSEOFResult: + """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``. + """ + 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) + return FSEOFResult(scan=result.scan, enforced=result.enforced, targets=targets) + + +def _is_transport(reaction) -> bool: + names = [m.name for m in reaction.metabolites if not (m.name or "").startswith("prot_")] + return len(names) != len(set(names)) + + def _drop_usage_prot(result: FSEOFResult) -> FSEOFResult: """Strip ``usage_prot_*`` rows from raven's ``scan`` and ``targets``. diff --git a/src/geckopy/utilities/get_subset_ec_model.py b/src/geckopy/utilities/get_subset_ec_model.py index 4f602eb..3e2c019 100644 --- a/src/geckopy/utilities/get_subset_ec_model.py +++ b/src/geckopy/utilities/get_subset_ec_model.py @@ -5,12 +5,15 @@ """ from __future__ import annotations +import contextlib import copy import logging from typing import TYPE_CHECKING import numpy as np +from cobra.manipulation import remove_genes +from ..solver_support import edit_on_fast_solver from ..ec_model.constants import ( POOL_EXCHANGE_ID, PROT_PREFIX, @@ -51,16 +54,20 @@ def get_subset_ec_model( 2. Compute genes to remove (`big.genes` minus `small_gem.genes`, preserving the ``standard`` pseudo-gene from ``assign_standard_kcat``). - 3. Remove reactions blocked by that gene removal (cobra's - ``GPR.eval(knockouts=...)``). - 4. Remove the genes from `model.genes`. - 5. Trim the corresponding columns of + 3. Remove the genes with ``cobra.manipulation.remove_genes``: the + reactions the removal blocks go, and the genes leave the rule of + every reaction that survives (an isozyme remains). + 4. Trim the corresponding columns of ``model.ec.genes/enzymes/mw/sequence/concs/rxn_enz_mat``. - 6. Remove reactions whose canonicalised id is not in + 5. Remove reactions whose canonicalised id is not in ``small_gem.reactions``, except ``usage_prot_*`` and ``prot_pool_exchange`` which are always preserved. + 6. Remove the metabolites and genes that no reaction uses, as + MATLAB's ``removeGenes`` and ``removeReactions`` do. 7. Trim the corresponding rows of - ``model.ec.rxns/kcat/source/notes/eccodes/rxn_enz_mat``. + ``model.ec.rxns/kcat/source/notes/eccodes/rxn_enz_mat``. In a + gecko-light model the ec reaction ids carry a ``001_`` prefix that + is not part of the reaction id. Ported from GECKO MATLAB: src/geckomat/utilities/getSubsetEcModel.m. @@ -124,6 +131,11 @@ def get_subset_ec_model( # deep-copy it explicitly to avoid mutating the input model. small = big_ec_model.copy() small.ec = copy.deepcopy(big_ec_model.ec) + # The removals below run on glpk when the model is on scipy and the model is + # converted back at the end (the scipy interface rewrites every constraint + # expression for a removed variable). + solver_edit = contextlib.ExitStack() + solver_edit.enter_context(edit_on_fast_solver(small)) # --- 4. Compute genes to remove --- small_gem_gene_ids = {g.id for g in small_gem.genes} @@ -131,26 +143,13 @@ def get_subset_ec_model( genes_to_remove = ( big_gene_ids - small_gem_gene_ids - {_STANDARD_GENE} ) - if not genes_to_remove and not _need_rxn_trim(small, small_gem): - return small # nothing to do - # --- 5. Remove reactions blocked by the gene removal --- - rxns_to_remove_for_genes = [] - for rxn in list(small.reactions): - if not rxn.gene_reaction_rule: - continue - rxn_gene_ids = {g.id for g in rxn.genes} - if not (rxn_gene_ids & genes_to_remove): - continue - if not rxn.gpr.eval(knockouts=genes_to_remove): - rxns_to_remove_for_genes.append(rxn) - if rxns_to_remove_for_genes: - small.remove_reactions(rxns_to_remove_for_genes, remove_orphans=False) - - # --- 6. Remove the genes from model.genes --- - for gene_id in list(genes_to_remove): - if gene_id in {g.id for g in small.genes}: - small.genes.remove(small.genes.get_by_id(gene_id)) + # --- 5. Remove the genes, and the reactions they block --- + # A reaction that survives (an isozyme remains) loses the gene from its + # rule, so no kept reaction refers to a gene that is no longer in + # model.genes. A complex that loses a subunit is dropped whole. + if genes_to_remove: + remove_genes(small, sorted(genes_to_remove), remove_reactions=True) # --- 7. Trim ec per-enzyme fields --- if genes_to_remove: @@ -213,10 +212,28 @@ def get_subset_ec_model( rxns_to_remove_for_subset, remove_orphans=True, ) + # --- 8b. Remove the metabolites and genes no reaction uses --- + # getSubsetEcModel's removeGenes and removeReactions calls drop every + # metabolite without a reaction and every gene without a reaction, also + # those that were unused in the big model (such as the ``standard`` + # pseudo-gene of a light model). + unused_mets = [m for m in small.metabolites if not m.reactions] + if unused_mets: + small.remove_metabolites(unused_mets, destructive=False) + unused_genes = [g for g in small.genes if not g.reactions] + if unused_genes: + remove_genes(small, unused_genes, remove_reactions=False) + + solver_edit.close() + # --- 9. Trim ec per-rxn fields --- cobra_rxn_ids = {r.id for r in small.reactions} keep_rxn_mask = np.array( - [rid in cobra_rxn_ids for rid in small.ec.rxns], dtype=bool, + [ + (rid[4:] if small.ec.gecko_light else rid) in cobra_rxn_ids + for rid in small.ec.rxns + ], + dtype=bool, ) if not keep_rxn_mask.all(): small.ec.rxns = [ @@ -250,15 +267,3 @@ def get_subset_ec_model( def _canonical(rxn_id: str) -> str: """Strip ``_REV`` and ``_EXP_`` suffixes.""" return canonicalize_rxn_id(rxn_id)[0] - - -def _need_rxn_trim(big: "EcModel", small_gem: "cobra.Model") -> bool: - """True if there's at least one reaction in big that's not in - small_gem and not part of the protein machinery.""" - small_gem_rxn_ids = {r.id for r in small_gem.reactions} - for rxn in big.reactions: - if rxn.id.startswith(USAGE_PREFIX) or rxn.id == POOL_EXCHANGE_ID: - continue - if _canonical(rxn.id) not in small_gem_rxn_ids: - return True - return False diff --git a/src/geckopy/utilities/load_conventional_gem.py b/src/geckopy/utilities/load_conventional_gem.py index 420e8d0..812d8b7 100644 --- a/src/geckopy/utilities/load_conventional_gem.py +++ b/src/geckopy/utilities/load_conventional_gem.py @@ -23,6 +23,8 @@ import cobra from raven_toolbox.io import read_yaml_model +from ..solver_support import fast_solver_while_building, switch_to_scipy + if TYPE_CHECKING: from ..adapter import ModelAdapter @@ -58,10 +60,17 @@ def load_conventional_gem(adapter: "ModelAdapter") -> cobra.Model: if not path.is_file(): raise FileNotFoundError(f"conv_gem file not found: {path}") suffix = path.suffix.lower() - if suffix in _YAML_SUFFIXES: - return read_yaml_model(path) - if suffix in _JSON_SUFFIXES: - return cobra.io.load_json_model(str(path)) - if suffix in _MAT_SUFFIXES: - return cobra.io.load_matlab_model(str(path)) - return cobra.io.read_sbml_model(str(path)) + # With the scipy solver as cobra's default the model is read on glpk and converted: optlang + # grows the scipy matrix densely, which takes hours for a genome-scale model. + with fast_solver_while_building() as on_fast_solver: + if suffix in _YAML_SUFFIXES: + model = read_yaml_model(path) + elif suffix in _JSON_SUFFIXES: + model = cobra.io.load_json_model(str(path)) + elif suffix in _MAT_SUFFIXES: + model = cobra.io.load_matlab_model(str(path)) + else: + model = cobra.io.read_sbml_model(str(path)) + if on_fast_solver: + switch_to_scipy(model) + return model diff --git a/src/geckopy/utilities/load_ec_model.py b/src/geckopy/utilities/load_ec_model.py index 8853935..18cc53a 100644 --- a/src/geckopy/utilities/load_ec_model.py +++ b/src/geckopy/utilities/load_ec_model.py @@ -20,7 +20,16 @@ clear message if the YAML lacks the ec sections, and attach the adapter to ``model.adapter`` so downstream geckopy code finds it. -4. **ecModel format** — a file whose protein usage and pool reactions +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 + 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 + removed and the quotes inside stay. + +5. **ecModel format**: a file whose protein usage and pool reactions raven-toolbox turned forward from the GECKO 3 direction gives a warning that names :func:`~geckopy.upgrade_project`; an ecModel format (``gecko_format`` in the file's ``metaData``, see @@ -34,7 +43,9 @@ """ from __future__ import annotations +import json import re +import tempfile import threading import warnings from pathlib import Path @@ -43,6 +54,7 @@ from raven_toolbox.io import EcData, read_yaml_model from ..ec_model.ec_model import EcModel +from ..solver_support import fast_solver_while_building, switch_to_scipy from ..upgrade import upgrade_ec_model if TYPE_CHECKING: @@ -108,7 +120,11 @@ def count_turned(message, category, filename, lineno, file=None, line=None): with _HOOK_LOCK, warnings.catch_warnings(): warnings.filterwarnings("always", message=r"ecModel uses the older reverse-direction convention") warnings.showwarning = count_turned - cobra_model = read_yaml_model(path) + # With the scipy solver as cobra's default the file is read on glpk and the + # ecModel is converted: optlang grows the scipy matrix densely, which takes + # hours for a genome-scale model. + with fast_solver_while_building() as on_fast_solver: + cobra_model = _read_yaml(path) meta = (cobra_model.notes or {}).get("metaData") or {} if _format_number(meta.get("gecko_format")) > 4: raise ValueError( @@ -127,6 +143,8 @@ def count_turned(message, category, filename, lineno, file=None, line=None): cobra_model, adapter=adapter, gecko_light=ec.gecko_light, ) 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. @@ -141,6 +159,70 @@ def count_turned(message, category, filename, lineno, file=None, line=None): return ec_model +_EMPTY_OMAP = re.compile( + # A key whose value is the tag `!!omap` and nothing else, when the next + # non-blank line is not indented deeper than the line of the key. + r"^( *)((?:- )?[^\n:]*: )!!omap[ \t]*$(?!\n(?:[ \t]*\n)*\1 +\S)", + re.MULTILINE, +) +_QUOTED_VALUE = re.compile( + r'^( *-(?: +[^\s":][^":\n]*:)? +)"(.*)"[ \t]*$', + re.MULTILINE, +) + + +def _read_yaml(path: Path): + """``read_yaml_model``; a file it rejects is read again after + :func:`_repair_yaml_text`, and the first error is raised when that + changes nothing.""" + try: + return read_yaml_model(path) + except Exception: + text = path.read_text(encoding="utf-8") + repaired = _repair_yaml_text(text) + if repaired == text: + raise + with tempfile.TemporaryDirectory() as folder: + fixed = Path(folder) / "ecModel.yml" + fixed.write_text(repaired, encoding="utf-8") + return read_yaml_model(fixed) + + +def _repair_yaml_text(text: str) -> str: + """Write the two forms that MATLAB reads and YAML does not as YAML. + + ``key: !!omap`` without entries becomes ``key: !!omap []``. A value in + double quotes that holds a double quote which is not escaped is read as + RAVEN's ``readYAMLmodel`` reads it: without its outer pair of quotes and + with an escaped backslash or quote as the character itself, and written + back as a YAML string. + """ + text = _EMPTY_OMAP.sub(r"\1\2!!omap []", text) + + def requote(match: re.Match) -> str: + inner = match.group(2) + if not _has_unescaped_quote(inner): + return match.group(0) + value = inner.replace("\\\\", "\0").replace('\\"', '"').replace("\0", "\\") + return match.group(1) + json.dumps(value, ensure_ascii=False) + + return _QUOTED_VALUE.sub(requote, text) + + +def _has_unescaped_quote(inner: str) -> bool: + """Whether ``inner``, the text between the outer quotes of a value, has a + double quote that no backslash escapes.""" + i = 0 + while i < len(inner): + if inner[i] == "\\": + i += 2 + continue + if inner[i] == '"': + return True + i += 1 + return False + + def _format_number(value) -> float: """``gecko_format`` as a number; a value that is not a number counts as no format, as GECKO's ``loadEcModel`` reads it with ``str2double``.""" diff --git a/src/geckopy/utilities/pfba_enzymes.py b/src/geckopy/utilities/pfba_enzymes.py index b4c3197..39581f7 100644 --- a/src/geckopy/utilities/pfba_enzymes.py +++ b/src/geckopy/utilities/pfba_enzymes.py @@ -33,12 +33,17 @@ _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 + def pfba_enzymes( model: "EcModel", objective: Union[cobra.Reaction, str, None] = None, *, - fraction_of_optimum: float = 1.0, + fraction_of_optimum: float = _PFBA_FRACTION, ) -> cobra.Solution: """Parsimonious FBA minimising the L1 norm of ``usage_prot_*`` fluxes. @@ -64,9 +69,10 @@ 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. ``1.0`` (default) - requires the biomass-optimal solution; a lower value trades - some objective value for potentially lower enzyme usage. + 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. Returns ------- diff --git a/tests/test_adapter.py b/tests/test_adapter.py index f7244e2..13d4537 100644 --- a/tests/test_adapter.py +++ b/tests/test_adapter.py @@ -251,3 +251,19 @@ def test_adapter_parameter_table_matches_schema(): expected = tomllib.loads(f"x = {default_text}")["x"] default = field.get_default(call_default_factory=True) assert default == expected, key + + +def test_uniprot_table_lookup_reads_quoted_cells_and_first_duplicate(tmp_path): + folder = _write_adapter_folder(tmp_path, MINIMAL_TOML) + data_dir = folder / "data" + data_dir.mkdir() + (data_dir / "uniprotConversion.tsv").write_text( + "gene_id\tuniprot_id\n" + '"YAL001C"\t"P12345"\n' + "YAL001C\tP00000\n" + 'YAL002W\t"P67890"\r\n' + ) + + adapter = ModelAdapter.from_folder(folder) + result = adapter.get_uniprot_ids_from_table(["YAL001C", "YAL002W"]) + assert result == ["P12345", "P67890"] diff --git a/tests/test_ec_fseof.py b/tests/test_ec_fseof.py index b22268b..79048fa 100644 --- a/tests/test_ec_fseof.py +++ b/tests/test_ec_fseof.py @@ -359,3 +359,59 @@ def test_knockout_slice_works(tmp_path): ko = result.knockout # raven calls this slice "knockout" but it includes knockdown assert isinstance(ko, pd.DataFrame) assert "R_BIO" in list(ko["reaction"]) + + +# --------------------------------------------------------------------------- # +# ecFSEOF.m's target order and fixed target flux +# --------------------------------------------------------------------------- # + +def test_targets_list_reaction_targets_before_transport_targets(): + """A reaction whose metabolite names repeat (not counting prot_*) is a transport target; + the score order is kept inside each group.""" + import importlib + + import pandas as pd + + module = importlib.import_module("geckopy.utilities.ec_fseof") + model = cobra.Model("m") + mets = { + "a_c": cobra.Metabolite("a_c", name="a", compartment="c"), + "a_e": cobra.Metabolite("a_e", name="a", compartment="e"), + "b_c": cobra.Metabolite("b_c", name="b", compartment="c"), + "p1": cobra.Metabolite("prot_E1", name="prot_E1", compartment="c"), + "p2": cobra.Metabolite("prot_E2", name="prot_E2", compartment="c"), + } + transport = cobra.Reaction("TRANSPORT") + transport.add_metabolites({mets["a_e"]: -1.0, mets["a_c"]: 1.0, mets["p1"]: -0.1}) + first = cobra.Reaction("FIRST") + first.add_metabolites({mets["a_c"]: -1.0, mets["b_c"]: 1.0, mets["p2"]: -0.1}) + second = cobra.Reaction("SECOND") + second.add_metabolites({mets["b_c"]: -1.0, mets["p1"]: -0.1, mets["p2"]: -0.1}) + model.add_reactions([transport, first, second]) + result = FSEOFResult( + scan=pd.DataFrame(), + enforced=[1.0], + targets=pd.DataFrame({"reaction": ["TRANSPORT", "FIRST", "SECOND"], "score": [3.0, 2.0, 1.0]}), + ) + + grouped = module._group_transport_targets(result, model).targets + + assert list(grouped["reaction"]) == ["FIRST", "SECOND", "TRANSPORT"] + assert list(grouped["transport"]) == [False, False, True] + + +def test_ec_fseof_fixes_the_target_flux_at_each_level(tmp_path, monkeypatch): + import importlib + + module = importlib.import_module("geckopy.utilities.ec_fseof") + seen = {} + real = module.fseof + + def spy(*args, **kwargs): + seen.update(kwargs) + return real(*args, **kwargs) + + monkeypatch.setattr(module, "fseof", spy) + 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 diff --git a/tests/test_flexibilize_enz_concs.py b/tests/test_flexibilize_enz_concs.py index 457fbb1..935bf58 100644 --- a/tests/test_flexibilize_enz_concs.py +++ b/tests/test_flexibilize_enz_concs.py @@ -241,3 +241,45 @@ def test_flex_concs_greater_than_old_for_returned_enzymes(tmp_path): result = flexibilize_enz_concs(model) for old, flex in zip(result.old_concs, result.flex_concs): assert flex >= old + + +# --------------------------------------------------------------------------- # +# Infeasible start: the same on every solver +# --------------------------------------------------------------------------- # + +def _usable_solvers(): + solvers = [] + for name in ("glpk", "glpk_exact", "gurobi", "scipy"): + try: + cobra.Model("probe").solver = name + except Exception: + continue + solvers.append(name) + return solvers + + +@pytest.fixture(params=_usable_solvers()) +def solver(request): + """Solver name; a solver that imports but is not licensed is skipped.""" + model = cobra.Model("probe") + try: + model.solver = request.param + probe = cobra.Reaction("probe_rxn") + probe.bounds = (0, 1) + model.add_reactions([probe]) + model.optimize() + except Exception as exc: + pytest.skip(f"{request.param} unusable: {exc}") + return request.param + + +def test_infeasible_start_returns_empty_result_on_every_solver(tmp_path, solver): + adapter = _adapter(tmp_path, gr_exp=500.0, bio_rxn="biomass") + model = _build_toy(adapter, enzyme_ub=1.0, measured_conc=1.0) + model.solver = solver + model.reactions.get_by_id("R_AB").lower_bound = 50.0 + assert model.optimize().status == "infeasible" + result = flexibilize_enz_concs(model) + assert result.uniprot_ids == [] + assert model.reactions.get_by_id("usage_prot_E").upper_bound == 1.0 + assert model.reactions.get_by_id("R_AB").lower_bound == 50.0 diff --git a/tests/test_get_standard_kcat.py b/tests/test_get_standard_kcat.py index b50da76..8b10c63 100644 --- a/tests/test_get_standard_kcat.py +++ b/tests/test_get_standard_kcat.py @@ -596,3 +596,63 @@ def test_fill_zero_kcat_false_leaves_zeros_alone(tmp_path): assert model.ec.kcat[r2_idx] == 0.0 + + +# --------------------------------------------------------------------------- # +# Subsystem names compare without regard to letter case +# --------------------------------------------------------------------------- # + +def _model_with_subsystems(tmp_path, subsystems: dict[str, str], missing_subsystem): + """12 reactions with a kcat (6 with 2.0 and 6 with 4.0) named in + ``subsystems`` (by index), 5 reactions of subsystem ``Other`` with 100.0, + and ``missing_r`` without genes in ``missing_subsystem``. The global + standard kcat is 4.0.""" + adapter = _adapter(tmp_path) + n = 17 + rxn_specs = [ + (f"r{i}", [(f"m{i}_in", -1.0, f"m{i}_in", "c"), (f"m{i}_out", 1.0, f"m{i}_out", "c")], f"g{i}") + for i in range(n) + ] + rxn_specs.append( + ("missing_r", [("X_in", -1.0, "X_in", "c"), ("X_out", 1.0, "X_out", "c")], None) + ) + kcats = [2.0] * 6 + [4.0] * 6 + [100.0] * 5 + model = _build_model_with_pool( + adapter, rxn_specs, + ec_rxns=[f"r{i}" for i in range(n)], + ec_kcats=kcats, + ec_genes=[f"g{i}" for i in range(n)], + ec_enzymes=[f"g{i}" for i in range(n)], + ec_mws=[100.0] * n, + rxn_to_ec_genes=[[i] for i in range(n)], + ) + for i in range(n): + model.reactions.get_by_id(f"r{i}").subsystem = ( + subsystems.get(i, "Other") if i < 12 else "Other" + ) + model.reactions.get_by_id("missing_r").subsystem = missing_subsystem + return model + + +def test_subsystem_names_differing_in_case_are_one_subsystem(tmp_path): + """"Alpha" and "alpha" have 6 reactions each, below the threshold of 10 on + their own; together they reach it and give their joint median, 3.0.""" + subsystems = {i: "Alpha" if i < 6 else "alpha" for i in range(12)} + model = _model_with_subsystems(tmp_path, subsystems, "ALPHA") + assign_standard_kcat(model, _uniprot([100.0]), threshold=10) + assert model.ec.kcat[model.ec.rxns.index("missing_r")] == pytest.approx(3.0) + + +def test_empty_subsystem_gets_the_standard_kcat(tmp_path): + subsystems = {i: "Alpha" for i in range(12)} + model = _model_with_subsystems(tmp_path, subsystems, "") + assign_standard_kcat(model, _uniprot([100.0]), threshold=10) + assert model.ec.kcat[model.ec.rxns.index("missing_r")] == pytest.approx(4.0) + + +def test_empty_first_subsystem_of_a_list_gets_the_standard_kcat(tmp_path): + subsystems = {i: "Alpha" for i in range(12)} + model = _model_with_subsystems(tmp_path, subsystems, "") + model.reactions.get_by_id("missing_r").subsystem = ["", "Alpha"] + assign_standard_kcat(model, _uniprot([100.0]), threshold=10) + assert model.ec.kcat[model.ec.rxns.index("missing_r")] == pytest.approx(4.0) diff --git a/tests/test_get_subset_ec_model.py b/tests/test_get_subset_ec_model.py index 5914598..7b7ab1b 100644 --- a/tests/test_get_subset_ec_model.py +++ b/tests/test_get_subset_ec_model.py @@ -421,3 +421,147 @@ def test_default_usage_bounds_emit_no_warning(tmp_path, caplog): "protein-concentration constraints" in rec.message for rec in caplog.records ) + + +# --------------------------------------------------------------------------- # +# Reaction with an isozyme rule that survives the removal of one of its genes +# --------------------------------------------------------------------------- # + +def _build_light_big_ec_model(adapter: ModelAdapter) -> EcModel: + """Light ecModel: R1 (rule g1) and R2 (rule ``g1 and g2 or g3``), whose + ec rows carry the ``001_`` prefix of a gecko-light model.""" + model = EcModel("big_light", adapter=adapter, gecko_light=True) + mets = {n: cobra.Metabolite(f"{n}_c", compartment="c") for n in "ABC"} + pool = cobra.Metabolite("prot_pool", compartment="c") + model.add_metabolites([*mets.values(), pool]) + + r1 = cobra.Reaction("R1") + r1.add_metabolites({mets["A"]: -1.0, mets["B"]: 1.0}) + r1.gene_reaction_rule = "g1" + r2 = cobra.Reaction("R2") + r2.add_metabolites({mets["B"]: -1.0, mets["C"]: 1.0}) + r2.gene_reaction_rule = "g1 and g2 or g3" + pool_ex = cobra.Reaction("prot_pool_exchange") + pool_ex.add_metabolites({pool: 1.0}) + pool_ex.upper_bound = 1000.0 + model.add_reactions([r1, r2, pool_ex]) + + mat = sparse.lil_matrix((2, 3), dtype=float) + mat[0, 0] = 1.0 + mat[1, 1] = 1.0 + mat[1, 2] = 1.0 + model.ec = EcData( + rxns=["001_R1", "001_R2"], + kcat=np.array([1.0, 2.0]), + source=["initial", "initial"], + notes=["", ""], + eccodes=["1.1.1.1", "2.2.2.2"], + genes=["g1", "g2", "g3"], + enzymes=["E1", "E2", "E3"], + mw=np.array([100.0, 200.0, 300.0]), + sequence=["A", "B", "C"], + concs=np.array([np.nan, np.nan, np.nan]), + rxn_enz_mat=mat.tocsr(), + gecko_light=True, + ) + return model + + +def _small_gem_without_r2_and_g3() -> cobra.Model: + """Context model with R1 and the genes g1 and g2, lacking R2 and g3.""" + small = cobra.Model("small") + r1 = cobra.Reaction("R1") + r1.gene_reaction_rule = "g1" + small.add_reactions([r1]) + small.genes.append(cobra.Gene("g2")) + return small + + +def test_rxn_with_surviving_rule_and_missing_in_small_gem(tmp_path): + """R2 survives the removal of g3 (g1 and g2 satisfy its rule) and is then + removed because the context model lacks it; the gene g3 is not removed twice.""" + big = _build_light_big_ec_model(_adapter(tmp_path)) + result = get_subset_ec_model(big, _small_gem_without_r2_and_g3()) + + assert {r.id for r in result.reactions} == {"R1", "prot_pool_exchange"} + assert {g.id for g in result.genes} == {"g1"} + assert result.ec.genes == ["g1", "g2"] + assert result.ec.enzymes == ["E1", "E2"] + + +def test_surviving_rxn_loses_removed_gene_from_rule(tmp_path): + """A kept reaction no longer refers to a gene that left model.genes.""" + big = _build_light_big_ec_model(_adapter(tmp_path)) + small = _small_gem_without_r2_and_g3() + small.add_reactions([cobra.Reaction("R2")]) + result = get_subset_ec_model(big, small) + + r2 = result.reactions.get_by_id("R2") + assert {g.id for g in r2.genes} == {"g1", "g2"} + assert "g3" not in r2.gene_reaction_rule + assert {g.id for g in result.genes} == {"g1", "g2"} + + +def test_light_model_keeps_ec_rows_of_kept_reactions(tmp_path): + """The ``001_`` prefix of the ec rows of a light model is not part of the + reaction id the rows are matched against.""" + big = _build_light_big_ec_model(_adapter(tmp_path)) + result = get_subset_ec_model(big, _small_gem_without_r2_and_g3()) + + assert result.ec.rxns == ["001_R1"] + np.testing.assert_array_equal(result.ec.kcat, [1.0]) + assert result.ec.eccodes == ["1.1.1.1"] + assert result.ec.rxn_enz_mat.shape == (1, 2) + + +def test_identity_subset_of_light_model_is_a_copy(tmp_path): + big = _build_light_big_ec_model(_adapter(tmp_path)) + small = cobra.Model("small") + for rid, rule in (("R1", "g1"), ("R2", "g1 and g2 or g3")): + rxn = cobra.Reaction(rid) + rxn.gene_reaction_rule = rule + small.add_reactions([rxn]) + result = get_subset_ec_model(big, small) + + assert result is not big + assert sorted(r.id for r in result.reactions) == sorted(r.id for r in big.reactions) + assert result.ec.rxns == ["001_R1", "001_R2"] + + +def test_light_model_ec_rows_without_gene_removal(tmp_path): + big = _build_light_big_ec_model(_adapter(tmp_path)) + small = cobra.Model("small") + r1 = cobra.Reaction("R1") + r1.gene_reaction_rule = "g1" + small.add_reactions([r1]) + small.genes.extend([cobra.Gene("g2"), cobra.Gene("g3")]) + result = get_subset_ec_model(big, small) + + assert [r.id for r in result.reactions if r.id.startswith("R")] == ["R1"] + assert result.ec.rxns == ["001_R1"] + assert result.ec.rxn_enz_mat.shape == (1, 3) + + +# --------------------------------------------------------------------------- # +# Metabolites and genes that no reaction uses +# --------------------------------------------------------------------------- # + +def test_unused_metabolites_and_genes_are_removed(tmp_path): + """MATLAB's removeGenes and removeReactions drop every metabolite and gene + without a reaction, also those that were unused in the big model.""" + big = _build_light_big_ec_model(_adapter(tmp_path)) + big.add_metabolites([cobra.Metabolite("X_c", compartment="c")]) + big.genes.append(cobra.Gene("standard")) + small = cobra.Model("small") + for rid, rule in (("R1", "g1"), ("R2", "g1 and g2 or g3")): + rxn = cobra.Reaction(rid) + rxn.gene_reaction_rule = rule + small.add_reactions([rxn]) + result = get_subset_ec_model(big, small) + + assert not result.metabolites.has_id("X_c") + assert not result.genes.has_id("standard") + assert {g.id for g in result.genes} == {"g1", "g2", "g3"} + assert all(m.reactions for m in result.metabolites) + assert big.metabolites.has_id("X_c") + assert big.genes.has_id("standard") diff --git a/tests/test_load_ec_model.py b/tests/test_load_ec_model.py index 7cb875b..ad9d17a 100644 --- a/tests/test_load_ec_model.py +++ b/tests/test_load_ec_model.py @@ -548,3 +548,110 @@ def read_with_warning(path): monkeypatch.setattr(lem, "read_yaml_model", read_with_warning) with pytest.warns(UserWarning, match="unrelated reader warning"): load_ec_model("ecModel.yml", adapter=adapter) + + +# --------------------------------------------------------------------------- # +# 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` +# 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 = """\ +!!omap +- metaData: !!omap + - id: "t" + - name: "t" +- metabolites: + - !!omap + - id: "A_c" + - name: "A" + - compartment: "c" +- reactions: + - !!omap + - id: "R1" + - name: "r" + - metabolites: !!omap + - A_c: -1 + - lower_bound: 0 + - upper_bound: 1000 + - gene_reaction_rule: "g1" +- genes: + - !!omap + - id: "g1" +- compartments: !!omap + - c: "cyto" +- gecko_light: false +- ec-rxns: + - !!omap + - id: "R1_EXP_1" + - kcat: 1.5 + - source: "standard" + - notes: ""Putative x. And y."" + - enzymes: !!omap + - !!omap + - id: "R2" + - kcat: 2 + - notes: "say \\"hi\\" now" + - enzymes: !!omap + - P1: 1 + - !!omap + - id: "R3" + - kcat: 3 + - notes: "say "hi" now \\\\ ok" + - enzymes: !!omap +- ec-enzymes: + - !!omap + - genes: "g1" + - enzymes: "P1" + - mw: 10 + - sequence: "ABC" +""" + + +def test_enzymes_omap_without_entries_and_unescaped_quotes_load(tmp_path): + """The values are those MATLAB's readYAMLmodel gives: no enzymes for the rows + without entries, the outer pair of quotes of a note removed and the inner ones kept.""" + path = tmp_path / "gecko.yml" + path.write_text(_GECKO_FORM, encoding="utf-8") + model = load_ec_model(path) + + assert model.ec.rxns == ["R1_EXP_1", "R2", "R3"] + np.testing.assert_array_equal(model.ec.kcat, [1.5, 2.0, 3.0]) + assert model.ec.source == ["standard", "", ""] + assert model.ec.notes == ['"Putative x. And y."', 'say "hi" now', 'say "hi" now \\ ok'] + assert model.ec.enzymes == ["P1"] + np.testing.assert_array_equal(model.ec.rxn_enz_mat.toarray(), [[0.0], [1.0], [0.0]]) + + +def test_valid_file_is_not_rewritten(tmp_path): + """A quoted value with escaped quotes stays as it is.""" + from geckopy.utilities.load_ec_model import _repair_yaml_text + + text = ' - notes: "say \\"hi\\" now"\n - enzymes: !!omap\n - P1: 1\n - notes: ""\n' + assert _repair_yaml_text(text) == text + + +def test_file_that_stays_invalid_raises_its_first_error(tmp_path): + from raven_toolbox.io import read_yaml_model + + path = tmp_path / "bad.yml" + path.write_text("!!omap\n- id: [unclosed\n", encoding="utf-8") + with pytest.raises(Exception) as direct: + read_yaml_model(path) + with pytest.raises(type(direct.value)): + load_ec_model(path) + + +def test_enzymes_omap_without_entries_loads(tmp_path): + text = "".join( + line for line in _GECKO_FORM.splitlines(keepends=True) + if "Putative" not in line and "ok" not in line + ) + path = tmp_path / "gecko.yml" + path.write_text(text, encoding="utf-8") + model = load_ec_model(path) + + assert model.ec.rxns == ["R1_EXP_1", "R2", "R3"] + assert model.ec.notes == ["", 'say "hi" now', ""] + np.testing.assert_array_equal(model.ec.rxn_enz_mat.toarray(), [[0.0], [1.0], [0.0]]) diff --git a/tests/test_pfba_enzymes.py b/tests/test_pfba_enzymes.py index 3c0d3d1..fd5aa2f 100644 --- a/tests/test_pfba_enzymes.py +++ b/tests/test_pfba_enzymes.py @@ -48,10 +48,19 @@ def test_returns_solution(): def test_objective_preserved_at_full_fraction(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() - sol = pfba_enzymes(model) + sol = pfba_enzymes(model, fraction_of_optimum=1.0) assert sol.fluxes["R3"] == pytest.approx(fba_optimum, abs=1e-6) +def test_default_fixes_the_objective_at_0_999999_of_the_optimum(): + """MATLAB's getPfbaEnzymes (RAVEN's solveLP) keeps 0.999999 of the optimum.""" + model = _ectestgem_ec_model_with_kcats() + fba_optimum = model.slim_optimize() + sol = pfba_enzymes(model) + ratio = sol.fluxes["R3"] / fba_optimum + assert ratio == pytest.approx(0.999999, abs=2e-7) + + def test_objective_preserved_at_partial_fraction(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() diff --git a/tests/test_review_fixes_r1p.py b/tests/test_review_fixes_r1p.py index 8d52768..43203b4 100644 --- a/tests/test_review_fixes_r1p.py +++ b/tests/test_review_fixes_r1p.py @@ -9,7 +9,7 @@ from test_get_conc_control_coeffs import _build_toy_ec_model from test_sigma_fitter import _adapter, _build_pool_only_model -SOLVERS = ["glpk", "glpk_exact", "gurobi"] +SOLVERS = ["glpk", "glpk_exact", "gurobi", "scipy"] def _with_solver(model, solver): diff --git a/tests/test_review_fixes_r2p.py b/tests/test_review_fixes_r2p.py index 906d9ea..0faab13 100644 --- a/tests/test_review_fixes_r2p.py +++ b/tests/test_review_fixes_r2p.py @@ -15,7 +15,7 @@ from test_get_conc_control_coeffs import _build_toy_ec_model -SOLVERS = ["glpk", "glpk_exact", "gurobi"] +SOLVERS = ["glpk", "glpk_exact", "gurobi", "scipy"] MODULE = "geckopy.limit_proteins.get_conc_control_coeffs" diff --git a/tests/test_scipy_solver_support.py b/tests/test_scipy_solver_support.py new file mode 100644 index 0000000..b0a155f --- /dev/null +++ b/tests/test_scipy_solver_support.py @@ -0,0 +1,507 @@ +"""The scipy (HiGHS) solver gives the results of glpk and gurobi. + +optlang's scipy interface renumbers its variable columns wrongly after a variable is removed, +which makes cobra read fluxes from the wrong columns, and it raises for duals. geckopy repairs +the first and reads the HiGHS marginals for the second (see ``geckopy.solver_support``). +""" +from pathlib import Path + +import cobra +import numpy as np +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.ec_model.pipeline import apply_kcat_constraints +from geckopy.limit_proteins import constrain_enz_concs, relax_proteomics_greedy +from geckopy.solver_support import ( + constraint_shadow_prices, + edit_on_fast_solver, + repair_scipy_problem, + switch_to_scipy, +) +from geckopy.utilities import ec_fva, get_enzyme_bottlenecks, get_subset_ec_model + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" +REFERENCE_SOLVERS = ["glpk", "gurobi"] + + +@pytest.fixture +def scipy_default(): + """cobra builds every new model on scipy while the test runs.""" + config = cobra.Configuration() + previous = config.solver + config.solver = "scipy" + yield + config.solver = previous + + +def _build(solver="scipy", objective="R3"): + """ecTestGEM as an ecModel built on ``solver``, kcat 10, with ``objective``; also the + conventional model it was built from.""" + config = cobra.Configuration() + previous = config.solver + try: + config.solver = solver + except Exception as exc: # solver not installed or not licensed + pytest.skip(f"solver {solver} unavailable: {exc}") + try: + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + conv = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + model = make_ec_model(conv, adapter) + finally: + config.solver = previous + model.ec.kcat[:] = 10.0 + apply_kcat_constraints(model) + model.objective = objective + return model, conv + + +def _on(solver, tight=False, **kwargs): + """Like ``_build``, solved on ``solver``. ``tight`` caps every enzyme at 0.01, which makes + the optimum and its duals unique; without it the LP has several optimal duals and the + solvers pick different ones.""" + model, conv = _build(solver, **kwargs) + if tight: + model.ec.concs[:] = 0.01 + constrain_enz_concs(model) + try: + model.solver = solver + model.optimize() + except Exception as exc: + pytest.skip(f"solver {solver} unavailable: {exc}") + return model, conv + + +def _columns_are_consistent(model): + problem = model.solver.problem + return list(problem._variables.values()) == list(range(len(problem._variables))) + + +def test_columns_are_consistent_after_make_ec_model(scipy_default): + model, _ = _build("scipy") + assert _columns_are_consistent(model) + + +def test_fluxes_after_make_ec_model_equal_glpk(scipy_default): + scipy_model, _ = _build("scipy") + glpk_model, _ = _build("glpk") + scipy_sol = scipy_model.optimize() + glpk_sol = glpk_model.optimize() + assert scipy_sol.objective_value == pytest.approx(glpk_sol.objective_value, rel=1e-7) + assert scipy_sol.fluxes["R3"] == pytest.approx(scipy_sol.objective_value, rel=1e-7) + assert scipy_model.reactions.get_by_id("R3").flux == pytest.approx(scipy_sol.objective_value, rel=1e-7) + + +def test_removing_reactions_keeps_the_fluxes_right(scipy_default): + model, _ = _build("scipy") + reference, _ = _build("glpk") + victims = [r.id for r in model.reactions if r.id.startswith("usage_prot_")][:2] + for rid in victims: + model.reactions.get_by_id(rid).delete() + reference.reactions.get_by_id(rid).delete() + assert model.slim_optimize() == pytest.approx(reference.slim_optimize(), rel=1e-7) + assert _columns_are_consistent(model) + solution = model.optimize() + assert solution.fluxes["R3"] == pytest.approx(solution.objective_value, rel=1e-7) + assert model.reactions.get_by_id("R3").flux == pytest.approx(solution.objective_value, rel=1e-7) + + +def _ring(solver): + """A plain cobra model on ``solver`` with more variables than constraints.""" + config = cobra.Configuration() + previous = config.solver + config.solver = solver + try: + model = cobra.Model("plain") + mets = [cobra.Metabolite(f"m{i}") for i in range(3)] + rxns = [] + for i in range(6): + rxn = cobra.Reaction(f"r{i}") + rxn.bounds = (0, 10) + rxn.add_metabolites({mets[i % 3]: 1.0, mets[(i + 1) % 3]: -1.0} if i else {mets[0]: 1.0}) + rxns.append(rxn) + model.add_reactions(rxns) + finally: + config.solver = previous + model.objective = "r5" + return model + + +def test_a_plain_cobra_model_is_repaired_on_request(): + model, reference = _ring("scipy"), _ring("glpk") + for m in (model, reference): + m.remove_reactions(["r1", "r2"]) + m.reactions.get_by_id("r4").bounds = (0, 3) + assert repair_scipy_problem(model) + assert _columns_are_consistent(model) + assert model.slim_optimize() == pytest.approx(reference.slim_optimize()) + assert model.reactions.get_by_id("r5").flux == pytest.approx(reference.reactions.get_by_id("r5").flux) + + +def test_repair_does_nothing_on_other_solvers(): + model = cobra.Model("m") + model.solver = "glpk" + assert repair_scipy_problem(model) is False + + +@pytest.mark.parametrize("reference", REFERENCE_SOLVERS) +def test_ec_fva_on_scipy_equals_the_reference(reference, scipy_default): + ref_model, ref_conv = _on(reference) + scipy_model, scipy_conv = _on("scipy") + expected = ec_fva(ref_model, ref_conv, progress=False, n_proc=1) + actual = ec_fva(scipy_model, scipy_conv, progress=False, n_proc=1) + assert list(actual.index) == list(expected.index) + assert actual.to_numpy(dtype=float) == pytest.approx(expected.to_numpy(dtype=float), abs=1e-6) + assert np.abs(actual.to_numpy(dtype=float)).max() > 0 + + +@pytest.mark.parametrize("reference", REFERENCE_SOLVERS) +def test_get_enzyme_bottlenecks_on_scipy_equals_the_reference(reference, scipy_default): + expected = get_enzyme_bottlenecks(_on(reference, tight=True)[0], top=999) + scipy_model = _on("scipy", tight=True)[0] + actual = get_enzyme_bottlenecks(scipy_model, top=999) + assert (expected["shadow_price"].abs() > 1e-9).any() + assert list(actual.index[:1]) == list(expected.index[:1]) + assert set(actual.index) == set(expected.index) + aligned = actual.loc[expected.index] + assert aligned["shadow_price"].to_numpy() == pytest.approx( + expected["shadow_price"].to_numpy(), abs=1e-6 + ) + assert actual["shadow_price"].abs().is_monotonic_decreasing + for uniprot, row in actual.iterrows(): + assert row["flux"] == pytest.approx(scipy_model.enzymes.get_by_id(uniprot).flux) + + +@pytest.mark.parametrize("reference", REFERENCE_SOLVERS) +def test_enzyme_shadow_price_on_scipy_equals_the_reference(reference, scipy_default): + ref_model, _ = _on(reference, tight=True) + scipy_model, _ = _on("scipy", tight=True) + for enz in ref_model.enzymes: + assert scipy_model.enzymes.get_by_id(enz.id).shadow_price == pytest.approx( + enz.shadow_price, abs=1e-6 + ) + + +def _tight(model): + model.ec.concs[:] = 1e-6 + constrain_enz_concs(model) + return model + + +@pytest.mark.parametrize("reference", REFERENCE_SOLVERS) +def test_relax_proteomics_greedy_on_scipy_equals_the_reference(reference, scipy_default): + ref_model, _ = _on(reference) + scipy_model, _ = _on("scipy") + expected = relax_proteomics_greedy(_tight(ref_model), minimal_growth=0.001, max_iterations=20) + actual = relax_proteomics_greedy(_tight(scipy_model), minimal_growth=0.001, max_iterations=20) + assert expected.converged and actual.converged + assert [s.relaxed_uniprot for s in actual.trace] == [s.relaxed_uniprot for s in expected.trace] + assert [s.shadow_price for s in actual.trace] == pytest.approx( + [s.shadow_price for s in expected.trace], abs=1e-6 + ) + assert actual.final_growth == pytest.approx(expected.final_growth, rel=1e-7) + + +def test_a_stored_result_from_before_a_change_is_not_read(scipy_default): + model, _ = _on("scipy", tight=True) + ids = [enz.prot_metabolite_id for enz in model.enzymes] + before = constraint_shadow_prices(model, ids) + usage = next(r for r in model.reactions if r.id.startswith("usage_prot_")) + usage.upper_bound = usage.upper_bound * 10 + 1.0 + after = constraint_shadow_prices(model, ids) + reference, _ = _on("glpk", tight=True) + reference.reactions.get_by_id(usage.id).upper_bound = usage.upper_bound + reference.optimize() + expected = np.array([float(reference.constraints[i].dual) for i in ids]) + assert after == pytest.approx(expected, abs=1e-6) + assert before.shape == after.shape + + +def _columns_ring(solver): + """Three metabolites, twelve reactions: removals leave late variables with stale columns.""" + config = cobra.Configuration() + previous = config.solver + config.solver = solver + try: + model = cobra.Model("ring") + mets = [cobra.Metabolite(f"m{i}") for i in range(3)] + rxns = [] + for i in range(12): + rxn = cobra.Reaction(f"r{i}") + rxn.bounds = (0, 10 + i) + if i == 0: + rxn.add_metabolites({mets[0]: 1.0}) + else: + rxn.add_metabolites({mets[i % 3]: -1.0, mets[(i + 1) % 3]: 0.5 + 0.1 * i}) + rxns.append(rxn) + model.add_reactions(rxns) + finally: + config.solver = previous + model.objective = "r11" + return model + + +def test_removals_after_the_install_keep_every_column_right(): + model, reference = _columns_ring("scipy"), _columns_ring("glpk") + assert repair_scipy_problem(model) + for rid in ("r3", "r7", "r9", "r10"): + model.reactions.get_by_id(rid).delete() + reference.reactions.get_by_id(rid).delete() + assert _columns_are_consistent(model) + assert model.slim_optimize() == pytest.approx(reference.slim_optimize()) + solution = model.optimize() + expected = reference.optimize() + assert solution.objective_value == pytest.approx(expected.objective_value) + for rid in solution.fluxes.index: + assert model.reactions.get_by_id(rid).flux == pytest.approx(solution.fluxes[rid]) + + +def test_a_map_that_is_already_stale_is_rebuilt_from_the_expressions(): + model, reference = _columns_ring("scipy"), _columns_ring("glpk") + model.reactions.get_by_id("r5").delete() + reference.reactions.get_by_id("r5").delete() + model.solver.update() + assert not _columns_are_consistent(model) + assert repair_scipy_problem(model) + assert _columns_are_consistent(model) + assert model.slim_optimize() == pytest.approx(reference.slim_optimize()) + + +def test_the_sparse_problem_hands_highs_the_matrix_of_the_dense_problem(): + model = _columns_ring("scipy") + dense = model.solver.problem.A.copy() + assert repair_scipy_problem(model) + assert (model.solver.problem.A.toarray() == dense).all() + + +def test_switch_to_scipy_keeps_the_objective_and_its_direction(): + model, reference = _columns_ring("glpk"), _columns_ring("glpk") + for m in (model, reference): + m.objective = "r8" + m.objective_direction = "min" + m.reactions.get_by_id("r1").lower_bound = 2.0 + switch_to_scipy(model) + assert model.solver.interface.__name__.endswith("scipy_interface") + assert model.objective_direction == "min" + assert model.slim_optimize() == pytest.approx(reference.slim_optimize()) + assert model.reactions.get_by_id("r8").flux == pytest.approx(model.slim_optimize()) + + +def test_load_conventional_gem_and_make_ec_model_end_on_scipy(scipy_default): + from geckopy import load_conventional_gem + + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + conv = load_conventional_gem(adapter) + assert conv.solver.interface.__name__.endswith("scipy_interface") + reference = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + for m in (conv, reference): + m.reactions.get_by_id("S1").bounds = (-10.0, 10.0) + m.objective = "R5" + assert conv.slim_optimize() == pytest.approx(reference.slim_optimize()) + assert reference.slim_optimize() > 0 + model = make_ec_model(conv, adapter) + assert model.solver.interface.__name__.endswith("scipy_interface") + assert cobra.Configuration().solver.__name__.endswith("scipy_interface") + + +def test_shadow_prices_need_an_optimal_solution(scipy_default): + model, _ = _on("scipy", tight=True) + model.reactions.get_by_id("R3").lower_bound = 1e6 + with pytest.raises(RuntimeError, match="optimal"): + constraint_shadow_prices(model, [model.enzymes.get_by_id(model.ec.enzymes[0]).prot_metabolite_id]) + + +@pytest.mark.parametrize("status_case", ["infeasible", "unbounded"]) +def test_a_solve_without_optimum_returns_a_solution_with_its_status(status_case): + config = cobra.Configuration() + previous = config.solver + config.solver = "scipy" + try: + model = cobra.Model("chain") + met = cobra.Metabolite("m") + source, sink = cobra.Reaction("source"), cobra.Reaction("sink") + source.add_metabolites({met: 1.0}) + sink.add_metabolites({met: -1.0}) + source.bounds = (0, 10) if status_case == "infeasible" else (0, 1000.0) + sink.bounds = (20, 30) if status_case == "infeasible" else (0, 1000.0) + model.add_reactions([source, sink]) + finally: + config.solver = previous + model.objective = "sink" + if status_case == "unbounded": + source.upper_bound = float("inf") + sink.upper_bound = float("inf") + assert repair_scipy_problem(model) + if status_case == "unbounded": + # cobra raises for this status on every solver + with pytest.raises(cobra.exceptions.OptimizationError, match="unbounded"): + model.optimize(raise_error=False) + return + solution = model.optimize(raise_error=False) + assert solution.status == status_case + assert np.isnan(solution.objective_value) + assert solution.fluxes.isna().all() + + +def test_setting_the_solver_of_an_ec_model_to_scipy_gives_the_sparse_problem(): + from geckopy.solver_support import _SparseProblem + + model, _ = _build("glpk") + expected = model.slim_optimize() + model.solver = "scipy" + assert isinstance(model._solver.problem, _SparseProblem) + assert model.slim_optimize() == pytest.approx(expected) + + +def test_load_ec_model_ends_on_the_sparse_scipy_problem(scipy_default, tmp_path): + from geckopy import load_ec_model, save_ec_model + from geckopy.solver_support import _SparseProblem + + model, _ = _build("glpk") + expected = model.slim_optimize() + save_ec_model(model, tmp_path / "ecModel.yml") + loaded = load_ec_model(tmp_path / "ecModel.yml") + assert isinstance(loaded._solver.problem, _SparseProblem) + assert cobra.Configuration().solver.__name__.endswith("scipy_interface") + loaded.objective = "R3" + assert loaded.slim_optimize() == pytest.approx(expected) + + +@pytest.mark.parametrize("how", ["pickle", "deepcopy", "copy"]) +def test_a_scipy_ec_model_pickles_and_copies_into_the_sparse_problem(how): + import copy + import pickle + + from geckopy.solver_support import _SparseProblem + + model, _ = _build("glpk") + expected = model.slim_optimize() + model.solver = "scipy" + model.slim_optimize() + clone = { + "pickle": lambda: pickle.loads(pickle.dumps(model)), + "deepcopy": lambda: copy.deepcopy(model), + "copy": model.copy, + }[how]() + assert isinstance(clone._solver.problem, _SparseProblem) + assert clone.solver.interface.__name__.endswith("scipy_interface") + assert clone.slim_optimize() == pytest.approx(expected) + assert clone.reactions.get_by_id("R3").flux == pytest.approx(model.reactions.get_by_id("R3").flux) + + +def test_model_tolerance_reaches_highs_and_survives_a_copy(monkeypatch): + import pickle + + from scipy import optimize + + model, _ = _build("scipy") + model.tolerance = 1e-9 + assert model.solver.configuration.tolerances.feasibility == 1e-9 + assert model.solver.configuration.tolerances.optimality == 1e-9 + + seen = [] + real = optimize.linprog + + def spy(*args, **kwargs): + seen.append(kwargs["options"]) + return real(*args, **kwargs) + + for m in (model, pickle.loads(pickle.dumps(model)), model.copy()): + monkeypatch.setattr("geckopy.solver_support.linprog", spy) + m.slim_optimize() + assert len(seen) == 3 + assert all(o["primal_feasibility_tolerance"] == 1e-9 for o in seen) + assert all(o["dual_feasibility_tolerance"] == 1e-9 for o in seen) + + +def test_a_solve_with_numerical_difficulties_is_repeated_at_the_default_tolerance(monkeypatch): + from scipy import optimize + + model, _ = _build("scipy") + expected = model.slim_optimize() + model.tolerance = 1e-9 + real = optimize.linprog + calls = [] + + def numerical_trouble_when_tight(*args, **kwargs): + calls.append(kwargs["options"]) + if "primal_feasibility_tolerance" in kwargs["options"]: + return optimize.OptimizeResult(status=4, fun=None, x=None, slack=None) + return real(*args, **kwargs) + + monkeypatch.setattr("geckopy.solver_support.linprog", numerical_trouble_when_tight) + assert model.slim_optimize() == pytest.approx(expected) + assert len(calls) == 2 + + def always_numerical_trouble(*args, **kwargs): + return optimize.OptimizeResult(status=4, fun=None, x=None, slack=None) + + monkeypatch.setattr("geckopy.solver_support.linprog", always_numerical_trouble) + assert model.optimize(raise_error=False).status == "numeric" + + +def _subset(solver, light): + config = cobra.Configuration() + adapter = ModelAdapter.from_folder(EXAMPLE_DIR) + try: + config.solver = solver + except Exception as exc: + pytest.skip(f"solver {solver} unavailable: {exc}") + try: + conv = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + big = make_ec_model(conv, adapter, gecko_light=light) + small = conv.copy() + small.remove_reactions(["R2"], remove_orphans=True) + sub = get_subset_ec_model(big, small) + finally: + config.solver = "glpk" + big.ec.kcat[:] = 10.0 + sub.ec.kcat[:] = 10.0 + apply_kcat_constraints(sub) + sub.objective = "R3" + return sub + + +@pytest.mark.parametrize("light", [False, True]) +def test_get_subset_ec_model_on_scipy_equals_glpk(light): + on_scipy, on_glpk = _subset("scipy", light), _subset("glpk", light) + assert on_scipy.solver.interface.__name__.endswith("scipy_interface") + assert {r.id for r in on_scipy.reactions} == {r.id for r in on_glpk.reactions} + assert {m.id for m in on_scipy.metabolites} == {m.id for m in on_glpk.metabolites} + assert {g.id for g in on_scipy.genes} == {g.id for g in on_glpk.genes} + assert len(on_scipy.solver.constraints) == len(on_scipy.metabolites) + assert on_scipy.slim_optimize() == pytest.approx(on_glpk.slim_optimize()) + assert on_scipy.slim_optimize() > 0 + + +def test_removing_several_constraints_at_once_on_the_sparse_problem(): + model = _columns_ring("scipy") + repair_scipy_problem(model) + ids = [m.id for m in list(model.metabolites)[:3]] + model.remove_cons_vars([model.constraints[i] for i in ids]) + model.solver.update() + assert len(model.solver.problem._constraints) == 2 * len(model.constraints) + assert not set(ids) & set(model.constraints) + + +def test_edit_on_fast_solver_returns_to_scipy_with_the_same_problem(): + model = _columns_ring("scipy") + repair_scipy_problem(model) + before = model.slim_optimize() + with edit_on_fast_solver(model) as fast: + assert fast + assert model.solver.interface.__name__.endswith("glpk_interface") + model.reactions.get_by_id("r1").bounds = (0.0, 5.0) + assert model.solver.interface.__name__.endswith("scipy_interface") + assert model.reactions.get_by_id("r1").upper_bound == 5.0 + assert model.slim_optimize() <= before + reference = _columns_ring("glpk") + reference.reactions.get_by_id("r1").bounds = (0.0, 5.0) + assert model.slim_optimize() == pytest.approx(reference.slim_optimize()) + + +def test_edit_on_fast_solver_changes_nothing_on_glpk(): + model = _columns_ring("glpk") + with edit_on_fast_solver(model) as fast: + assert not fast + assert model.solver.interface.__name__.endswith("glpk_interface") diff --git a/tests/test_sensitivity_tuning.py b/tests/test_sensitivity_tuning.py index 294aa8f..5d8da43 100644 --- a/tests/test_sensitivity_tuning.py +++ b/tests/test_sensitivity_tuning.py @@ -251,3 +251,41 @@ def test_result_reports_its_own_corrections(): assert rows[0].name == "drifting enzyme" assert rows[1].source == "brenda" assert rows[1].fold_change == pytest.approx(10.0) + + +# --------------------------------------------------------------------------- # +# Infeasible start: the same on every solver +# --------------------------------------------------------------------------- # + +def _usable_solvers(): + solvers = [] + for name in ("glpk", "glpk_exact", "gurobi"): + try: + cobra.Model("probe").solver = name + except Exception: + continue + solvers.append(name) + return solvers + + +@pytest.fixture(params=_usable_solvers()) +def solver(request): + """Solver name; a solver that imports but is not licensed is skipped.""" + model = cobra.Model("probe") + try: + model.solver = request.param + model.add_metabolites([cobra.Metabolite("m")]) + model.optimize() + except Exception as exc: + pytest.skip(f"{request.param} unusable: {exc}") + return request.param + + +def test_infeasible_start_raises_on_every_solver(tmp_path, solver): + adapter = _adapter(tmp_path, gr_exp=300.0) + model = _build_toy(adapter, initial_kcat=1.0, enzyme_mw=100.0) + model.solver = solver + model.reactions.get_by_id("biomass").lower_bound = 500.0 + assert model.optimize().status == "infeasible" + with pytest.raises(ValueError, match="not optimal"): + sensitivity_tuning(model) diff --git a/tests/test_uniprot_loader.py b/tests/test_uniprot_loader.py index 7497864..e41496c 100644 --- a/tests/test_uniprot_loader.py +++ b/tests/test_uniprot_loader.py @@ -234,3 +234,79 @@ def test_find_by_id_returns_index(tmp_path): db = load_uniprot_tsv(path) assert db.find_by_id("P2") == 1 assert db.find_by_id("Pxxx") is None + + +# --------------------------------------------------------------------------- # +# Quoted cells, read as MATLAB's textscan with %q reads them +# --------------------------------------------------------------------------- # + +T = "\t" + + +def test_quoted_cell_is_given_without_its_quotes(tmp_path): + path = _write_tsv(tmp_path / "u.tsv", [ + 'P1\tG1\t"1.1.1.25; 2.5.1.19"\t10000\tMKAL', + ]) + db = load_uniprot_tsv(path) + assert db.eccodes == ["1.1.1.25; 2.5.1.19"] + assert db.genes == ["G1"] + assert db.sequences == ["MKAL"] + + +def test_quoted_gene_cell_is_split_without_quotes(tmp_path): + path = _write_tsv(tmp_path / "u.tsv", [ + 'P1\t"G1 G2"\t\t10000\tMKAL', + ]) + db = load_uniprot_tsv(path, split_gene_cells=True) + assert db.genes == ["G1", "G2"] + + +def test_quoted_entry_ids_count_as_duplicates(tmp_path): + path = _write_tsv(tmp_path / "u.tsv", [ + 'P1\tG1\t\t10000\tMKAL', + '"P1"\tG2\t\t10000\tMKAL', + ]) + with pytest.raises(ValueError, match="Duplicate"): + load_uniprot_tsv(path) + + +# The expected rows are what MATLAB's textscan(fid, '%q %q %q %q %q', 'Delimiter', '\t', +# 'HeaderLines', 1) returns for the same text (MATLAB R2024b). +@pytest.mark.parametrize( + "text, rows", + [ + (f'P1{T}g1{T}"1.1.1.25; 2.5.1.19"{T}100{T}SEQ\n', [["P1", "g1", "1.1.1.25; 2.5.1.19", "100", "SEQ"]]), + (f'P2{T}"a ""b"" c"{T}{T}10{T}S\n', [["P2", 'a "b" c', "", "10", "S"]]), + (f'P3{T}"has{T}tab"{T}"x"{T}5{T}S\n', [["P3", f"has{T}tab", "x", "5", "S"]]), + (f'P4{T}ab"cd"{T}"x"y{T}7{T}S\n', [["P4", 'ab"cd"', "xy", "7", "S"]]), + (f'P6{T}" spaced "{T} lead {T}4{T}S\n', [["P6", " spaced ", "lead ", "4", "S"]]), + (f"P7{T}g\n", [["P7", "g", "", "", ""]]), + (f'P12{T}"""a"""{T}x{T}1{T}S\n', [["P12", '"a"', "x", "1", "S"]]), + (f'P13{T}a""b{T}"x\ny"{T}1{T}S\n', [["P13", 'a""b', "x\ny", "1", "S"]]), + (f'P10{T}"unterminated{T}x{T}1{T}S\n', [["P10", f"unterminated{T}x{T}1{T}S\n", "", "", ""]]), + (f'A{T}g{T}e{T}5{T}S\nB{T}h\nC{T}i{T}f{T}6{T}T\n', + [["A", "g", "e", "5", "S"], ["B", "h", "", "", ""], ["C", "i", "f", "6", "T"]]), + (f"A{T}g{T}e{T}5{T}S{T}extra{T}more\nB{T}h{T}e{T}6{T}T\n", + [["A", "g", "e", "5", "S"], ["extra", "more", "", "", ""], ["B", "h", "e", "6", "T"]]), + (f'A{T}"x y"{T}e{T}5{T}S\r\nB{T}h{T}e{T}6{T}T\r\n', + [["A", "x y", "e", "5", "S"], ["B", "h", "e", "6", "T"]]), + (f"A{T}g{T}{T}{T}\nB{T}h{T}{T}6{T}\n", + [["A", "g", "", "", ""], ["B", "h", "", "6", ""]]), + (f"A{T}g{T}e{T}5{T}S\n \nB{T}h{T}e{T}6{T}T\n", + [["A", "g", "e", "5", "S"], ["B", "h", "e", "6", "T"]]), + (f'A{T}"x" {T}e{T}5{T}S\n', [["A", "x ", "e", "5", "S"]]), + (f'A{T} "x y"{T}e{T}5{T}S\n', [["A", "x y", "e", "5", "S"]]), + (f"A{T}g{T}e{T}5{T}S", [["A", "g", "e", "5", "S"]]), + (f"A{T}'x y'{T}e{T}5{T}S\n", [["A", "'x y'", "e", "5", "S"]]), + ], +) +def test_read_quoted_tsv_matches_matlab_textscan(text, rows): + from geckopy.databases.uniprot_loader import read_quoted_tsv + + assert read_quoted_tsv(text, 5) == rows + + +def test_read_quoted_tsv_skips_header_lines(): + from geckopy.databases.uniprot_loader import read_quoted_tsv + + assert read_quoted_tsv(f"a{T}b\r\nx{T}y\n", 2, header_lines=1) == [["x", "y"]]