Same results as GECKO for subsets, subsystems, uniprot cells, FSEOF and pFBA; scipy solver support

.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

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)

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]
 

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

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,