read_yaml_model reads empty enzymes maps and quotes inside values; fseof follows ecFSEOF

CHANGELOG.md

@@ -6,6 +6,14 @@ Milestones in the raven-toolbox port. For function-level status see
 
 ## Unreleased
 
+* `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.
+* `analysis.fseof` gained `equality`: `True` fixes the target flux at each level as GECKO's `ecFSEOF` does, `False` (the default) keeps the lower bound of RAVEN's `FSEOF`.
+* `io.read_yaml_model` reads an ec reaction written as `- enzymes: !!omap` with nothing under the tag, the form GECKO and RAVEN write for a reaction without enzymes, as a reaction with no enzymes, as `readYAMLmodel` does. It raised `AttributeError: 'ScalarNode' object has no attribute 'flow_style'`.
+* `io.read_yaml_model` reads a scalar written as `""text""` (invalid YAML, written by older RAVEN and GECKO for a note that carried its own quotes) as `"text"`, as `readYAMLmodel` does. Files that parse as YAML are read unchanged.
+* `io.read_yaml_model` reads a double-quoted scalar with an unescaped quote inside, such as `- notes: "a "b" c"`, as `a "b" c`: it drops the outer pair of quotes and unescapes `\\` and `\"`, as `readYAMLmodel` does. This repair runs only when the strict parse fails.
+* `io.read_yaml_model` reads a single-quoted scalar with an apostrophe inside, such as `- name: 'it's'`, as `it's`: it drops the outer pair of quotes and reads each `''` as `'`, as `readYAMLmodel` does. This repair runs only when the strict parse fails; a comment after the closing quote is not handled.
 * `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux
   capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as
   `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its

src/raven_toolbox/analysis/fseof.py

@@ -26,6 +26,12 @@
 from cobra.flux_analysis import pfba
 from scipy.stats import linregress
 
+# pFBA keeps the objective at this fraction of the optimum found at the same level, as
+# RAVEN's solveLP does when it minimises the flux sum (``solution.f*0.999999``). A solver
+# reports its optimum to within its feasibility tolerance, so a second solve that demands
+# the full value can be infeasible on that solver (HiGHS at the top level of ecYeastGEM).
+_PFBA_FRACTION = 0.999999
+
 
 @dataclass
 class FSEOFResult:
@@ -81,11 +87,13 @@ def fseof(
     correlation_threshold: float = 0.9,
     flux_eps: float = 1e-6,
     min_target: float | None = None,
+    equality: bool = False,
 ) -> FSEOFResult:
     """Run FSEOF for over-production of ``target_rxn``'s product.
 
     Enforces target flux from ``max_fraction/n_steps`` up to ``max_fraction`` of the
-    theoretical maximum in ``n_steps`` steps, maximising growth (``biomass_rxn`` or the
+    theoretical maximum (the target flux of the pFBA solution that maximises it) in
+    ``n_steps`` steps, maximising growth (``biomass_rxn`` or the
     model's current objective) with pFBA at each step. Returns an :class:`FSEOFResult`.
 
     Parameters
@@ -102,24 +110,39 @@ def fseof(
         it is ``knockout`` when the slope is negative and the last absolute flux is at most
         1e-8, ``knockdown`` for a negative slope otherwise.
 
+        Each level is solved in two steps as in RAVEN: the objective is maximised, then the
+        flux sum is minimised with the objective held at 0.999999 of its optimum at that
+        level. A level at which the first step is infeasible ends the scan.
+
+        The top of the scan is ``max_fraction`` of the target flux in the pFBA solution that
+        maximises the target, as in RAVEN's ``FSEOF`` and GECKO's ``ecFSEOF``.
+
         Differences from ``FSEOF.m``: the scan stops at the first infeasible enforced level,
         so the last absolute flux used for the knockout label is that of the last feasible
-        level; ``targets`` is sorted by ``|slope| * |r|`` and not grouped by label; and
+        level; ``targets`` is sorted by ``|slope| * |r|`` (equal scores in model order) and is
+        neither grouped by label nor split into reaction and transport targets; and
         ``flux_eps`` only gates the check of the theoretical maximum.
     min_target:
         Floor of the enforced-flux range, evenly spaced up to the scan's ceiling
         (``max_fraction`` of the theoretical maximum) in ``n_steps`` steps.
         Defaults to ``None``, which floors the range at
         ``max_fraction of the theoretical maximum / n_steps`` instead. Pass e.g.
         the target reaction's flux at growth-optimum to scan from there.
+    equality:
+        How the enforced flux is imposed. ``False`` (the default) sets it as the lower bound
+        of the target reaction, as RAVEN's ``FSEOF`` does; ``True`` fixes the target flux at
+        the level (lower and upper bound), as GECKO's ``ecFSEOF`` does.
     """
     with model:  # find the theoretical maximum target flux
         model.objective = target_rxn
         target_opt = model.slim_optimize()
     # slim_optimize returns NaN on an infeasible model; np.isfinite catches that too.
     if target_opt is None or not np.isfinite(target_opt) or target_opt <= flux_eps:
         raise ValueError(f"{target_rxn!r} cannot carry positive flux; nothing to scan.")
