fix: align raven-toolbox with RAVEN (parity review, round 3)

CHANGELOG.md

@@ -6,6 +6,70 @@ Milestones in the raven-toolbox port. For function-level status see
 
 ## Unreleased
 
+* `manipulation.remove_genes` defaults to `blocked_reactions="constrain"` (was: `"remove"`), as
+  RAVEN's `removeGenes` does: a reaction whose rule loses its last enzyme stays, with the rule
+  dropped and bounds `(0, 0)`. Pass `blocked_reactions="remove"` to delete such reactions.
+* `manipulation.remove_genes` also removes every gene that no remaining reaction named in its
+  rule before the removal, also when not requested (was: only the requested genes), as RAVEN's
+  `removeGenes` does, when at least one of the requested genes is in the model; when none is,
+  the model is not changed. A gene kept only by a constrained reaction stays.
+* `init.merge_linear` numbers the merged groups 1, 2, ... in the order of the first reaction
+  of each group in the original model, as RAVEN's `mergeLinear` does, so `group_ids` no
+  longer depends on the merge order.
+* `manipulation.merge_compartments` gives a merged metabolite the id of the copy already in
+  the merged compartment, if any, otherwise the id that sorts first (`ACCOA_c` for `ACCOA_c`
+  and `ACCOA_m`), as RAVEN's `mergeCompartments` does. It rewrote the first member's
+  compartment suffix to the merged compartment (`ACCOA_s`). The merged metabolite also keeps
+  the annotation of that copy (was: dropped), and `drop_single_metabolite_reactions` defaults to
+  `False`, as RAVEN's `deleteRxnsWithOneMet` (was: `True`).
+* `manipulation.merge_compartments` stores the base id of each merged metabolite in
+  `notes["merged_base_id"]` (`ACCOA` for the kept `ACCOA_c`; the own id when two merged
+  metabolites would get the same base id). The default metabolite key of
+  `localization.assign_compartments`, `apply_assignment`, `relocate_reactions` and
+  `curation_priority` is that stored base id, otherwise the id without the suffix of the
+  metabolite's own compartment, so the merged `ACCOA_c` in `s` and a universal `ACCOA_m` share
+  the key `ACCOA`. The four functions raise `ValueError` when two metabolites in one
+  compartment have the same key (was: one of them stood for both).
+* `manipulation.merge_compartments` documents where `notes["merged_base_id"]` survives: YAML
+  written and read by raven-toolbox, `copy()`, cobra JSON and SBML. RAVEN's `readYAMLmodel`
+  shows it in `metNotes` (`merged_base_id: ACCOA note: free text note`) and no RAVEN function
+  uses it. A round trip through RAVEN YAML or a .mat file loses it, and the localization
+  functions then key a merged `ACCOA_c` in `s` as `ACCOA_c`, without an error.
+* `io.read_yaml_model` puts the notes that `write_yaml_model` writes as a `notes` mapping (notes
+  without a RAVEN field) back into `notes` (was: nested under `notes["note"]`).
+* `reconstruction.homology.run_blast`, `run_diamond` and `run_supporth`, and
+  `reconstruction.kegg.run_hmmsearch`, `get_kegg_model_from_sequences`, `build_ko_hmm` and
+  `build_hmm_library`, use every CPU the process may run on when `threads` is not given, as
+  RAVEN's `getBlast`, `getDiamond`, `runSuppOrth` and `getKEGGModelForOrganism` do (was: the
+  machine's CPU count minus one).
+* `gapfilling.analyse_topology` takes RAVEN's `gapFillTopological` defaults. Seeds: the
+  metabolites of every reaction with only substrates and `lower_bound < 0` or only products
+  and `upper_bound > 0`, whatever its number of metabolites (was: one-metabolite reactions with
+  `lower_bound < 0`); a `UserWarning` when there is none (was: a message only when `verbose`).
+  Targets: the metabolites of the first objective reaction, substrates and products (was: every
+  metabolite of the model); `ValueError` when the model has no objective.
+* New: `utils.parse_gr_rule` reads a grRule in the grammar of RAVEN's `parseGrRule`
+  (`and`/`or` in any case, `&`, `&&`, `|`, `||`, operators as whole words). `manipulation.gpr_to_dnf`
+  now also takes a rule string and reads it with this parser, as RAVEN's `grRuleToDNF` does;
+  `utils.is_dnf` reads a string with it too, and both raise `ValueError` with the message of
+  RAVEN's `RAVEN:badGrRule` for a malformed rule. `G1 && G2`, `G1 || G2` and `G1 Or G2` are
+  read as rules, `G1&G2` is one gene id. Whitespace is what MATLAB's `isspace` counts, so a
+  no-break space (U+00A0), U+0085, U+2007 and U+202F are part of a word.
+* `curation.batch_curate` (and `batch_curate_from_tsv`) stores a gene rule as `parse_gr_rule`
+  reads it, so the reaction gets the genes that were checked: `G1&G2 or G3` has the genes
+  `G1&G2` and `G3`, as in RAVEN's `curateModelFromTables` (was: `G1`, `G2` and `G3`).
+* `io.read_yaml_model` reads gene rules with `parse_gr_rule`, so `G1&G2` written by
+  `write_yaml_model` is read back as one gene (was: `G1` and `G2`, cobra's reading). A rule that
+  `parse_gr_rule` refuses is read by cobra as before.
+* Functions that copy or build a gene rule keep RAVEN's reading of it, so `G1&G2` stays one
+  gene (was: the rule was copied as a string and cobra split `G1&G2` into `G1` and `G2`):
+  `manipulation.merge_compartments`, `copy_to_compartment`, `convert_to_irreversible`,
+  `merge_models`, `add_reactions_from_model`, `expand_model`, `change_gene_reaction_rules`,
+  `add_reactions_from_equations`, `localization.apply_assignment`, `init.fill_tasks`,
+  `init.remove_low_score_genes`, `reconstruction.homology.get_model_from_homology`,
+  `reconstruction.kegg.assemble_model_from_ko_genes` and `io.read_matlab_model`. A rule string
+  is read with `parse_gr_rule`; one that `parse_gr_rule` refuses is read by cobra. cobra's JSON
+  reader (`cobra.io.load_json_model`) still splits `G1&G2`; cobra's SBML reader keeps it.
 * `analysis.fseof` (and so `geckopy.ec_fseof`) solves each level with pFBA that holds the objective at 0.999999 of its optimum at that level, as RAVEN's `solveLP` does, where it held the full value. On the scipy solver (HiGHS) the optimum is only accurate to the solver's feasibility tolerance, so the second step was infeasible at the top level of ecYeastGEM and the scan returned 9 of 10 levels; it now returns all 10, as on Gurobi and GLPK. A level at which the biomass maximisation is infeasible still ends the scan.
 * `analysis.fseof` takes the top of the scan as `max_fraction` of the target flux in the pFBA solution that maximises the target, as RAVEN's `FSEOF` and GECKO's `ecFSEOF` do, where it took the plain maximum (about 1e-6 relative difference).
 * `analysis.fseof` lists targets of equal score in model order (a stable sort), as MATLAB's `sort` does, where their order was arbitrary.
@@ -23,6 +87,25 @@ Milestones in the raven-toolbox port. For function-level status see
   counts the reaction as essential only if that solve is not feasible either, as RAVEN's
   `solveLP` does. Gurobi could report a feasible knock-out as infeasible when started from
   the previous knock-out's basis, which added reactions to the essential set.
+* `init.ftinit` and `init.run_ftinit` with `resolve_ties` take each step's fluxes from the
+  least-total-flux solution for the tie-broken on set, so the direction and magnitude at
+  which the next step forces a reaction no longer depend on the solver's vertex.
+  That phase stops at a relative or absolute gap of 1e-6 and has the full `time_limit`, as
+  in RAVEN.
+* `init.ftinit` and `init.run_ftinit` with `resolve_ties`: among the least-total-flux
+  solutions, each step takes its fluxes from the one with the least id-rank-weighted backward
+  flux, so the direction of a flux that carries the same total either way no longer depends
+  on the solver. The id-rank phases rank a merged reaction by the smallest original reaction
+  id in its group (new `run_ftinit` argument `rank_keys`), so they no longer depend on which
+  member names it.
+* `init.ftinit` and `init.run_ftinit` with `resolve_ties` raise a `UserWarning` naming the
+  phase (`phase2c-minflux`, `phase2d-direction`) when the minimal-flux or the direction phase
+  finds no solution, as RAVEN does, and keep the fluxes of the phase before it. A phase that
+  stops at the time limit is adopted with a `UserWarning` that its fluxes are an unproven
+  incumbent.
+* `init.ftinit` with `resolve_ties`: the last step skips the minimal-flux and direction phases,
+  as no step reads its fluxes, and reports the fluxes of the id-rank phase. The reaction set
+  does not change.
 
 ## 3.0.0b1 — 2026-09-07
 

scripts/benchmark_certified_multiorg.py

@@ -35,6 +35,7 @@
     load_deeploc,
 )
 from raven_toolbox.manipulation.compartments import merge_compartments  # noqa: E402
