fix: four writer, task and gene-pruning problems from #178

.github/workflows/ci.yml

@@ -22,7 +22,7 @@ jobs:
           python-version: "3.12"
           cache: pip
           cache-dependency-path: pyproject.toml
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       # Pinned so the lint result is reproducible (unpinned ruff silently
       # adds rules over time). Keep in sync with the ``dev`` extra in pyproject.
       - run: pip install ruff==0.15.15
@@ -38,7 +38,7 @@ jobs:
           python-version: "3.12"
           cache: pip
           cache-dependency-path: pyproject.toml
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       # Installs runtime deps + the dev extra (mypy). Third-party stubs are
       # absent for cobra/optlang/scipy; ignore_missing_imports handles that.
       - run: pip install -e ".[dev]"
@@ -65,7 +65,7 @@ jobs:
           python-version: ${{ matrix.python }}
           cache: pip
           cache-dependency-path: pyproject.toml
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       # ``-e .[dev,excel]`` so every optional extra is exercised.
       # Gurobi is not installable on free runners; the relevant tests
       # skip themselves when ``optlang.gurobi_interface`` cannot import.
@@ -98,7 +98,7 @@ jobs:
           python-version: "3.12"
           cache: pip
           cache-dependency-path: pyproject.toml
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       # gurobipy ships a size-limited licence (2000 variables / 2000
       # constraints). That is useless at genome scale but ample for the small
       # fixtures here, and it gives the cross-solver check a second solver to

src/raven_toolbox/init/genes.py

@@ -19,14 +19,18 @@
 from raven_toolbox.utils.gpr import resolve_aggregators
 
 
-def _prune(node, scores, iso, cplx) -> tuple[str | None, float | None]:
-    """Return (pruned GPR string, aggregate score) for an AST node, or (None, None)."""
+def _prune(node, scores, iso, cplx, ties: list) -> tuple[str | None, float | None]:
+    """Return (pruned GPR string, aggregate score) for an AST node, or (None, None).
+
+    Appends to ``ties`` each time an OR node with only negative alternatives keeps
+    one of several alternatives that share the best score.
+    """
     if isinstance(node, ast.Name):
         return node.id, scores.get(node.id)  # None = unscored (NaN: never removed)
     if not isinstance(node, ast.BoolOp):
         return None, None
 
-    pruned = [_prune(v, scores, iso, cplx) for v in node.values]
+    pruned = [_prune(v, scores, iso, cplx, ties) for v in node.values]
     children: list[tuple[str, float | None]] = [(s, sc) for s, sc in pruned if s is not None]
 
     if isinstance(node.op, ast.And):  # complex: keep every subunit, prune nested ORs
@@ -35,6 +39,8 @@ def _prune(node, scores, iso, cplx) -> tuple[str | None, float | None]:
         kept = [(s, sc) for s, sc in children if sc is None or sc >= 0]
         if not kept:  # all negative → keep the least-negative (every sc is non-None here)
             kept = [max(children, key=lambda c: c[1] if c[1] is not None else float("-inf"))]
+            if sum(sc == kept[0][1] for _, sc in children) > 1:
+                ties.append(node)
 
     parts = [s for s, _ in kept]
     score_vals = [sc for _, sc in kept if sc is not None]
@@ -51,17 +57,35 @@ def remove_low_score_genes(
     *,
     isozyme_scoring: str = "max",
     complex_scoring: str = "min",
+    tied_reactions: list[str] | None = None,
 ) -> tuple[cobra.Model, list[str]]:
     """Remove negative-scoring genes from GPRs (RAVEN ``removeLowScoreGenes``).
 
-    ``gene_scores`` maps gene id → score; genes absent from it are treated as unscored
-    (never removed). Returns ``(new_model, removed_gene_ids)`` — genes dropped from
-    *every* rule they were in (and thus from the model). ``isozyme_scoring`` /
-    ``complex_scoring`` aggregate alternative/subunit scores (``max``/``min`` default).
+    Genes joined by OR are alternatives: the negative-scoring ones are removed, and when
+    all alternatives are negative the least-negative one is kept. Genes joined by AND
+    are subunits of a complex and are not removed individually.
+
+    Parameters
+    ----------
+    model : cobra.Model
+        The model to prune. It is not changed.
+    gene_scores : Mapping[str, float]
+        Gene id to score. Genes absent from it are unscored and never removed.
+    isozyme_scoring, complex_scoring : str
+        How the scores of OR alternatives and of complex subunits are aggregated
+        (``"max"`` and ``"min"`` by default).
+    tied_reactions : list of str, optional
+        When a list is given, the ids of the reactions in which an OR rule with only
+        negative alternatives was cut to one of several alternatives with the same
+        best score are appended to it, once per reaction. In such a rule the first of
+        the tied alternatives in the rule is kept, so the result depends on the order
+        of the rule; RAVEN picks one of them at random.
 
-    When all isozyme alternatives are negative the least-negative one is kept
-    **deterministically** (first on a tie), unlike RAVEN's random tie-break — same
-    quality, reproducible.
+    Returns
+    -------
+    tuple of (cobra.Model, list of str)
+        The pruned copy of the model, and the sorted ids of the genes removed from
+        every rule they were in, and thus from the model.
     """
     iso, cplx = resolve_aggregators(isozyme_scoring, complex_scoring)
 