-    target_max = target_opt * max_fraction
+    with model:  # the maximum is read from the pFBA solution, as RAVEN's solveLP(model, 1) gives it
+        model.objective = target_rxn
+        target_pfba = pfba(model, fraction_of_optimum=_PFBA_FRACTION).fluxes[target_rxn]
+    target_max = float(target_pfba) * max_fraction
     if min_target is None:
         levels = [target_max * (i + 1) / n_steps for i in range(n_steps)]
     else:
@@ -131,9 +154,13 @@ def fseof(
         with model:
             if biomass_rxn is not None:
                 model.objective = biomass_rxn
-            model.reactions.get_by_id(target_rxn).lower_bound = level
+            target = model.reactions.get_by_id(target_rxn)
+            if equality:
+                target.bounds = (level, level)
+            else:
+                target.lower_bound = level
             try:
-                columns[level] = pfba(model).fluxes
+                columns[level] = pfba(model, fraction_of_optimum=_PFBA_FRACTION).fluxes
             except OptimizationError:
                 break  # enforced flux became infeasible: stop scanning
             enforced.append(level)
@@ -187,4 +214,5 @@ def _classify(model, scan, enforced, corr_threshold) -> pd.DataFrame:
         "reaction", "name", "subsystem", "gene_reaction_rule", "genes",
         "target_type", "slope", "correlation", "initial_flux", "final_flux", "score",
     ])
-    return table.sort_values("score", ascending=False, ignore_index=True)
+    # A stable sort keeps reactions with equal scores in model order, as MATLAB's sort does.
+    return table.sort_values("score", ascending=False, ignore_index=True, kind="stable")

src/raven_toolbox/io/yaml.py

@@ -49,8 +49,10 @@
 
 import cobra
 from cobra.io.dict import model_from_dict, model_to_dict
-from cobra.io.yaml import yaml as _cobra_yaml  # ruamel round-trip YAML (handles !!omap)
 from ruamel.yaml import YAML
+from ruamel.yaml.constructor import RoundTripConstructor
+from ruamel.yaml.error import YAMLError
+from ruamel.yaml.nodes import ScalarNode
 from ruamel.yaml.scalarstring import DoubleQuotedScalarString
 
 from raven_toolbox.io.ec_data import (
@@ -59,11 +61,10 @@
     ec_data_to_yaml_sections,
 )
 
-# A dedicated instance for *writing* only (reading still goes through
-# _cobra_yaml, whose parsing is unaffected by any of this). Mutating
-# _cobra_yaml's own width would also change cobra.io.save_yaml_model for
-# the rest of the process, so this module gets its own ruamel round-trip
-# instance instead. `width` is set far past any real line length so
+# A dedicated instance for *writing* only (reading uses _read_yaml below).
+# Mutating cobra's own YAML instance would change its width for
+# cobra.io.save_yaml_model in the rest of the process, so this module gets its own
+# ruamel round-trip instance instead. `width` is set far past any real line length so
 # ruamel never folds a long scalar onto a continuation line, matching
 # writeYAMLmodel.m. `indent` is pinned explicitly (to the values ruamel
 # already used implicitly, verified byte-identical) rather than left for
@@ -83,6 +84,112 @@
 _write_yaml.indent(mapping=2, sequence=2, offset=0)
 
 