+from raven_toolbox.utils.gpr import set_gene_rule  # noqa: E402
 
 
 def _fix_notes_gpr(model):
@@ -44,7 +45,7 @@ def _fix_notes_gpr(model):
         ga = r.notes.get("GENE_ASSOCIATION")
         if ga:
             rule = re.sub(r"\bOR\b", "or", re.sub(r"\bAND\b", "and", ga)).strip().strip("()").strip()
-            r.gene_reaction_rule = rule
+            set_gene_rule(r, rule)
     cobra.manipulation.remove_genes(model, [g for g in model.genes if not g.reactions], remove_reactions=False)
 
 

scripts/benchmark_replicate_yeast_gem.py

@@ -122,10 +122,10 @@ def build_draft(yeast: cobra.Model) -> tuple[cobra.Model, str]:
     """Flatten curated yeast-GEM to one compartment via the real `merge_compartments`, and return the
     reaction to use as the biomass/growth objective on the flattened draft.
 
-    `drop_single_metabolite_reactions=False` is essential here: exchange reactions are
-    single-metabolite, so the default would delete every nutrient uptake and the flattened draft
-    couldn't grow (it keeps the collapsed one-metabolite transports too, but those are negligible —
-    the only material effect is preserving the medium). With them kept, yeast-GEM's objective
+    `drop_single_metabolite_reactions=False` (the default, passed explicitly) is essential here:
+    exchange reactions are single-metabolite, so `True` would delete every nutrient uptake and the
+    flattened draft couldn't grow (it keeps the collapsed one-metabolite transports too, but those
+    are negligible; the only material effect is preserving the medium). With them kept, yeast-GEM's objective
     ``r_2111`` "growth" (a ``biomass ->`` drain) survives directly. The `add_boundary` branch is a
     defensive fallback for models whose objective still doesn't survive: re-add an equivalent demand
     sink on the merged biomass metabolite so its producer isn't mass-balance-blocked."""