@@ -70,7 +94,10 @@ def remove_low_score_genes(
         body = rxn.gpr.body
         if body is None or not rxn.genes:
             continue
-        pruned, _ = _prune(body, gene_scores, iso, cplx)
+        ties: list = []
+        pruned, _ = _prune(body, gene_scores, iso, cplx, ties)
+        if ties and tied_reactions is not None:
+            tied_reactions.append(rxn.id)
         if pruned is not None:
             rxn.gene_reaction_rule = pruned
 

src/raven_toolbox/io/git.py

@@ -18,6 +18,7 @@
 from pathlib import Path
 
 import cobra
+import libsbml
 from cobra.core import Group
 
 from raven_toolbox.io.excel import _equation, export_to_excel
@@ -127,6 +128,23 @@ def _for_sbml(model: cobra.Model) -> cobra.Model:
     return out
 
 
+def _write_sbml(model: cobra.Model, path: Path) -> None:
+    """Write ``model`` as SBML with each gene product's fbc:label set to the gene name.
+
+    cobra's writer sets fbc:label to the gene product id (``G_<id>``), while RAVEN's
+    exportModel writes the gene short name there, and RAVEN's importModel reads
+    geneShortNames from fbc:label before fbc:name. cobra writes the gene name as
+    fbc:name, or the id when the gene has no name, so copying fbc:name into fbc:label
+    keeps the id as label for unnamed genes. cobra reads gene ids from fbc:id and
+    names from fbc:name, so the label does not change what cobra reads back.
+    """
+    cobra.io.write_sbml_model(_for_sbml(model), str(path))
+    doc = libsbml.readSBMLFromFile(str(path))
+    for gene_product in doc.getModel().getPlugin("fbc").getListOfGeneProducts():
+        gene_product.setLabel(gene_product.getName())
+    libsbml.writeSBMLToFile(doc, str(path))
+
+
 def _write_txt(model: cobra.Model, path: Path) -> None:
     """Single-file, human-readable reaction table (RAVEN exportForGit txt)."""
     with open(path, "w", encoding="utf-8") as fh:
@@ -198,7 +216,7 @@ def target(fmt: str) -> Path:
     if "yml" in formats:
         write_yaml_model(model, target("yml"))
     if "xml" in formats:
-        cobra.io.write_sbml_model(_for_sbml(model), str(target("xml")))
+        _write_sbml(model, target("xml"))
     if "mat" in formats:
         write_matlab_model(model, target("mat"), varname=varname)
     if "xlsx" in formats:

src/raven_toolbox/io/yaml.py

@@ -531,6 +531,7 @@ def _coerce_floats(obj):
 
 _QUOTE_RESOLVER_PATTERN = re.compile(
     r"^(true|True|TRUE|false|False|FALSE)$"
+    r"|^(?i:y|n|yes|no|on|off)$"
     r"|^[-+]?[0-9][0-9_]*\.[0-9_]*([eE][-+]?[0-9]+)?$"
     r"|^[-+]?[0-9][0-9_]*[eE][-+]?[0-9]+$"
     r"|^[-+]?\.[0-9_]+([eE][-+][0-9]+)?$"
@@ -557,13 +558,15 @@ def _is_ws(ch: str) -> bool:
 def _needs_quote(s: str) -> bool:
     """Whether this scalar needs quoting when written as YAML.
 
-    Ports writeYAMLmodel.m's needsQuote so both writers quote exactly the
-    same set of values, for the same two reasons: a plain (bare) reading
-    would resolve to a non-string YAML type (bool/int/float/null/
-    timestamp — neither writer emits explicit type tags), or the text
+    Ports writeYAMLmodel.m's needsQuote, for the same two reasons: a plain
+    (bare) reading would resolve to a non-string YAML type (bool/int/float/
+    null/timestamp, as neither writer emits explicit type tags), or the text
     itself is not valid as a plain block scalar (leading indicator
     character, a mid-string ": "/" #", leading/trailing whitespace, or an
-    embedded line break).
+    embedded line break). Besides the YAML 1.2 booleans it quotes the YAML
+    1.1 ones (y, n, yes, no, on, off, in any case), which a YAML 1.1 loader
+    such as PyYAML would otherwise read as booleans: ``name: NO`` (nitric
+    oxide) is written ``name: "NO"``.
     """
     if not s:
         return True

src/raven_toolbox/tasks/check.py

@@ -79,16 +79,32 @@ def _metabolite_bounds(
 
     Bulk tokens (ALLMETS / ALLMETSIN) are applied before specific metabolites, as
     RAVEN does. Returns the bounds and a list of unresolved tokens (→ task error).
+
+    Raises
+    ------
+    ValueError
+        If a metabolite is listed more than once among the inputs or among the
+        outputs, or if it is both an input and an output with a nonzero lower bound
+        on both sides. The messages are those of RAVEN's ``checkTasks``.
     """
     bounds: dict[str, list[float]] = {}
     missing: list[str] = []
+    label = f'"[{task.id}] {task.description}"'
 
     def touch(mid: str) -> list[float]:
         return bounds.setdefault(mid, [0.0, 0.0])
 
     for entries, is_input in ((task.inputs, True), (task.outputs, False)):
         bulk = [(t, lb, ub) for (t, lb, ub) in entries if _classify(t)[0] != "met"]
         specific = [(t, lb, ub) for (t, lb, ub) in entries if _classify(t)[0] == "met"]
+        # As in RAVEN, only the tokens that match a model metabolite are compared.
+        listed = [_classify(t)[1] for t, _, _ in specific]
+        listed = [arg for arg in listed if arg in name_to_ids]
+        if len(listed) != len(set(listed)):
+            side = "input(s)" if is_input else "output(s)"
+            raise ValueError(
+                f"The constraints on some {side} in {label} are defined more than one time"
+            )
         for token, lb, ub in bulk + specific:
             kind, arg = _classify(token)
             if kind == "all":
@@ -107,6 +123,11 @@ def touch(mid: str) -> list[float]:
                     if kind == "met":
                         b[1] = -lb
                 else:
+                    if kind == "met" and lb > 0 and b[1] < 0:
+                        raise ValueError(
+                            f"The IN LB and OUT LB in {label} cannot be nonzero for the "
+                            "same metabolite"
+                        )
                     b[1] = ub  # allow net production up to UB
                     if kind == "met" and lb > 0:
                         b[0] = lb  # require at least LB produced

tests/test_init_genes.py

@@ -69,3 +69,43 @@ def test_removed_genes_reported_and_pruned():
     out, removed = remove_low_score_genes(_model("G1 or G2"), {"G1": -1.0, "G2": 1.0})
     assert removed == ["G1"]
     assert "G1" not in {g.id for g in out.genes}
+
+
+def test_tie_among_negative_isozymes_reported():
+    """All four isozymes negative with the same score: the first is kept and the
+    reaction is reported."""
+    tied: list[str] = []
+    scores = {"CYP2C9": -1.0, "CYP2C8": -1.0, "CYP3A4": -1.0, "CYP26C1": -1.0}
+    out, removed = remove_low_score_genes(
+        _model("CYP2C9 or CYP2C8 or CYP3A4 or CYP26C1"), scores, tied_reactions=tied
+    )
+    assert out.reactions.R.gene_reaction_rule == "CYP2C9"
+    assert removed == ["CYP26C1", "CYP2C8", "CYP3A4"]
+    assert tied == ["R"]
+
+
+def test_tie_in_nested_or_reported_once():
+    tied: list[str] = []
+    scores = {"G1": -1.0, "G2": -1.0, "G3": -2.0, "G4": -2.0, "G5": 1.0}
+    out, _ = remove_low_score_genes(
+        _model("((G1 or G2) and G5) or ((G3 or G4) and G5)"), scores, tied_reactions=tied
+    )
+    assert out.reactions.R.gene_reaction_rule == _norm("G1 and G5")
+    assert tied == ["R"]
+
+
+def test_distinct_negative_scores_not_reported():
+    tied: list[str] = []
+    out, _ = remove_low_score_genes(
+        _model("G1 or G2"), {"G1": -2.0, "G2": -1.0}, tied_reactions=tied
+    )
+    assert out.reactions.R.gene_reaction_rule == "G2"
+    assert tied == []
+
+
+def test_tie_with_a_positive_alternative_not_reported():
+    tied: list[str] = []
+    remove_low_score_genes(
+        _model("G1 or G2 or G3"), {"G1": -1.0, "G2": -1.0, "G3": 1.0}, tied_reactions=tied
+    )
+    assert tied == []

tests/test_io_git.py

@@ -1,5 +1,6 @@
 """Tests for raven_toolbox.io.git (exportForGit port)."""
 import cobra
+import libsbml
 import pytest
 from cobra.core import Group
 
@@ -174,3 +175,35 @@ def test_sbml_notes_keep_the_model_note():
     m = _sbml_model()
     set_model_metadata(m, note="a generic human cell")
     assert _for_sbml(m).notes["note"] == "a generic human cell"
+
+
+def _consistency(path):
+    doc = libsbml.readSBMLFromFile(str(path))
+    doc.checkConsistency()
+    return sorted(
+        (doc.getError(i).getSeverityAsString(), doc.getError(i).getErrorId())
+        for i in range(doc.getNumErrors())
+    )
+
+
+def test_sbml_gene_label_is_the_gene_name(model, tmp_path):
+    """RAVEN's importModel reads geneShortNames from fbc:label before fbc:name."""
+    model.reactions.R1.gene_reaction_rule = "ENSG1 or ENSG2"
+    model.genes.ENSG1.name = "DPM1"
+    root = export_for_git(model, tmp_path, formats=("xml",))
+    path = root / "xml" / "model.xml"
+
+    doc = libsbml.readSBMLFromFile(str(path))
+    products = doc.getModel().getPlugin("fbc").getListOfGeneProducts()
+    labels = {gp.getId(): (gp.getName(), gp.getLabel()) for gp in products}
+    assert labels == {"G_ENSG1": ("DPM1", "DPM1"), "G_ENSG2": ("G_ENSG2", "G_ENSG2")}
+
+    plain = tmp_path / "plain.xml"
+    cobra.io.write_sbml_model(_for_sbml(model), str(plain))
+    assert _consistency(path) == _consistency(plain)
+
+    back, plain_back = cobra.io.read_sbml_model(str(path)), cobra.io.read_sbml_model(str(plain))
+    assert sorted((g.id, g.name) for g in back.genes) == sorted(
+        (g.id, g.name) for g in plain_back.genes
+    )
+    assert back.reactions.R1.gene_reaction_rule == "ENSG1 or ENSG2"

tests/test_io_yaml.py

@@ -390,3 +390,45 @@ def test_sbo_term_written_as_a_list_is_read_as_a_string(tmp_path):
     )
     model = read_yaml_model(path)
     assert model.reactions.get_by_id("r1").annotation["sbo"] == "SBO:0000627"
+
+
+def test_write_quotes_yaml_1_1_booleans(tmp_path):
+    """Nitric oxide (name and formula NO), formula Y and compartment n are
+    YAML 1.1 booleans, which PyYAML reads as True or False when unquoted."""
+    m = cobra.Model("t")
+    m.compartments = {"n": "nucleus"}
+    no = cobra.Metabolite("no_n", name="NO", formula="NO", compartment="n")
+    y = cobra.Metabolite("y_n", name="on", formula="Y", compartment="n")
+    m.add_metabolites([no, y])
+    r = cobra.Reaction("R1", name="yes")
+    r.add_metabolites({no: -1, y: 1})
+    m.add_reactions([r])
+
+    out = tmp_path / "m.yml"
+    write_yaml_model(m, out)
+    text = out.read_text()
+    for line in ('name: "NO"', 'formula: "NO"', 'formula: "Y"', 'compartment: "n"',
+                 'name: "on"', 'name: "yes"'):
+        assert line in text
+    assert "compartment: n\n" not in text
+
+    back = read_yaml_model(out)
+    assert back.metabolites.no_n.name == "NO"
+    assert back.metabolites.no_n.formula == "NO"
+    assert back.metabolites.y_n.formula == "Y"
+    assert back.metabolites.y_n.compartment == "n"
+    assert back.reactions.R1.name == "yes"
+
+
+@pytest.mark.parametrize("value", ["y", "N", "Yes", "nO", "ON", "Off"])
+def test_needs_quote_yaml_1_1_booleans_any_case(value):
+    from raven_toolbox.io.yaml import _needs_quote
+
+    assert _needs_quote(value)
+
+
+@pytest.mark.parametrize("value", ["NOS1", "Yn", "onion", "ny", "off-target"])
+def test_needs_quote_leaves_longer_words_plain(value):
+    from raven_toolbox.io.yaml import _needs_quote
+
+    assert not _needs_quote(value)