+class _OmapConstructor(RoundTripConstructor):
+    """Round-trip constructor that reads a content-less ``!!omap`` as an empty map.
+
+    GECKO and RAVEN write an ec reaction without enzymes as ``- enzymes: !!omap``
+    with nothing under the tag. YAML parses that as a null scalar carrying the omap
+    tag, which ruamel's own constructor rejects; RAVEN's ``readYAMLmodel.m`` reads
+    it as an empty enzyme list.
+    """
+
+    def construct_yaml_omap(self, node):
+        if isinstance(node, ScalarNode) and not node.value:
+            yield OrderedDict()
+            return
+        yield from super().construct_yaml_omap(node)
+
+
+_OmapConstructor.add_constructor("tag:yaml.org,2002:omap", _OmapConstructor.construct_yaml_omap)
+
+# A dedicated instance for *reading*, so the constructor above does not change
+# cobra.io.load_yaml_model for the rest of the process.
+_read_yaml = YAML(typ="rt")
+_read_yaml.Constructor = _OmapConstructor
+
+# A double-quoted scalar that holds a double quote no backslash escapes, such as
+# `""text""` (the shape older RAVEN and GECKO writers gave to a note that already
+# carried its own quotes) or `"a "b" c"`. It is not valid YAML. `readYAMLmodel.m` strips
+# the outer pair of quotes of the line and unescapes `\\` and `\"` in what is left, so the
+# values are `"text"` and `a "b" c`.
+_QUOTED_VALUE = re.compile(
+    r"""^(?P<head>\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)"(?P<body>.*)"\s*$"""
+)
+
+
+# A single-quoted scalar that holds an apostrophe no second apostrophe pairs with, such as
+# `'it's'`. It is not valid YAML. `readYAMLmodel.m` strips the outer pair of quotes of the line
+# and turns each `''` in what is left into `'`, so the value is `it's`.
+_SINGLE_QUOTED_VALUE = re.compile(
+    r"""^(?P<head>\s*-\s+|\s*(?:-\s+)?[^\s:\-][^:]*:\s+)'(?P<body>.*)'\s*$"""
+)
+
+
+def _has_lone_apostrophe(body: str) -> bool:
+    """Whether ``body``, the text between the outer quotes, has a ``'`` that is not half of a ``''``."""
+    i = 0
+    while i < len(body):
+        if body[i] != "'":
+            i += 1
+        elif body[i : i + 2] == "''":
+            i += 2
+        else:
+            return True
+    return False
+
+
+def _has_unescaped_quote(body: str) -> bool:
+    """Whether ``body``, the text between the outer quotes, has a ``"`` that no backslash escapes."""
+    i = 0
+    while i < len(body):
+        if body[i] == "\\":
+            i += 2
+        elif body[i] == '"':
+            return True
+        else:
+            i += 1
+    return False
+
+
+def _fix_doubled_quotes(text: str) -> str:
+    """Rewrite each quoted scalar with an unpaired interior quote as a valid single-quoted scalar.
+
+    A double-quoted scalar with an unescaped ``"`` takes the text between the outer quotes with
+    escaped backslashes and quotes unescaped as its value. A single-quoted scalar with an
+    apostrophe that no second apostrophe pairs with takes the text between the outer quotes
+    with each ``''`` read as ``'``. Both are the values ``readYAMLmodel.m`` reads; quoted
+    scalars without such a quote are left alone.
+    """
+
+    def fix(match: re.Match) -> str:
+        body = match["body"]
+        if not _has_unescaped_quote(body):
+            return match[0]
+        value = re.sub(r'\\(["\\])', r"\1", body)
+        return f"{match['head']}'{value.replace(chr(39), chr(39) * 2)}'"
+
+    def fix_single(match: re.Match) -> str:
+        body = match["body"]
+        if not _has_lone_apostrophe(body):
+            return match[0]
+        return f"{match['head']}'{body.replace(chr(39) * 2, chr(39)).replace(chr(39), chr(39) * 2)}'"
+
+    return "\n".join(
+        _SINGLE_QUOTED_VALUE.sub(fix_single, _QUOTED_VALUE.sub(fix, line)) for line in text.split("\n")
+    )
+
+
+def _load_yaml_text(text: str):
+    """Parse ``text`` strictly; when that fails, retry with unpaired quotes inside quoted scalars repaired."""
+    try:
+        return _read_yaml.load(text)
+    except YAMLError:
+        repaired = _fix_doubled_quotes(text)
+        if repaired == text:
+            raise
+        return _read_yaml.load(repaired)
+
+
 def _open_text(path: str | Path, mode: str):
     """Open ``path`` as a text handle, transparently gzipping when it ends ``.gz``."""
     if str(path).endswith(".gz"):