src/raven_toolbox/curation/batch.py

@@ -12,7 +12,9 @@
 import cobra
 import numpy as np
 import pandas as pd
+from cobra.core.gene import GPR
 
+from raven_toolbox.utils.gpr import parse_gr_rule, rule_tokens
 from raven_toolbox.utils.ids import next_ids
 from raven_toolbox.utils.parse import subsystem_to_str
 
@@ -177,6 +179,10 @@ def batch_curate(
     called without permission to add genes and stops with an error of its own. The rule of an
     existing reaction may name a gene the model does not have. A gene id that ``genes_df``
     lists more than once raises :class:`ValueError`, whether the model has the gene or not.
+    The reaction gets the rule as this grammar reads it, so ``G1&G2`` is one gene, and it stays
+    one gene through :func:`~raven_toolbox.io.write_yaml_model` and
+    :func:`~raven_toolbox.io.read_yaml_model` and through cobra's SBML files. cobra's JSON
+    reader reads the rule with cobra's own grammar and splits ``G1&G2`` into ``G1`` and ``G2``.
 
     A reaction table with a header and no row is empty: when both reaction tables are empty
     nothing happens, and a coefficient row without a reaction row (or the reverse) raises
@@ -434,75 +440,14 @@ def _row_has_content(row, columns=None) -> bool:
     return any(_has_value(v) for v in values)
 
 
-_RULE_TOKEN = re.compile(r"\(|\)|[^\s()]+")
-_AND_WORDS = frozenset({"and", "&", "&&"})
-_OR_WORDS = frozenset({"or", "|", "||"})
-#: Spellings of the operators that cobra reads as they are.
-_COBRA_OPERATORS = frozenset({"and", "AND", "&", "or", "OR", "|"})
-
-
-def _rule_tokens(rule: str) -> list[tuple[str, str]]:
-    """Split a gene rule into ``(kind, text)`` pairs; kind is ``(``, ``)``, ``and``, ``or`` or ``gene``.
-
-    The operators are ``and``/``or`` in any case and the symbols ``&``, ``&&``, ``|`` and ``||``,
-    as RAVEN's ``parseGrRule`` reads them; every other word is a gene id.
-    """
-    tokens = []
-    for match in _RULE_TOKEN.finditer(rule):
-        word = match.group()
-        low = word.lower()
-        if word in ("(", ")"):
-            tokens.append((word, word))
-        elif low in _AND_WORDS:
-            tokens.append(("and", word))
-        elif low in _OR_WORDS:
-            tokens.append(("or", word))
-        else:
-            tokens.append(("gene", word))
-    return tokens
-
-
 def _rule_genes(rule: str) -> list[str]:
-    """Gene ids of a gene rule; :class:`ValueError` when the rule is not well formed.
+    """Gene ids of a gene rule, in order; :class:`ValueError` when the rule is not well formed.
 
-    The grammar is that of RAVEN's ``parseGrRule``: ``expr := term (or term)*``,
-    ``term := factor (and factor)*``, ``factor := gene | ( expr )``.
+    The rule is read with :func:`~raven_toolbox.utils.parse_gr_rule`, the grammar of RAVEN's
+    ``parseGrRule``.
     """
-    tokens = _rule_tokens(rule)
-    genes: list[str] = []
-
-    def factor(pos: int) -> int:
-        if pos >= len(tokens):
-            raise ValueError("the rule ends after an operator")
-        kind, text = tokens[pos]
-        if kind == "(":
-            pos = expr(pos + 1)
-            if pos >= len(tokens) or tokens[pos][0] != ")":
-                raise ValueError("unbalanced parentheses")
-            return pos + 1
-        if kind == "gene":
-            genes.append(text)
-            return pos + 1
-        raise ValueError(f'expected a gene but found "{text}"')
-
-    def term(pos: int) -> int:
-        pos = factor(pos)
-        while pos < len(tokens) and tokens[pos][0] == "and":
-            pos = factor(pos + 1)
-        return pos
-
-    def expr(pos: int) -> int:
-        pos = term(pos)
-        while pos < len(tokens) and tokens[pos][0] == "or":
-            pos = term(pos + 1)
-        return pos
-
-    if not tokens:
-        return genes
-    end = expr(0)
-    if end < len(tokens):
-        raise ValueError(f'unexpected "{tokens[end][1]}"')
-    return genes
+    parse_gr_rule(rule)
+    return [text for kind, text in rule_tokens(rule) if kind == "gene"]
 
 
 def _row_rule_genes(row: pd.Series) -> list[str]:
@@ -518,28 +463,12 @@ def _row_rule_genes(row: pd.Series) -> list[str]:
         ) from None
 
 
-def _stored_gene_rule(rule) -> str:
-    """The rule as cobra reads it: operator spellings cobra does not know (``And``, ``&&``,
-    ``||``) are written as ``and``/``or``; everything else is kept as it is."""
-    text = str(rule)
-    pieces: list[str] = []
-    last = 0
-    for match in _RULE_TOKEN.finditer(text):
-        word = match.group()
-        low = word.lower()
-        if word in _COBRA_OPERATORS or word in ("(", ")"):
-            continue
-        if low in _AND_WORDS:
-            replacement = "and"
-        elif low in _OR_WORDS:
-            replacement = "or"
-        else:
-            continue
-        pieces.append(text[last:match.start()])
-        pieces.append(replacement)
-        last = match.end()
-    pieces.append(text[last:])
-    return "".join(pieces)
+def _stored_gpr(rule) -> GPR:
+    """The rule as :func:`~raven_toolbox.utils.parse_gr_rule` reads it, so the reaction gets the
+    genes that were validated: ``G1&G2`` stays one gene id, as in RAVEN, where cobra reading the
+    string would split it."""
+    parsed = parse_gr_rule(str(rule))
+    return parsed if parsed is not None else GPR()
 
 
 def _check_gene_rules(rxns_df: pd.DataFrame) -> None:
@@ -909,7 +838,7 @@ def _update_reaction(
     if _is_written(row, "rxnNames", overwrite_with_empty):
         rxn.name = _text(row["rxnNames"])
     if _has_value(row.get("grRules")):
-        rxn.gene_reaction_rule = _stored_gene_rule(row["grRules"])
+        rxn.gpr = _stored_gpr(row["grRules"])
     if _has_value(row.get("lb")) or _has_value(row.get("ub")):
         # Both bounds are set together, so that a pair that is valid as a whole is accepted
         # whatever the stored bounds are.
@@ -965,7 +894,7 @@ def _build_new_reactions(
         )
         rxn.add_metabolites(stoich)
         if _has_value(row.get("grRules")):
-            rxn.gene_reaction_rule = _stored_gene_rule(row["grRules"])
+            rxn.gpr = _stored_gpr(row["grRules"])
         if _has_value(row.get("subSystems")):
             rxn.subsystem = subsystem_to_str(row["subSystems"])
         if _has_value(row.get("eccodes")):

src/raven_toolbox/gapfilling/topological.py

@@ -15,6 +15,7 @@
 """
 from __future__ import annotations
 
+import warnings
 from collections import defaultdict, deque
 from collections.abc import Iterable
 from dataclasses import dataclass, field
@@ -147,40 +148,59 @@ def analyse_topology(
         Universal database model(s) containing candidate reactions.
         Metabolite identifiers must match those in *model*.
     seeds:
-        Metabolite IDs available from the medium. Default: metabolites
-        involved in exchange reactions (one-metabolite reactions) where
-        ``lower_bound < 0`` (uptake is allowed).
+        Metabolite IDs available from the medium. Default: the metabolites of
+        every exchange reaction whose bounds let it supply them, as RAVEN's
+        ``gapFillTopological`` picks them: a reaction with only substrates
+        (``met -->``) and ``lower_bound < 0``, or one with only products
+        (``--> met``) and ``upper_bound > 0``. An exchange reaction may have
+        more than one metabolite. A :class:`UserWarning` is issued when there
+        is none.
     targets:
-        Metabolite IDs that should be reachable. Default: all metabolites
-        in *model*.
+        Metabolite IDs that should be reachable. Default: every metabolite of
+        the first objective reaction (in model order), substrates and products,
+        as in RAVEN's ``gapFillTopological``.
     verbose:
         Print a summary of the results.
 
     Returns
     -------
     TopologicalAnalysisResult
+
+    Raises
+    ------
+    ValueError
+        If *targets* is not given and the model has no objective reaction.
     """
     templates = _as_models(templates)
 
     # ---- Seeds ----
     if seeds is None:
         seed_ids: set[str] = set()
         for rxn in model.reactions:
-            mets = list(rxn.metabolites)
-            if len(mets) == 1 and rxn.lower_bound < 0:
-                seed_ids.add(mets[0].id)
-        if not seed_ids and verbose:
-            print(
-                "analyse_topology: no uptake exchange reactions found — "
-                "supply seeds explicitly via the 'seeds' argument."
+            coeffs = list(rxn.metabolites.values())
+            no_products = all(c < 0 for c in coeffs)
+            no_substrates = all(c > 0 for c in coeffs)
+            if (no_products and rxn.lower_bound < 0) or (no_substrates and rxn.upper_bound > 0):
+                seed_ids.update(m.id for m in rxn.metabolites)
+        if not seed_ids:
+            warnings.warn(
+                "analyse_topology: no uptake exchange reactions found. "
+                "Provide seeds manually via the 'seeds' argument.",
+                stacklevel=2,
             )
     else:
         seed_ids = set(seeds)
 
     # ---- Targets ----
     all_model_met_ids = {m.id for m in model.metabolites}
     if targets is None:
-        target_ids = all_model_met_ids
+        objective = [r for r in model.reactions if r.objective_coefficient != 0]
+        if not objective:
+            raise ValueError(
+                "analyse_topology: no objective function defined. "
+                "Provide targets manually via the 'targets' argument."
+            )
+        target_ids = {m.id for m, c in objective[0].metabolites.items() if c != 0}
     else:
         target_ids = set(targets) & all_model_met_ids