diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d8082da..d949202 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/src/raven_toolbox/init/genes.py b/src/raven_toolbox/init/genes.py index cbee6fd..794514f 100644 --- a/src/raven_toolbox/init/genes.py +++ b/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 diff --git a/src/raven_toolbox/io/git.py b/src/raven_toolbox/io/git.py index 560f19b..15c83fb 100644 --- a/src/raven_toolbox/io/git.py +++ b/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_``), 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: diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 2be733a..92b9e17 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/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 diff --git a/src/raven_toolbox/tasks/check.py b/src/raven_toolbox/tasks/check.py index 51807b4..05eb7ea 100644 --- a/src/raven_toolbox/tasks/check.py +++ b/src/raven_toolbox/tasks/check.py @@ -79,9 +79,17 @@ 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]) @@ -89,6 +97,14 @@ def touch(mid: str) -> list[float]: 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 diff --git a/tests/test_init_genes.py b/tests/test_init_genes.py index 0a9bd09..1ec0590 100644 --- a/tests/test_init_genes.py +++ b/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 == [] diff --git a/tests/test_io_git.py b/tests/test_io_git.py index 2fdcb67..8e59cdb 100644 --- a/tests/test_io_git.py +++ b/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" diff --git a/tests/test_io_yaml.py b/tests/test_io_yaml.py index 54e4f74..c71ef17 100644 --- a/tests/test_io_yaml.py +++ b/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) diff --git a/tests/test_tasks.py b/tests/test_tasks.py index 394fcda..8c83732 100644 --- a/tests/test_tasks.py +++ b/tests/test_tasks.py @@ -1,4 +1,6 @@ """Tests for metabolic tasks (Phase 4a): parse_task_list + check_tasks.""" +import re + import cobra import pytest @@ -169,6 +171,28 @@ def test_unknown_metabolite_reported(model): assert not res.passed and "unknown metabolite" in res.error +def test_in_and_out_with_nonzero_lower_bounds_raises(model): + task = Task("21", "protons", inputs=[("A[c]", 1, 1)], outputs=[("A[c]", 1, 1)]) + with pytest.raises(ValueError, match=r'The IN LB and OUT LB in "\[21\] protons" cannot be ' + "nonzero for the same metabolite"): + check_tasks(model, [task]) + + +def test_in_and_out_with_one_zero_lower_bound_runs(model): + task = Task("io", inputs=[("A[c]", 0, 1000)], outputs=[("A[c]", 1, 1000)]) + (res,) = check_tasks(model, [task]) + assert res.error is None + + +@pytest.mark.parametrize("side", ["inputs", "outputs"]) +def test_metabolite_listed_twice_on_one_side_raises(model, side): + task = Task("dup", "twice", **{side: [("A[c]", 0, 10), ("a[C]", 0, 20)]}) + word = "input(s)" if side == "inputs" else "output(s)" + with pytest.raises(ValueError, match=rf'some {re.escape(word)} in "\[dup\] twice" are ' + "defined more than one time"): + check_tasks(model, [task]) + + def test_open_exchange_is_closed_so_task_controls_io(model): # An open sink would let B leave for free, but check_tasks closes all boundary # reactions first. The task defines input A but no output for B, so B must still