@@ -166,7 +273,7 @@ def read_yaml_model(path: str | Path) -> cobra.Model:
     latter is merged into one mapping before parsing.
     """
     with _open_text(path, "r") as handle:
-        raw = _to_plain(_cobra_yaml.load(handle))
+        raw = _to_plain(_load_yaml_text(handle.read()))
 
     if isinstance(raw, list):
         # Very old RAVEN files: a sequence of one-key mappings instead of

tests/test_analysis_fseof.py

@@ -171,3 +171,94 @@ def test_knockout_threshold_is_1e_minus_8():
     table = _classify(m, scan, enforced, corr_threshold=0.5).set_index("reaction")
     assert table.loc["r_down", "target_type"] == "knockdown"
     assert table.loc["r_out", "target_type"] == "knockout"
+
+
+def test_pfba_holds_objective_at_matlab_fraction(model, monkeypatch):
+    # RAVEN's solveLP minimises the flux sum with the objective fixed at 0.999999 of its optimum,
+    # so a solver whose optimum overshoots by its tolerance still has a feasible second step.
+    import importlib
+
+    module = importlib.import_module("raven_toolbox.analysis.fseof")
+
+    fractions = []
+    real = module.pfba
+
+    def spy(m, **kwargs):
+        fractions.append(kwargs.get("fraction_of_optimum"))
+        return real(m, **kwargs)
+
+    monkeypatch.setattr(module, "pfba", spy)
+    res = fseof(model, "EX_P", n_steps=4)
+    assert len(fractions) >= len(res.enforced) >= 2
+    assert set(fractions) == {0.999999}
+
+
+def test_level_with_infeasible_maximisation_ends_the_scan(model):
+    # The theoretical maximum of EX_P is 10; levels above it cannot be enforced and are not
+    # reported, however the second pFBA step is tolerated.
+    res = fseof(model, "EX_P", n_steps=4, max_fraction=2.0)
+    assert res.enforced == pytest.approx([5.0, 10.0])
+
+
+def test_equality_fixes_the_target_flux_at_each_level(model):
+    res = fseof(model, "EX_P", n_steps=4, equality=True)
+    assert res.scan.loc["EX_P"].to_numpy() == pytest.approx(res.enforced)
+    # the default lower bound lets a level exceed the enforced flux when that suits growth;
+    # here the product branch is the only sink of the shared intermediate at the top level.
+    lower = fseof(model, "EX_P", n_steps=4)
+    assert (lower.scan.loc["EX_P"].to_numpy() >= lower.enforced).all()
+
+
+def test_equality_differs_from_lower_bound_when_the_target_is_coupled_to_growth():
+    # S -> B + P is the only route to biomass, so growth alone exports 10 of P; a lower bound
+    # below that leaves the flux at 10, while an equality holds it at the level.
+    m = cobra.Model("m")
+    S, P, B = (cobra.Metabolite(x, compartment="c") for x in ("S", "P", "B"))
+    m.add_metabolites([S, P, B])
+    sup = cobra.Reaction("sup", lower_bound=0, upper_bound=10)
+    sup.add_metabolites({S: 1})
+    v1 = cobra.Reaction("v1", lower_bound=0, upper_bound=1000)
+    v1.add_metabolites({S: -1, B: 1, P: 1})
+    ex_p = cobra.Reaction("EX_P", lower_bound=0, upper_bound=1000)
+    ex_p.add_metabolites({P: -1})
+    ex_b = cobra.Reaction("EX_B", lower_bound=0, upper_bound=1000)
+    ex_b.add_metabolites({B: -1})
+    m.add_reactions([sup, v1, ex_p, ex_b])
+    m.objective = "EX_B"
+    lower = fseof(m, "EX_P", n_steps=4, min_target=2.0)
+    fixed = fseof(m, "EX_P", n_steps=4, min_target=2.0, equality=True)
+    assert lower.scan.loc["EX_P"].to_numpy() == pytest.approx(10.0)
+    assert fixed.scan.loc["EX_P"].to_numpy() == pytest.approx(fixed.enforced)
+
+
+def test_top_level_is_a_fraction_of_the_pfba_maximum(model):
+    # pFBA holds the target at 0.999999 of its maximum while it minimises the flux sum,
+    # so the top level is 0.9 of that flux and not of the plain maximum.
+    from cobra.flux_analysis import pfba
+
+    with model:
+        model.objective = "EX_P"
+        plain = model.slim_optimize()
+        in_pfba = pfba(model, fraction_of_optimum=0.999999).fluxes["EX_P"]
+    assert in_pfba < plain
+    res = fseof(model, "EX_P", n_steps=4)
+    assert res.enforced[-1] == pytest.approx(0.9 * in_pfba, rel=1e-8)
+    assert res.enforced[-1] != pytest.approx(0.9 * plain, rel=1e-8)
+
+
+def test_equal_scores_keep_model_order():
+    # Reactions of three score classes are interleaved in the model; within a class the
+    # scores tie exactly, and the table lists each class in model order, highest score first.
+    import numpy as np
+    import pandas as pd
+
+    from raven_toolbox.analysis.fseof import _classify
+
+    m = cobra.Model("m")
+    ids = [f"r{i:03d}" for i in range(600)]
+    m.add_reactions([cobra.Reaction(rid) for rid in ids])
+    enforced = np.array([1.0, 2.0, 3.0, 4.0])
+    scan = pd.DataFrame([enforced * (1 + i % 3) for i in range(600)], index=ids, columns=enforced)
+    table = _classify(m, scan, enforced, 0.9)
+    expected = [ids[i] for c in (2, 1, 0) for i in range(600) if i % 3 == c]
+    assert list(table["reaction"]) == expected