From 0922ae4b147e1dbb6b23543dc08b0801f413abff Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 11:17:25 +0100 Subject: [PATCH 01/35] ftinit: resolve_ties takes each step's fluxes from a least-total-flux solution After the id-rank phase, a further phase holds the on/off state and minimises the summed absolute flux, so the direction and magnitude at which the next step forces a reaction do not depend on the solver's vertex. --- CHANGELOG.md | 3 ++ src/raven_toolbox/init/ftinit.py | 76 +++++++++++++++++++++++++++++--- tests/test_init_ftinit.py | 39 ++++++++++++++++ 3 files changed, 113 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b629dd2..bfd078e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,9 @@ 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. ## 3.0.0b1 — 2026-09-07 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index bd67126..8d3b062 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -187,7 +187,8 @@ def run_ftinit( exhaust ``time_limit``; when that happens the (still-adopted) incumbent is an unproven tie-break and a warning is raised. The tie-break is skipped when the primary solve itself stops at ``time_limit``, since there is then no proven optimum to hold. - See :func:`_resolve_ties`. + After it, the reported ``fluxes`` are a least-total-flux solution for the selected on + set (:func:`_min_abs_flux`). See :func:`_resolve_ties`. ``metabolomics`` (already in **this model's** id space — :func:`ftinit` resolves a caller-facing metabolite-name list into this once) maps a detected-metabolite label @@ -507,6 +508,10 @@ def _read_solution(): elif resolve_ties and indicators and _resolve_ties(opt, prob, obj_expr, indicators, primary_obj, time_limit): on, fluxes = _read_solution() + # Like phases 2a/2b, only when there is a removable (negative-score) reaction. + if any(score < 0 for _, score in indicators.values()) \ + and _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): + fluxes = _read_solution()[1] kept = free_or_essential | on deleted = [r.id for r in model.reactions if r.id not in kept] @@ -556,9 +561,9 @@ def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit) -> bool: indicators are continuous and pinned near 1 by the score objective), so the two phases run over those: their count and id-sum are integers, provable with a cheap absolute gap below 1 rather than the near-full proof a tiny relative gap would need. This resolves - the removable-reaction choices — the genuine seed-fragile degeneracy; a residual - flux-distribution degeneracy (which only feeds the next step's small ``ess_force`` - clamp) is left to the solver. + the removable-reaction choices, the genuine seed-fragile degeneracy. The flux + distribution, which the staged schedule reads for the next step's forced directions and + ``ess_force`` clamp, is pinned afterwards by :func:`_min_abs_flux`. Caveat, measured on genome-scale Human-GEM (not visible at toy/unit-test scale, where both phases prove in milliseconds): the phases can themselves exhaust ``time_limit``, @@ -629,6 +634,64 @@ def _phase(objective, label: str) -> bool: return ok +def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: + """Phase 2c of ``resolve_ties``: the least total absolute flux for the tie-broken on set. + + Phases 2a/2b fix which reactions are on but leave the flux distribution free, and the + staged schedule reads that distribution: the next step forces a reversible reaction in + the direction of its flux and at ``min(0.99·|flux|, force_on)``. This phase holds the + on/off state of every indicator (a binary at its value; a continuous indicator that is + on keeps at least its value, one that is off at most its value) and minimises + ``Σ|v|`` over all reactions, so the flux read-out comes from a minimal-flux solution + rather than from whichever vertex the solver happened to reach. Directions of + reversible scored reactions stay free, so the phase is still a MILP; its objective is + proven to an absolute gap of 1e-6. + + Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), + ``False`` otherwise (the caller keeps the phase-2b fluxes). It never fails the + extraction; a phase ending at the time limit is adopted with a warning. + """ + # Read every value before changing a bound: a bound change discards the solution. + values = [(ind, ind.primal or 0.0) for ind, _ in indicators.values()] + for ind, val in values: + if ind.type == "binary": + ind.set_bounds(round(val), round(val)) + elif val >= 0.5: + ind.lb = min(max(val, 0.0), 1.0) + else: + ind.ub = min(max(val, 0.0), 1.0) + abs_terms, extra = [], [] + for rid, terms in flux_terms.items(): + if len(terms) == 1 and terms[0][0].lb is not None and terms[0][0].lb >= 0: + abs_terms.append(mul([Real(terms[0][1]), terms[0][0]])) + elif len(terms) == 1 and terms[0][0].ub is not None and terms[0][0].ub <= 0: + abs_terms.append(mul([Real(-terms[0][1]), terms[0][0]])) + else: + net = add([mul([Real(sign), var]) for var, sign in terms]) + t = prob.Variable(f"_abs_{rid}", lb=0.0) + extra += [t, prob.Constraint(t - net, lb=0.0, name=f"_absp_{rid}"), + prob.Constraint(t + net, lb=0.0, name=f"_absn_{rid}")] + abs_terms.append(t) + opt.add(extra) + opt.objective = prob.Objective(add(abs_terms), direction="min") + if _is_gurobi(prob): + opt.problem.Params.MIPGap = 0.0 + opt.problem.Params.MIPGapAbs = 1e-6 + if time_limit is not None: + opt.configuration.timeout = max(1, int(time_limit * _TIE_BREAK_TIME_FRACTION)) + opt.optimize() + _dbg(f"[ftinit] tie-break phase2c-minflux: status={opt.status} obj={opt.objective.value}") + if opt.status == "time_limit": + warnings.warn( + "ftINIT minimal-flux phase (resolve_ties) stopped at the time limit: the flux " + "distribution, and hence the next step's forced directions, is an unproven " + "incumbent. Raise time_limit for a proven minimal-flux solution.", + stacklevel=3, + ) + return (opt.status in ("optimal", "feasible", "suboptimal", "time_limit") + and _has_solution(opt)) + + def _nudge_scores(rxn_scores: Mapping[str, float]) -> dict[str, float]: """Push tiny reaction scores off zero (RAVEN ``ftINIT.m:160-161``). @@ -860,7 +923,10 @@ def ftinit( * ``resolve_ties`` adds a lexicographic phase 2 that selects the sparsest (then lowest-id) optimum, so the degenerate choice is pinned rather than left to the - solver — applied both to each extraction step and to the task gap-fill. + solver — applied both to each extraction step and to the task gap-fill. In each + extraction step it then takes the fluxes from the least-total-flux solution for + that on set, so the directions and forced magnitudes the next step reads from + them do not depend on the solver's vertex either. * ``prove_abs_gap`` replaces the loose relative-gap escalation with a single solve per step proven to this fixed *absolute* gap, so the kept-set objective is optimal at any scale — not an arbitrary within-gap incumbent. **Try 1.0**: on Human-GEM/DLD1 diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 581f5db..42d4c0f 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -220,6 +220,45 @@ def test_resolve_ties_safe_on_unique_optimum(): assert res.objective == pytest.approx(8.0, abs=1e-6) +def _direction_model(): + """A positive reversible R3 (m <=> q) that can be on in either direction. + + Run forward it needs m from a via the free R1; run backward it imports q and feeds the + essential E. Both directions reach the same score, so the flux of R3 (and the direction + the next ftINIT step would force it in) is a tie that only the flux distribution + resolves. The free loop L1/L2 can carry any flux. R2 is a negative-score dead end. + """ + m = cobra.Model("direction") + a, mm, p, q, n, z = (cobra.Metabolite(x, name=x, compartment="s") + for x in ("a", "m", "p", "q", "n", "z")) + rxns = { + "R1": ({a: -1, mm: 1}, 0), "R2": ({a: -1, z: 1}, 0), "R3": ({mm: -1, q: 1}, -1000), + "E": ({mm: -1, p: 1}, 0), "L1": ({mm: -1, n: 1}, 0), "L2": ({n: -1, mm: 1}, 0), + "EX_a": ({a: -1}, -1000), "EX_p": ({p: -1}, -1000), "EX_q": ({q: -1}, -1000), + } + for rid, (mets, lb) in rxns.items(): + r = cobra.Reaction(rid, lower_bound=lb, upper_bound=1000) + m.add_reactions([r]) + r.add_metabolites(mets) + return m + + +def test_resolve_ties_takes_fluxes_from_a_minimal_flux_solution(): + """resolve_ties reports the least-total-flux distribution for the tie-broken on set. + + The backward R3 route (import q, feed E) carries 0.4 in total, the forward one (a via + R1 into both R3 and E) 0.8, and the free loop L1/L2 carries nothing. + """ + res = run_ftinit(_direction_model(), {"R2": -1.0, "R3": 1.0}, essential_rxns=["E"], + resolve_ties=True) + assert "R3" in res.on_reactions and "R2" not in res.kept_reactions + expected = {"R1": 0.0, "R2": 0.0, "R3": -0.1, "E": 0.1, "L1": 0.0, "L2": 0.0, + "EX_a": 0.0, "EX_p": 0.1, "EX_q": -0.1} + for rid, v in expected.items(): + assert res.fluxes[rid] == pytest.approx(v, abs=1e-6), rid + assert res.objective == pytest.approx(1.0, abs=1e-6) + + # --------------------------------------------------------------------------- # # seed as an explicit parameter, and the unproven-incumbent warning. # --------------------------------------------------------------------------- # From b0fc04c4214a67170f337cd86f10753e002f3a2c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 18:16:42 +0100 Subject: [PATCH 02/35] ftinit: the resolve_ties minimal-flux phase stops at a 1e-6 relative or absolute gap and has the full time_limit --- CHANGELOG.md | 2 ++ src/raven_toolbox/init/ftinit.py | 9 +++++---- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bfd078e..67a14ea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,8 @@ Milestones in the raven-toolbox port. For function-level status see * `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. ## 3.0.0b1 — 2026-09-07 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 8d3b062..fad3ab0 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -644,8 +644,9 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: on keeps at least its value, one that is off at most its value) and minimises ``Σ|v|`` over all reactions, so the flux read-out comes from a minimal-flux solution rather than from whichever vertex the solver happened to reach. Directions of - reversible scored reactions stay free, so the phase is still a MILP; its objective is - proven to an absolute gap of 1e-6. + reversible scored reactions stay free, so the phase is still a MILP; it stops at a + relative or absolute gap of 1e-6, with the full ``time_limit`` as its budget (as RAVEN's + phase 2c). Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), ``False`` otherwise (the caller keeps the phase-2b fluxes). It never fails the @@ -675,10 +676,10 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: opt.add(extra) opt.objective = prob.Objective(add(abs_terms), direction="min") if _is_gurobi(prob): - opt.problem.Params.MIPGap = 0.0 + opt.problem.Params.MIPGap = 1e-6 opt.problem.Params.MIPGapAbs = 1e-6 if time_limit is not None: - opt.configuration.timeout = max(1, int(time_limit * _TIE_BREAK_TIME_FRACTION)) + opt.configuration.timeout = max(1, int(time_limit)) opt.optimize() _dbg(f"[ftinit] tie-break phase2c-minflux: status={opt.status} obj={opt.objective.value}") if opt.status == "time_limit": From d699d60e09772442494784838883a565bef566f6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 18:23:46 +0100 Subject: [PATCH 03/35] ftinit: resolve_ties fixes equal-total-flux directions by id rank and ranks a merged reaction by its smallest original id --- CHANGELOG.md | 6 ++ src/raven_toolbox/init/ftinit.py | 112 ++++++++++++++++++++++++++----- tests/test_init_ftinit.py | 38 ++++++++++- 3 files changed, 138 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 67a14ea..62373e5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,12 @@ Milestones in the raven-toolbox port. For function-level status see 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. ## 3.0.0b1 — 2026-09-07 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index fad3ab0..4f4defe 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -142,6 +142,7 @@ def run_ftinit( time_limit: float | None = None, prove_abs_gap: float | None = None, resolve_ties: bool = False, + rank_keys: Mapping[str, str] | None = None, metabolomics: Mapping[str, Iterable[str]] | None = None, prod_weight: float = 0.5, seed: int = _EXTRACT_SEED, @@ -188,7 +189,13 @@ def run_ftinit( unproven tie-break and a warning is raised. The tie-break is skipped when the primary solve itself stops at ``time_limit``, since there is then no proven optimum to hold. After it, the reported ``fluxes`` are a least-total-flux solution for the selected on - set (:func:`_min_abs_flux`). See :func:`_resolve_ties`. + set (:func:`_min_abs_flux`) and, among those, the one with the least id-rank-weighted + backward flux (:func:`_min_backward_flux`). See :func:`_resolve_ties`. + + ``rank_keys`` maps a reaction id to the key the id-rank phases of ``resolve_ties`` + sort it by (default: the id itself). :func:`ftinit` passes the smallest original + reaction id of each merge group, so the order does not depend on which member names + the merged reaction. ``metabolomics`` (already in **this model's** id space — :func:`ftinit` resolves a caller-facing metabolite-name list into this once) maps a detected-metabolite label @@ -503,15 +510,23 @@ def _read_solution(): # It is skipped when the primary solve stopped at the time limit: phase 2 would pin # an unproven incumbent's objective, so it cannot select among true optima, and at # genome scale it then only spends its own time budget without finding a solution. + keys = dict(rank_keys or {}) + + def rank_key(rid: str) -> str: + return keys.get(rid, rid) + if resolve_ties and indicators and _status == "time_limit": _dbg("[ftinit] tie-break skipped: the primary solve stopped at the time limit") elif resolve_ties and indicators and _resolve_ties(opt, prob, obj_expr, indicators, - primary_obj, time_limit): + primary_obj, time_limit, rank_key): on, fluxes = _read_solution() # Like phases 2a/2b, only when there is a removable (negative-score) reaction. - if any(score < 0 for _, score in indicators.values()) \ - and _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): + held = (_min_abs_flux(opt, prob, indicators, flux_terms, time_limit) + if any(score < 0 for _, score in indicators.values()) else None) + if held is not None: fluxes = _read_solution()[1] + if _min_backward_flux(opt, prob, flux_terms, *held, rank_key, time_limit): + fluxes = _read_solution()[1] kept = free_or_essential | on deleted = [r.id for r in model.reactions if r.id not in kept] @@ -541,7 +556,7 @@ def _has_solution(opt) -> bool: return True -def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit) -> bool: +def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit, rank_key) -> bool: """Pin ftINIT's degenerate optimum to a single, reproducible solution (on ``opt``). The ftINIT MILP is highly degenerate: many reaction subsets reach the same score @@ -551,8 +566,9 @@ def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit) -> bool: This runs a lexicographic phase 2, holding the score objective at its optimum with a floor constraint: first minimise the count of kept removable reactions (the sparsest - optimum), then, among the sparsest, minimise their summed id rank (prefer lower ids). - The result is a stable, near-unique optimum independent of seed/solver version *when + optimum), then, among the sparsest, minimise their summed id rank (prefer lower ids; + ``rank_key`` maps a reaction id to the key it is ranked by). The result is a stable, + near-unique optimum independent of seed/solver version *when both phases converge* — see the caveat below. ``opt`` is left holding the phase-2 solution, which the caller reads for the "on" set and fluxes; the reported objective stays the phase-1 ``primary`` value. @@ -618,7 +634,7 @@ def _phase(objective, label: str) -> bool: unproven.append("phase2a-objective-not-finite") # Phase 2b — among the sparsest, prefer lower reaction ids (deterministic tie-break). opt.add(prob.Constraint(count, ub=kmin + 0.5, name="_canon_count_cap")) - ranks = {rid: i for i, rid in enumerate(sorted(binaries))} + ranks = {rid: i for i, rid in enumerate(sorted(binaries, key=rank_key))} idsum = add([mul([Real(float(1 + ranks[rid])), ind]) for rid, ind in binaries.items()]) ok = _phase(prob.Objective(idsum, direction="min"), "phase2b-idrank") if ok and unproven: @@ -634,7 +650,7 @@ def _phase(objective, label: str) -> bool: return ok -def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: +def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): """Phase 2c of ``resolve_ties``: the least total absolute flux for the tie-broken on set. Phases 2a/2b fix which reactions are on but leave the flux distribution free, and the @@ -648,9 +664,11 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: relative or absolute gap of 1e-6, with the full ``time_limit`` as its budget (as RAVEN's phase 2c). - Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), - ``False`` otherwise (the caller keeps the phase-2b fluxes). It never fails the - extraction; a phase ending at the time limit is adopted with a warning. + Returns ``(abs_expr, total, free_sign)`` if it produced a usable solution (the caller + then reads its fluxes): the ``Σ|v|`` expression, its value and the ids of the + reactions whose flux can take either sign. Returns ``None`` otherwise (the caller keeps + the phase-2b fluxes). It never fails the extraction; a phase ending at the time limit + is adopted with a warning. """ # Read every value before changing a bound: a bound change discards the solution. values = [(ind, ind.primal or 0.0) for ind, _ in indicators.values()] @@ -661,7 +679,7 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: ind.lb = min(max(val, 0.0), 1.0) else: ind.ub = min(max(val, 0.0), 1.0) - abs_terms, extra = [], [] + abs_terms, extra, free_sign = [], [], [] for rid, terms in flux_terms.items(): if len(terms) == 1 and terms[0][0].lb is not None and terms[0][0].lb >= 0: abs_terms.append(mul([Real(terms[0][1]), terms[0][0]])) @@ -673,8 +691,10 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: extra += [t, prob.Constraint(t - net, lb=0.0, name=f"_absp_{rid}"), prob.Constraint(t + net, lb=0.0, name=f"_absn_{rid}")] abs_terms.append(t) + free_sign.append(rid) opt.add(extra) - opt.objective = prob.Objective(add(abs_terms), direction="min") + abs_expr = add(abs_terms) + opt.objective = prob.Objective(abs_expr, direction="min") if _is_gurobi(prob): opt.problem.Params.MIPGap = 1e-6 opt.problem.Params.MIPGapAbs = 1e-6 @@ -689,6 +709,56 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit) -> bool: "incumbent. Raise time_limit for a proven minimal-flux solution.", stacklevel=3, ) + if opt.status in ("optimal", "feasible", "suboptimal", "time_limit") and _has_solution(opt): + return abs_expr, float(opt.objective.value), free_sign + return None + + +def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_key, + time_limit) -> bool: + """Phase 2d of ``resolve_ties``: fix the directions left open by the minimal total flux. + + Phase 2c can leave the direction of a flux free when two directions carry the same total + flux, for example a pair of reversible reactions that close a loop in either direction. + This phase keeps the on/off state of phase 2c, holds ``Σ|v|`` at most at its phase-2c + value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -v)`` over the reactions + whose flux can take either sign. The weight is the reaction's id rank from the end: + sorted by ``rank_key``, the first of ``n`` reactions gets ``n`` and the last ``1``, so a + reaction with a lower id prefers its forward direction. It stops at a relative or + absolute gap of 1e-6, with the full ``time_limit`` as its budget. + + Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), + ``False`` otherwise (the caller keeps the phase-2c fluxes). It never fails the + extraction; a phase ending at the time limit is adopted with a warning. + """ + if not free_sign: + return False + tol = max(abs(total) * 1e-7, 1e-6) + n = len(free_sign) + weight = {rid: float(n - i) for i, rid in enumerate(sorted(free_sign, key=rank_key))} + extra, bwd_terms = [prob.Constraint(abs_expr, ub=total + tol, name="_minflux_cap")], [] + for rid in free_sign: + net = add([mul([Real(sign), var]) for var, sign in flux_terms[rid]]) + b = prob.Variable(f"_bwd_{rid}", lb=0.0) + extra += [b, prob.Constraint(b + net, lb=0.0, name=f"_bwdc_{rid}")] + bwd_terms.append(mul([Real(weight[rid]), b])) + opt.add(extra) + opt.objective = prob.Objective(add(bwd_terms), direction="min") + if _is_gurobi(prob): + opt.problem.Params.MIPGap = 1e-6 + opt.problem.Params.MIPGapAbs = 1e-6 + if time_limit is not None: + opt.configuration.timeout = max(1, int(time_limit)) + opt.optimize() + _dbg(f"[ftinit] tie-break phase2d-direction: status={opt.status} obj={opt.objective.value}") + if opt.status == "time_limit": + warnings.warn( + "ftINIT direction phase (resolve_ties) stopped at the time limit: the " + "directions of fluxes that carry the same total flux either way, and hence the " + "next step's forced directions, are an unproven incumbent. Raise time_limit " + "for a proven solution.", + stacklevel=3, + ) return (opt.status in ("optimal", "feasible", "suboptimal", "time_limit") and _has_solution(opt)) @@ -771,7 +841,7 @@ def _metabolomics_producers(prep, metabolomics: Iterable[str]) -> dict[str, set[ def _solve_step( min_model, scores, step, *, essential, directions, ess_force, force_on, big_m, mip_gap, mip_gap_abs, time_limit, prove_abs_gap=None, resolve_ties=False, - metabolomics=None, prod_weight=0.5, + rank_keys=None, metabolomics=None, prod_weight=0.5, seed=_EXTRACT_SEED, threads=_EXTRACT_THREADS, feas_tol=_EXTRACT_TOL, opt_tol=_EXTRACT_TOL, int_feas_tol=_EXTRACT_TOL, presolve=_EXTRACT_PRESOLVE, @@ -792,7 +862,7 @@ def _run(mg, mga, tl, prove=None): rem_pos_rev=step.pos_rev_off, ignore_mets=step.mets_to_ignore, force_on=force_on, force_on_ess=force_on, big_m=big_m, mip_gap=mg, mip_gap_abs=mga, time_limit=tl, - prove_abs_gap=prove, resolve_ties=resolve_ties, + prove_abs_gap=prove, resolve_ties=resolve_ties, rank_keys=rank_keys, metabolomics=metabolomics, prod_weight=prod_weight, seed=seed, threads=threads, feas_tol=feas_tol, opt_tol=opt_tol, int_feas_tol=int_feas_tol, presolve=presolve, @@ -992,6 +1062,14 @@ def ftinit( # it is never forced above what it last carried (ftINIT.m:172,248) — this applies to # the permanent (prep) essentials too, not only reactions turned on by a prior step. flux_of: dict[str, float] = {r.id: force_on for r in min_model.reactions} + # resolve_ties ranks a merged reaction by its group's smallest original id, which does + # not depend on which member names the merged reaction. + smallest: dict[int, str] = {} + for rid, gid in zip(prep.orig_rxn_ids, prep.group_ids, strict=True): + if gid and (gid not in smallest or rid < smallest[gid]): + smallest[gid] = rid + rank_keys = {r.id: smallest[group_of[r.id]] for r in min_model.reactions + if group_of.get(r.id, 0)} n_steps = len(steps) for i, step in enumerate(steps): if verbose: @@ -1011,7 +1089,7 @@ def ftinit( min_model, scores, step, essential=essential, directions=directions, ess_force=ess_force, force_on=force_on, big_m=big_m, mip_gap=mip_gap, mip_gap_abs=mip_gap_abs, time_limit=time_limit, - prove_abs_gap=prove_abs_gap, resolve_ties=resolve_ties, + prove_abs_gap=prove_abs_gap, resolve_ties=resolve_ties, rank_keys=rank_keys, metabolomics=met_producers, prod_weight=prod_weight, seed=seed, threads=threads, feas_tol=feas_tol, opt_tol=opt_tol, int_feas_tol=int_feas_tol, presolve=presolve, diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 42d4c0f..4585e25 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -259,6 +259,42 @@ def test_resolve_ties_takes_fluxes_from_a_minimal_flux_solution(): assert res.objective == pytest.approx(1.0, abs=1e-6) +def _two_route_model(): + """The essential E consumes y, which the reversible R1 (y <=> x) or R2 (y <=> z) can + supply when run backward from the imports IN_x or IN_z. + + Both routes carry the same total flux, so the least total flux leaves open which of R1 + and R2 runs backward. N is a negative-score dead end. + """ + m = cobra.Model("two_route") + x, y, z, p, d = (cobra.Metabolite(k, name=k, compartment="s") for k in "xyzpd") + rxns = { + "IN_x": ({x: 1}, 0), "IN_z": ({z: 1}, 0), "R1": ({y: -1, x: 1}, -1000), + "R2": ({y: -1, z: 1}, -1000), "E": ({y: -1, p: 1}, 0), "OUT_p": ({p: -1}, 0), + "N": ({x: -1, d: 1}, 0), + } + for rid, (mets, lb) in rxns.items(): + r = cobra.Reaction(rid, lower_bound=lb, upper_bound=1000) + m.add_reactions([r]) + r.add_metabolites(mets) + return m + + +def test_resolve_ties_lower_ranked_reaction_keeps_its_forward_direction(): + """Between equal-total-flux directions, resolve_ties runs the higher-ranked one backward. + + Ranked by id, R1 comes before R2, so R2 carries the backward flux; with ``rank_keys`` + ranking R1 last, R1 does. + """ + def solve(rank_keys=None): + res = run_ftinit(_two_route_model(), {"N": -1.0}, essential_rxns=["E"], + resolve_ties=True, rank_keys=rank_keys) + return res.fluxes["R1"], res.fluxes["R2"] + + assert solve() == pytest.approx((0.0, -0.1), abs=1e-6) + assert solve({"R1": "R9"}) == pytest.approx((-0.1, 0.0), abs=1e-6) + + # --------------------------------------------------------------------------- # # seed as an explicit parameter, and the unproven-incumbent warning. # --------------------------------------------------------------------------- # @@ -353,7 +389,7 @@ def optimize(self): try: with pytest.warns(UserWarning, match="tie resolution did not converge"): ok = ftinit_mod._resolve_ties(stub, interface, interface.Variable("objx"), - {"R1": (ind, -1.0)}, -1.0, 1.0) + {"R1": (ind, -1.0)}, -1.0, 1.0, str) finally: ftinit_mod._has_solution = monkey assert ok # the incumbent is still adopted, just no longer silently From d3ab2792a24c167ee89831e3ed5bccacca07c6a1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 18:26:46 +0100 Subject: [PATCH 04/35] ftinit: unpack the minimal-flux phase result before the direction phase --- src/raven_toolbox/init/ftinit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 4f4defe..2d0dd62 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -525,7 +525,9 @@ def rank_key(rid: str) -> str: if any(score < 0 for _, score in indicators.values()) else None) if held is not None: fluxes = _read_solution()[1] - if _min_backward_flux(opt, prob, flux_terms, *held, rank_key, time_limit): + abs_expr, total, free_sign = held + if _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, + rank_key, time_limit): fluxes = _read_solution()[1] kept = free_or_essential | on From ae8cf6e47841dc2e61a6ece61c278bce7e6f8790 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 20:03:38 +0100 Subject: [PATCH 05/35] ftinit: the resolve_ties direction phase starts from the minimal-flux solution --- src/raven_toolbox/init/ftinit.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 2d0dd62..6c3b497 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -726,8 +726,10 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -v)`` over the reactions whose flux can take either sign. The weight is the reaction's id rank from the end: sorted by ``rank_key``, the first of ``n`` reactions gets ``n`` and the last ``1``, so a - reaction with a lower id prefers its forward direction. It stops at a relative or - absolute gap of 1e-6, with the full ``time_limit`` as its budget. + reaction with a lower id prefers its forward direction. On Gurobi the phase-2c solution + is its MIP start, as the tight cap on ``Σ|v|`` leaves the solver's heuristics little + room to find a first solution. It stops at a relative or absolute gap of 1e-6, with + the full ``time_limit`` as its budget. Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), ``False`` otherwise (the caller keeps the phase-2c fluxes). It never fails the @@ -735,6 +737,8 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k """ if not free_sign: return False + # Read the phase-2c solution before the model changes: a change discards it. + start = {var.name: var.primal or 0.0 for var in opt.variables} if _is_gurobi(prob) else {} tol = max(abs(total) * 1e-7, 1e-6) n = len(free_sign) weight = {rid: float(n - i) for i, rid in enumerate(sorted(free_sign, key=rank_key))} @@ -749,6 +753,14 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k if _is_gurobi(prob): opt.problem.Params.MIPGap = 1e-6 opt.problem.Params.MIPGapAbs = 1e-6 + opt.update() + grb = opt.problem + grb.update() + for name, val in start.items(): + grb.getVarByName(name).Start = val + for rid in free_sign: + net = sum(sign * start[var.name] for var, sign in flux_terms[rid]) + grb.getVarByName(f"_bwd_{rid}").Start = max(0.0, -net) if time_limit is not None: opt.configuration.timeout = max(1, int(time_limit)) opt.optimize() From bd660c7f884932d8f64b0b25354af6487ec365a1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:11:34 +0100 Subject: [PATCH 06/35] merge_compartments: keep the id of the copy in the merged compartment or the one that sorts first The merged metabolite took the first member's id with its compartment suffix rewritten (ACCOA_c became ACCOA_s). It now keeps the id of the copy already in the merged compartment, if any, otherwise the id that sorts first, as RAVEN's mergeCompartments does. --- CHANGELOG.md | 4 ++ .../manipulation/compartments.py | 47 +++++------------ tests/test_manipulation_compartments.py | 50 ++++++++++++++++--- 3 files changed, 59 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b629dd2..caebf08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `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`). * `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 diff --git a/src/raven_toolbox/manipulation/compartments.py b/src/raven_toolbox/manipulation/compartments.py index 1911fe3..6cd405e 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -47,13 +47,13 @@ def merge_compartments( → ``glc__D``). Models that key the same species to a *different id per compartment* — e.g. yeast-GEM's ``s_####`` ids — need ``base_metabolite=lambda m: m.name`` (or an annotation key) so the compartments actually unify; pass the same callable to - :func:`~raven_toolbox.localization.assign_compartments`. The key **only decides grouping** — the - merged metabolite keeps an *existing* id (the copy already in ``merged_id``, or an arbitrary member - adapted to ``merged_id`` when none is there yet), never one minted from the key. So a name-based key - like ``"1-acyl-sn-glycerol 3-phosphate (16:0)"`` unifies the copies while the merged id stays the - model's own, already solver/SBML-safe identifier — no character sanitising. (cobra rejects an id - containing whitespace when a model is built or loaded, so the source ids are solver-safe to begin - with; this transform only ever inherits them.) + :func:`~raven_toolbox.localization.assign_compartments`. The key **only decides grouping**. The + merged metabolite keeps the id, name, formula, charge and notes of one of its copies, never an id + minted from the key: the copy already in ``merged_id`` if there is one, otherwise the copy whose id + sorts first (``ACCOA_c`` before ``ACCOA_m``); among several copies in ``merged_id`` the one whose id + sorts first. RAVEN's ``mergeCompartments`` applies the same rule. The id can therefore name a + compartment that no longer exists. A name-based key like ``"1-acyl-sn-glycerol 3-phosphate (16:0)"`` + unifies the copies while the merged id stays the model's own, already solver/SBML-safe identifier. Metabolites that share a key collapse into one entity in the merged compartment; their stoichiometric contributions are summed per reaction. Reactions that *become* trivial @@ -74,19 +74,14 @@ def merge_compartments( base = key(m) new_to_old.setdefault(base, []).append(m) - # 2. Build the merged metabolites and rewrite reactions. The merged metabolite inherits an existing - # id rather than one minted from the key, so it stays solver/SBML-safe by construction: keep the - # id of the copy already in the merged compartment (it already reflects it), or, when none is - # there yet, adapt an arbitrary member's id to the merged compartment. A name key thus never - # leaks whitespace/punctuation into a metabolite id. + # 2. Build the merged metabolites and rewrite reactions. Each merged metabolite is a copy of one + # member, chosen as the docstring states. Groups are disjoint and model ids are unique, so the + # chosen ids are unique too. canonical: dict[str, cobra.Metabolite] = {} - used_ids: set[str] = set() for base, mets in new_to_old.items(): - in_target = next((m for m in mets if m.compartment == merged_id), None) - proto = in_target if in_target is not None else mets[0] - mid = proto.id if in_target is not None else f"{_base_id(proto)}_{merged_id}" - mid = _unique_id(mid, used_ids) - used_ids.add(mid) + in_target = [m for m in mets if m.compartment == merged_id] + proto = min(in_target or mets, key=lambda m: m.id) + mid = proto.id new_met = cobra.Metabolite(mid, name=proto.name, compartment=merged_id, formula=proto.formula, charge=proto.charge) new_met.notes = dict(proto.notes or {}) @@ -231,22 +226,6 @@ def copy_to_compartment( # ----------------------------------------------------------------- helpers -def _unique_id(candidate: str, used: set[str]) -> str: - """Return ``candidate`` if free, else ``candidate_2``, ``candidate_3`` ... . - - Used only to break the rare collision between two *derived* merged ids (two groups whose arbitrary - members reduce to the same ``_``). No character sanitising: ``candidate`` is either - an existing metabolite id or built from one, so it is already solver- and SBML-safe (cobra rejects a - whitespace id at model construction, so the source ids are legal to begin with); this transform - never introduces an illegal id.""" - if candidate not in used: - return candidate - n = 2 - while f"{candidate}_{n}" in used: - n += 1 - return f"{candidate}_{n}" - - def _base_id(m: cobra.Metabolite) -> str: """Strip the trailing ``_`` suffix from a metabolite id (if present).""" if m.compartment and m.id.endswith(f"_{m.compartment}"): diff --git a/tests/test_manipulation_compartments.py b/tests/test_manipulation_compartments.py index 269f12e..113d7bc 100644 --- a/tests/test_manipulation_compartments.py +++ b/tests/test_manipulation_compartments.py @@ -67,11 +67,11 @@ def rxn(rid, lb, ub, mets, gpr=None): def test_merge_compartments_collapses_to_one(): """A_c + A_m collapse to one species; B_c + B_m to another; transport A_c↔A_m self-cancels and is - dropped. No copy lived in the (default, synthetic) merged compartment 's', so each merged id is an - existing id adapted to it: A_s, B_s.""" + dropped. No copy lived in the merged compartment 's', so each merged metabolite keeps the id that + sorts first: A_c, B_c.""" m = _two_compartment_model() merged, deleted, dupes = merge_compartments(m) - assert {x.id for x in merged.metabolites} == {"A_s", "B_s"} + assert {x.id for x in merged.metabolites} == {"A_c", "B_c"} # The transport reaction collapsed (A → A) and was deleted. assert "tr_A" in deleted # r_c and r_m are now both A → B; one of them gets deduplicated. @@ -159,8 +159,8 @@ def test_merge_compartments_default_key_cannot_unify_name_keyed_model(): suffix-stripping key leaves every metabolite distinct — nothing collapses. This is the failure the base_metabolite override exists to fix.""" merged, deleted, _ = merge_compartments(_name_keyed_model()) - # nothing unifies (4 distinct species); each existing id is adapted to the merged compartment 's'. - assert {x.id for x in merged.metabolites} == {"s_1_s", "s_2_s", "s_3_s", "s_4_s"} + # nothing unifies (4 distinct species); each metabolite keeps its own id. + assert {x.id for x in merged.metabolites} == {"s_1", "s_2", "s_3", "s_4"} assert deleted == [] # tr_A (s_1 -> s_2) still spans two distinct mets, not collapsed @@ -170,8 +170,8 @@ def test_merge_compartments_base_metabolite_unifies_by_name(): A->B reactions deduplicate — the same outcome as suffix-keyed ids get by default.""" merged, deleted, dupes = merge_compartments( _name_keyed_model(), base_metabolite=lambda m: m.name) - # merged id is inherited from an existing member (adapted to 's'), NOT synthesised from the name key. - assert {x.id for x in merged.metabolites} == {"s_1_s", "s_3_s"} + # merged id is the id of a member, NOT synthesised from the name key. + assert {x.id for x in merged.metabolites} == {"s_1", "s_3"} assert "tr_A" in deleted surviving = {r.id for r in merged.reactions} assert len(surviving & {"r_c", "r_m"}) == 1 @@ -199,7 +199,7 @@ def test_merge_by_name_inherits_legal_ids_despite_punctuated_names(): m.add_reactions([r, tr]) merged, deleted, _ = merge_compartments(m, base_metabolite=lambda mm: mm.name) ids = {x.id for x in merged.metabolites} - assert ids == {"x1_s", "x3_s"} # inherited from x1/x3, adapted to 's' — never the raw name + assert ids == {"x1", "x3"} # inherited from x1/x3, never the raw name assert "tr" in deleted # the two punctuated copies unified, so the transport cancels merged.slim_optimize() # must not raise while building the solver problem @@ -212,6 +212,40 @@ def test_merge_into_existing_compartment_keeps_that_compartments_ids(): assert {x.id for x in merged.metabolites} == {"A_c", "B_c"} # the 'c' copies' ids survive +def _acetyl_coa_model(first: tuple[str, str], second: tuple[str, str]) -> cobra.Model: + """Two copies of acetyl-CoA, given as (id, compartment) in model order, and X_c: + R1 X_c => first copy, R2 second copy => X_c.""" + m = cobra.Model("accoa") + a = cobra.Metabolite(first[0], name="acetyl-CoA", compartment=first[1]) + b = cobra.Metabolite(second[0], name="acetyl-CoA", compartment=second[1]) + x = cobra.Metabolite("X_c", name="X", compartment="c") + m.add_metabolites([a, b, x]) + r1 = cobra.Reaction("R1", lower_bound=0, upper_bound=1000) + r1.add_metabolites({x: -1, a: 1}) + r2 = cobra.Reaction("R2", lower_bound=0, upper_bound=1000) + r2.add_metabolites({b: -1, x: 1}) + m.add_reactions([r1, r2]) + return m + + +def test_merge_compartments_keeps_the_id_that_sorts_first(): + """The merged metabolite keeps the id that sorts first, whatever the order of the copies in + the model (RAVEN's mergeCompartments applies the same rule).""" + for order in ((("ACCOA_m", "m"), ("ACCOA_c", "c")), (("ACCOA_c", "c"), ("ACCOA_m", "m"))): + merged, _, _ = merge_compartments(_acetyl_coa_model(*order), + base_metabolite=lambda mm: mm.name) + assert {x.id for x in merged.metabolites} == {"ACCOA_c", "X_c"} + assert merged.metabolites.ACCOA_c.compartment == "s" + + +def test_merge_compartments_prefers_the_copy_in_the_merged_compartment(): + """A copy already in the merged compartment keeps its id, even when another copy's id sorts + before it.""" + merged, _, _ = merge_compartments(_acetyl_coa_model(("ACCOA_c", "c"), ("ACCOA_z", "s")), + base_metabolite=lambda mm: mm.name) + assert {x.id for x in merged.metabolites} == {"ACCOA_z", "X_c"} + + # ----------------------------------------------------------------- copy_to_compartment def test_copy_to_compartment_basic(): From dcfc3c3d5950f10fc088a6864ec4b3d621ad0198 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:02 +0100 Subject: [PATCH 07/35] merge_linear: number groups by their first reaction Group numbers followed the merge order. They are now 1, 2, ... in the order of the first reaction of each group in the original model, as RAVEN's mergeLinear numbers them. --- CHANGELOG.md | 3 +++ src/raven_toolbox/init/merge.py | 13 ++++++++++--- tests/test_init_merge.py | 20 ++++++++++++++++++++ 3 files changed, 33 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index caebf08..eebb605 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,9 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `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 diff --git a/src/raven_toolbox/init/merge.py b/src/raven_toolbox/init/merge.py index ce87fcf..670e0d4 100644 --- a/src/raven_toolbox/init/merge.py +++ b/src/raven_toolbox/init/merge.py @@ -12,7 +12,8 @@ * ``group_ids`` — one integer per original reaction; ``0`` = not merged, equal non-zero integers = merged into the same combined reaction (which keeps one - member's id). + member's id). The groups are numbered 1, 2, ... in the order of their first + reaction in the original model, as RAVEN's ``mergeLinear`` numbers them. * ``reversed_rxns`` — which originals were flipped (their stored direction negated) when oriented for merging; needed to map fluxes/directions back. @@ -158,8 +159,14 @@ def relabel(rx: _Rxn, grp: int) -> None: break rxns = [rx for rx in rxns if rx.coeffs] - return _build_model(model, rxns), orig_ids, [group_of[o] for o in orig_ids], \ - [reversed_of[o] for o in orig_ids] + # Number the groups in the order of their first original reaction, so the numbers do not + # depend on the order in which the reactions were merged. + number: dict[int, int] = {} + for o in orig_ids: + if group_of[o] and group_of[o] not in number: + number[group_of[o]] = len(number) + 1 + group_ids = [number.get(group_of[o], 0) for o in orig_ids] + return _build_model(model, rxns), orig_ids, group_ids, [reversed_of[o] for o in orig_ids] def _build_model(template: cobra.Model, rxns: list[_Rxn]) -> cobra.Model: diff --git a/tests/test_init_merge.py b/tests/test_init_merge.py index 028a820..29254f0 100644 --- a/tests/test_init_merge.py +++ b/tests/test_init_merge.py @@ -99,6 +99,26 @@ def test_multipass_chain_collapses_to_one_group(): assert {mt.id: c for mt, c in only.metabolites.items()} == {"A": -1.0, "Z": 1.0} +def test_groups_are_numbered_by_their_first_reaction(): + """Group numbers follow the first reaction of each group in the original model, not the + merge order. Here R3+R4 merge before R1+R2's group reaches them, and the merge of two + groups keeps the larger number, so numbering by merge order gives 2, 2, 2, 2, 3, 3.""" + m = cobra.Model("numbering") + mets = {i: cobra.Metabolite(i, name=i, compartment="c") for i in ("x", "a", "m", "b", "y", "z", "c", "w")} + m.add_metabolites(list(mets.values())) + chain = [("R1", "x", "a"), ("R2", "a", "m"), ("R3", "m", "b"), ("R4", "b", "y"), + ("R5", "z", "c"), ("R6", "c", "w")] + for rid, sub, prod in chain: + r = cobra.Reaction(rid, lower_bound=0, upper_bound=1000) + r.add_metabolites({mets[sub]: -1, mets[prod]: 1}) + m.add_reactions([r]) + + reduced, orig_ids, group_ids, _ = merge_linear(m) + assert orig_ids == ["R1", "R2", "R3", "R4", "R5", "R6"] + assert group_ids == [1, 1, 1, 1, 2, 2] + assert sorted(r.id for r in reduced.reactions) == ["R1", "R5"] + + def test_group_scores_zero_handling(): """Genuine-zero score → 0.01; a group cancelling to zero with nonzero members → 0.01.""" reduced, orig_ids, group_ids, _ = merge_linear(make_test_model()) From 3144fdb1e27e9ab6f4eee4a1894ac71789939336 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:43 +0100 Subject: [PATCH 08/35] gpr_to_dnf and is_dnf read a rule string with parse_gr_rule, the grammar of RAVEN's parseGrRule, and raise ValueError for a malformed rule --- CHANGELOG.md | 6 ++ src/raven_toolbox/curation/batch.py | 79 ++------------- src/raven_toolbox/manipulation/expand.py | 18 +++- src/raven_toolbox/utils/__init__.py | 3 +- src/raven_toolbox/utils/gpr.py | 122 ++++++++++++++++++----- tests/test_manipulation_expand.py | 20 ++++ tests/test_utils_gpr.py | 45 ++++++++- 7 files changed, 195 insertions(+), 98 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b629dd2..27e9f3f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* 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. * `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 diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index d12d8bf..d42565b 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -13,6 +13,7 @@ import numpy as np import pandas as pd +from raven_toolbox.utils.gpr import _RULE_TOKEN, parse_gr_rule, rule_tokens from raven_toolbox.utils.ids import next_ids from raven_toolbox.utils.parse import subsystem_to_str @@ -434,75 +435,18 @@ 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]: @@ -526,15 +470,10 @@ def _stored_gene_rule(rule) -> 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: + kind = rule_tokens(word)[0][0] + if word in _COBRA_OPERATORS or kind not in ("and", "or"): continue + replacement = kind pieces.append(text[last:match.start()]) pieces.append(replacement) last = match.end() diff --git a/src/raven_toolbox/manipulation/expand.py b/src/raven_toolbox/manipulation/expand.py index f509ea4..2bd98c1 100644 --- a/src/raven_toolbox/manipulation/expand.py +++ b/src/raven_toolbox/manipulation/expand.py @@ -15,10 +15,12 @@ import cobra from cobra.core.gene import GPR +from raven_toolbox.utils.gpr import parse_gr_rule + __all__ = ["expand_model", "gpr_to_dnf"] -def gpr_to_dnf(gpr: GPR) -> list[list[str]]: +def gpr_to_dnf(gpr: GPR | str | None) -> list[list[str]]: """Convert a GPR to disjunctive normal form (list of AND-clauses). An empty GPR yields an empty list. A single clause (no OR anywhere) @@ -27,7 +29,21 @@ def gpr_to_dnf(gpr: GPR) -> list[list[str]]: Handles distributivity: ``g1 and (g2 or g3)`` becomes ``[[g1, g2], [g1, g3]]``. + + Parameters + ---------- + gpr + A cobra :class:`~cobra.core.gene.GPR`, a grRule string or ``None``. A string is + read with :func:`raven_toolbox.utils.parse_gr_rule`, the grammar of RAVEN's + ``parseGrRule``, as RAVEN's ``grRuleToDNF`` reads it. + + Raises + ------ + ValueError + If a rule string is malformed in that grammar (RAVEN raises ``RAVEN:badGrRule``). """ + if isinstance(gpr, str): + gpr = parse_gr_rule(gpr) if gpr is None or gpr.body is None: return [] return _node_to_dnf(gpr.body) diff --git a/src/raven_toolbox/utils/__init__.py b/src/raven_toolbox/utils/__init__.py index 3919e0a..58a97bb 100644 --- a/src/raven_toolbox/utils/__init__.py +++ b/src/raven_toolbox/utils/__init__.py @@ -2,7 +2,7 @@ curation checks, id generation and sorting.""" from raven_toolbox.utils.balance import ElementalBalance, get_elemental_balance from raven_toolbox.utils.composition import GuessCompositionResult, guess_composition -from raven_toolbox.utils.gpr import GPRIssue, find_non_dnf_grrules, is_dnf +from raven_toolbox.utils.gpr import GPRIssue, find_non_dnf_grrules, is_dnf, parse_gr_rule from raven_toolbox.utils.ids import generate_new_ids from raven_toolbox.utils.sort import sort_identifiers, sort_model, subsystem_order_scores from raven_toolbox.utils.validate import ModelIssue, check_model @@ -18,6 +18,7 @@ "get_elemental_balance", "guess_composition", "is_dnf", + "parse_gr_rule", "sort_identifiers", "sort_model", "subsystem_order_scores", diff --git a/src/raven_toolbox/utils/gpr.py b/src/raven_toolbox/utils/gpr.py index 86e352b..c84f763 100644 --- a/src/raven_toolbox/utils/gpr.py +++ b/src/raven_toolbox/utils/gpr.py @@ -17,9 +17,8 @@ from __future__ import annotations import ast -import logging +import re import statistics -import warnings from collections.abc import Callable from dataclasses import dataclass @@ -88,13 +87,13 @@ def is_dnf(gpr: GPR | str | None) -> bool: ---------- gpr A cobra :class:`~cobra.core.gene.GPR`, a grRule string, or ``None``. - An empty/``None`` rule is trivially DNF. + A string is read with :func:`parse_gr_rule`. An empty/``None`` rule is trivially DNF. Raises ------ ValueError - If a rule string is not valid (unbalanced parentheses, a dangling - operator, two genes without an operator between them, ``()``). + If a rule string is not valid in the grammar of :func:`parse_gr_rule` (unbalanced + parentheses, a dangling operator, two genes without an operator between them, ``()``). Examples -------- @@ -103,33 +102,108 @@ def is_dnf(gpr: GPR | str | None) -> bool: >>> is_dnf("(G1 or G2) and G3") False """ - parsed = _parse_rule(gpr) if isinstance(gpr, str) else gpr + parsed = parse_gr_rule(gpr) if isinstance(gpr, str) else gpr if parsed is None: return True return _is_dnf_node(parsed.body) -def _parse_rule(rule: str) -> GPR | None: - """Parse a grRule string; ``None`` for a blank rule, ``ValueError`` if malformed. +_RULE_TOKEN = re.compile(r"\(|\)|[^\s()]+") +_AND_WORDS = frozenset({"and", "&", "&&"}) +_OR_WORDS = frozenset({"or", "|", "||"}) - ``GPR.from_string`` swallows a syntax error, logs it and returns an empty GPR, which - would make any malformed rule look trivially DNF. RAVEN's ``isDnfGrRule`` raises - ``RAVEN:badGrRule`` instead, so this does too. + +def rule_tokens(rule: str) -> list[tuple[str, str]]: + """Split a gene rule into ``(kind, text)`` pairs; kind is ``(``, ``)``, ``and``, ``or`` or ``gene``. + + A word is everything up to whitespace or a bracket. The operators are ``and``/``or`` in any + case and the symbols ``&``, ``&&``, ``|`` and ``||``, each as a whole word, as RAVEN's + ``parseGrRule`` reads them; every other word is a gene id, so ``G1&G2`` is one gene. + """ + 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 parse_gr_rule(rule: str | None) -> GPR | None: + """Parse a grRule string in the grammar of RAVEN's ``parseGrRule``. + + The grammar, with AND binding tighter than OR:: + + expr := term ( OR term )* + term := factor ( AND factor )* + factor := gene | '(' expr ')' + + AND is ``and`` in any case, ``&`` or ``&&``; OR is ``or`` in any case, ``|`` or ``||``. + Operators and brackets are separated from gene ids by whitespace or a bracket; any other + word is a gene id. ``curate_model_from_tables``/``batch_curate``, :func:`is_dnf` and + :func:`~raven_toolbox.manipulation.gpr_to_dnf` all read rules with this parser, so they + accept and refuse the same rules as RAVEN. + + Parameters + ---------- + rule + A grRule such as ``"(G1 and G2) or G3"``. + + Returns + ------- + GPR or None + The parsed rule; ``None`` for ``None``, an empty or an all-blank rule. + + Raises + ------ + ValueError + If the rule is malformed: an operator without an operand, unbalanced or empty + parentheses, or two genes without an operator between them. The messages are those + of RAVEN's ``RAVEN:badGrRule`` error. """ - if not rule.strip(): + if rule is None or not rule.strip(): return None - gene_logger = logging.getLogger("cobra.core.gene") - was_disabled = gene_logger.disabled - gene_logger.disabled = True - try: - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - parsed = GPR.from_string(rule) - finally: - gene_logger.disabled = was_disabled - if parsed.body is None: - raise ValueError(f"Malformed grRule: {rule!r}") - return parsed + tokens = rule_tokens(rule) + + def factor(pos: int) -> tuple[ast.expr, int]: + if pos >= len(tokens): + raise ValueError(f"grRule ends after an operator: {rule}") + kind, text = tokens[pos] + if kind == "(": + node, pos = expr(pos + 1) + if pos >= len(tokens) or tokens[pos][0] != ")": + raise ValueError(f"Unbalanced parentheses in grRule: {rule}") + return node, pos + 1 + if kind == "gene": + return ast.Name(id=text, ctx=ast.Load()), pos + 1 + raise ValueError(f'Expected a gene but found "{text}" in grRule: {rule}') + + def chain(pos: int, sub, kind: str, op: ast.boolop) -> tuple[ast.expr, int]: + node, pos = sub(pos) + values = [node] + while pos < len(tokens) and tokens[pos][0] == kind: + node, pos = sub(pos + 1) + values.append(node) + if len(values) == 1: + return values[0], pos + return ast.BoolOp(op=op, values=values), pos + + def term(pos: int) -> tuple[ast.expr, int]: + return chain(pos, factor, "and", ast.And()) + + def expr(pos: int) -> tuple[ast.expr, int]: + return chain(pos, term, "or", ast.Or()) + + body, end = expr(0) + if end < len(tokens): + raise ValueError(f'Unexpected "{tokens[end][1]}" in grRule: {rule}') + return GPR(ast.Expression(body=body)) @dataclass(frozen=True) diff --git a/tests/test_manipulation_expand.py b/tests/test_manipulation_expand.py index a692ef3..21c1290 100644 --- a/tests/test_manipulation_expand.py +++ b/tests/test_manipulation_expand.py @@ -1,5 +1,6 @@ """Tests for expand_model (RAVEN expandModel.m) — splitting isozymes into reactions.""" import cobra +import pytest from raven_toolbox.manipulation import expand_model from raven_toolbox.manipulation.expand import gpr_to_dnf @@ -297,3 +298,22 @@ def test_first_copy_keeps_the_original_position_and_other_copies_are_appended(): "r1_EXP_1", "r2", "r3_EXP_1", "r4", "r1_EXP_2", "r3_EXP_2", "r3_EXP_3", ] assert model.reactions.index("r4") == 3 + + +@pytest.mark.parametrize( + "rule, expected", + [ + ("G1 && (G2 || G3)", [["G1", "G2"], ["G1", "G3"]]), + ("G1 Or G2", [["G1"], ["G2"]]), + ("G1&G2 or G3", [["G1&G2"], ["G3"]]), + ("", []), + ], +) +def test_dnf_from_string_uses_the_parsegrrule_grammar(rule, expected): + assert gpr_to_dnf(rule) == expected + + +@pytest.mark.parametrize("rule", ["G1 and", "(G1 or G2", "G1 G2", "()", "G1 or or G2"]) +def test_dnf_from_string_raises_on_a_malformed_rule(rule): + with pytest.raises(ValueError, match="grRule"): + gpr_to_dnf(rule) diff --git a/tests/test_utils_gpr.py b/tests/test_utils_gpr.py index 68f00f7..e29e3b7 100644 --- a/tests/test_utils_gpr.py +++ b/tests/test_utils_gpr.py @@ -1,8 +1,10 @@ """Tests for raven_toolbox.utils.gpr (GPR linting).""" +import re + import cobra import pytest -from raven_toolbox.utils import GPRIssue, find_non_dnf_grrules, is_dnf +from raven_toolbox.utils import GPRIssue, find_non_dnf_grrules, is_dnf, parse_gr_rule @pytest.mark.parametrize( @@ -90,10 +92,49 @@ def test_find_non_dnf_grrules_empty_when_all_clean(): ) def test_is_dnf_raises_on_malformed_rule(rule): """RAVEN's isDnfGrRule errors on these (RAVEN:badGrRule); cobra alone returns an empty GPR.""" - with pytest.raises(ValueError, match="Malformed grRule"): + with pytest.raises(ValueError, match="grRule"): is_dnf(rule) @pytest.mark.parametrize("rule", ["", " ", None]) def test_is_dnf_blank_rule_is_still_true(rule): assert is_dnf(rule) is True + + +@pytest.mark.parametrize( + "rule, expected", + [ + ("G1 && G2", "G1 and G2"), + ("G1 || G2", "G1 or G2"), + ("G1 Or G2 aNd G3", "G1 or (G2 and G3)"), + ("G1 AND (G2 | G3)", "G1 and (G2 or G3)"), + ("10 and 100", "10 and 100"), + ], +) +def test_parse_gr_rule_takes_the_parsegrrule_operators(rule, expected): + assert parse_gr_rule(rule).to_string() == expected + + +def test_parse_gr_rule_reads_an_operator_only_as_a_whole_word(): + """``G1&G2`` is one gene id in RAVEN's parseGrRule; cobra alone would read it as an AND.""" + assert sorted(parse_gr_rule("G1&G2 or RAND1").genes) == ["G1&G2", "RAND1"] + + +@pytest.mark.parametrize("rule", ["", " ", None]) +def test_parse_gr_rule_blank_is_none(rule): + assert parse_gr_rule(rule) is None + + +@pytest.mark.parametrize( + "rule, message", + [ + ("G1 and", "grRule ends after an operator"), + ("(G1 or G2", "Unbalanced parentheses in grRule"), + ("G1 G2", 'Unexpected "G2" in grRule'), + ("()", 'Expected a gene but found ")" in grRule'), + ("G1 xor G2", 'Unexpected "xor" in grRule'), + ], +) +def test_parse_gr_rule_raises_ravens_message(rule, message): + with pytest.raises(ValueError, match=re.escape(message)): + parse_gr_rule(rule) From 497bc1e8c91ca225c7794fa68a03d8366aaeaf5c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:52 +0100 Subject: [PATCH 09/35] analyse_topology: default seeds and targets are those of RAVEN's gapFillTopological --- CHANGELOG.md | 6 +++ src/raven_toolbox/gapfilling/topological.py | 46 +++++++++++++++------ tests/test_gapfilling_topological.py | 43 ++++++++++++++++++- 3 files changed, 81 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 27e9f3f..899a74b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `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; diff --git a/src/raven_toolbox/gapfilling/topological.py b/src/raven_toolbox/gapfilling/topological.py index af3d8a8..6d14601 100644 --- a/src/raven_toolbox/gapfilling/topological.py +++ b/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,18 +148,28 @@ 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) @@ -166,13 +177,16 @@ def analyse_topology( 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) @@ -180,7 +194,13 @@ def analyse_topology( # ---- 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 diff --git a/tests/test_gapfilling_topological.py b/tests/test_gapfilling_topological.py index f827385..1266656 100644 --- a/tests/test_gapfilling_topological.py +++ b/tests/test_gapfilling_topological.py @@ -15,7 +15,7 @@ def chain_with_gap(): C has no consumer and B is a dead-end — nothing past B is reachable. Template supplies B -(r3)-> C -(r4)-> D and EX_D (exchange). - Biomass: B -> (biomass reaction). + Objective: r2, so the default targets are B and C. With seeds = {A}, the scope reaches A and B (via r1) but NOT C (r2 is blocked because C has no further consumer). If targets = {D}, D is unreachable. @@ -34,6 +34,7 @@ def chain_with_gap(): r2.add_metabolites({B: -1, C: 1}) # B -> C, but C has no consumer draft.add_reactions([exa, r1, r2]) + draft.objective = "r2" # Template supplies C consumer and D pathway template = cobra.Model("template") @@ -116,3 +117,43 @@ def test_reversible_reaction_allows_reverse_scope(chain_with_gap): # r2 is irreversible, so scope still can't reach B (and hence A) from seed # C_c even with r1 now reversible; this only exercises the reverse-scope path. assert isinstance(result, TopologicalAnalysisResult) + + +def test_default_targets_are_the_objective_metabolites(chain_with_gap): + """As in RAVEN: substrates and products of the first objective reaction.""" + draft, template = chain_with_gap + result = analyse_topology(draft, template, verbose=False) + assert result.blocked_metabolites == set() + draft.objective = "EX_A" + draft.reactions.EX_A.lower_bound = 0 + with pytest.warns(UserWarning, match="no uptake exchange"): + result = analyse_topology(draft, template, verbose=False) + assert result.blocked_metabolites == {"A_c"} + + +def test_default_targets_need_an_objective(chain_with_gap): + draft, template = chain_with_gap + draft.objective = {} + with pytest.raises(ValueError, match="no objective function"): + analyse_topology(draft, template, verbose=False) + + +def test_default_seeds_follow_ravens_exchange_rule(): + """Seeds: metabolites of an only-substrate reaction with lb < 0 or an only-product one + with ub > 0, whatever the number of metabolites; nothing else.""" + model = cobra.Model("m") + mets = {mid: _met(mid) for mid in ("A_c", "B_c", "C_c", "D_c", "E_c")} + model.add_metabolites(list(mets.values())) + rxns = [ + ("EX_AB", {"A_c": -1, "B_c": -1}, -10, 0), # two-metabolite uptake: seeds A and B + ("EX_C", {"C_c": -1}, 0, 1000), # excretion only: not a seed + ("R_DE", {"D_c": -1, "E_c": 1}, -10, 10), # internal, reversible: not a seed + ] + for rid, stoich, lb, ub in rxns: + rxn = cobra.Reaction(rid, lower_bound=lb, upper_bound=ub) + rxn.add_metabolites({mets[m]: c for m, c in stoich.items()}) + model.add_reactions([rxn]) + model.objective = "R_DE" + result = analyse_topology(model, cobra.Model("t"), verbose=False) + assert result.reachable_metabolites == {"A_c", "B_c"} + assert result.blocked_metabolites == {"D_c", "E_c"} From c0a2e7430ef046fc3472bf81180db959c3445867 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:52 +0100 Subject: [PATCH 10/35] run_blast, run_diamond: threads default to every CPU the process may use, as in getBlast and getDiamond --- CHANGELOG.md | 3 +++ .../reconstruction/homology/blast.py | 20 ++++++++++++---- tests/test_reconstruction_blast.py | 24 ++++++++++++++++++- 3 files changed, 42 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 899a74b..01b6892 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,9 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `reconstruction.homology.run_blast` and `run_diamond` use every CPU the process may run on + when `threads` is not given, as RAVEN's `getBlast` and `getDiamond` 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 diff --git a/src/raven_toolbox/reconstruction/homology/blast.py b/src/raven_toolbox/reconstruction/homology/blast.py index 8ffe1e0..afbbb7e 100644 --- a/src/raven_toolbox/reconstruction/homology/blast.py +++ b/src/raven_toolbox/reconstruction/homology/blast.py @@ -44,6 +44,14 @@ def _as_list(x): return [x] if isinstance(x, (str, Path)) else list(x) +def _available_cpus() -> int: + """Number of CPUs this process may run on, the count MATLAB's ``feature('numcores')`` + reports as assigned by the OS and RAVEN passes to the aligners.""" + if hasattr(os, "sched_getaffinity"): + return max(1, len(os.sched_getaffinity(0))) + return max(1, os.cpu_count() or 1) + + def _run(cmd: list[str]) -> str: proc = subprocess.run(cmd, capture_output=True, text=True) if proc.returncode != 0: @@ -58,7 +66,7 @@ def run_blast( ref_fastas: Sequence[str | Path], *, evalue: float = 1e-4, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, blastp: str | Path | None = None, makeblastdb: str | Path | None = None, ) -> pd.DataFrame: @@ -68,7 +76,8 @@ def run_blast( reference proteome for template ``model_ids[i]``. Returns the hits DataFrame (filtered at - ``evalue``). Requires BLAST+ (`blastp`, `makeblastdb`). + ``evalue``). Requires BLAST+ (`blastp`, `makeblastdb`). ``threads`` + defaults to every CPU the process may use, as in RAVEN's ``getBlast``. The default matches MATLAB RAVEN's ``getBlast``, which hardcodes ``-evalue 10e-5`` -- that is 1e-4, not 1e-5. Unlike RAVEN's, this one is an @@ -81,6 +90,7 @@ def run_blast( ref_fastas = _as_list(ref_fastas) if len(model_ids) != len(ref_fastas): raise ValueError("model_ids and ref_fastas must have the same length.") + threads = _available_cpus() if threads is None else threads blastp = resolve_binary("blastp", binary=blastp) makeblastdb = resolve_binary("makeblastdb", binary=makeblastdb) outfmt = "10 " + " ".join(_OUTFMT_FIELDS) # 10 = CSV @@ -112,7 +122,7 @@ def run_diamond( ref_fastas: Sequence[str | Path], *, evalue: float = 1e-3, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, sensitivity: str = "--more-sensitive", diamond: str | Path | None = None, ) -> pd.DataFrame: @@ -123,7 +133,8 @@ def run_diamond( passed straight through as a DIAMOND CLI flag (e.g. ``"--more-sensitive"``, ``"--sensitive"``, or ``""`` to use DIAMOND's default fast mode). - Returns the hits DataFrame. Requires DIAMOND. + Returns the hits DataFrame. Requires DIAMOND. ``threads`` defaults to + every CPU the process may use, as in RAVEN's ``getDiamond``. The default matches MATLAB RAVEN's ``getDiamond``, which passes no ``--evalue`` at all and so inherits DIAMOND's own default of 1e-3. Stated @@ -136,6 +147,7 @@ def run_diamond( ref_fastas = _as_list(ref_fastas) if len(model_ids) != len(ref_fastas): raise ValueError("model_ids and ref_fastas must have the same length.") + threads = _available_cpus() if threads is None else threads diamond = resolve_binary("diamond", binary=diamond) frames = [] diff --git a/tests/test_reconstruction_blast.py b/tests/test_reconstruction_blast.py index 184150d..30b8c23 100644 --- a/tests/test_reconstruction_blast.py +++ b/tests/test_reconstruction_blast.py @@ -4,7 +4,12 @@ import pandas as pd import pytest -from raven_toolbox.reconstruction.homology import HIT_COLUMNS, blast_from_table, run_blast +from raven_toolbox.reconstruction.homology import ( + HIT_COLUMNS, + blast_from_table, + run_blast, + run_diamond, +) from raven_toolbox.reconstruction.homology.blast import _parse_tabular _SEQ = ( @@ -76,3 +81,20 @@ def test_run_blast_integration(tmp_path): # the reciprocal pair tgene<->ngene is found fwd = hits[(hits.from_gene == "tgene") & (hits.to_gene == "ngene")] assert not fwd.empty + + +@pytest.mark.parametrize("runner", [run_blast, run_diamond]) +def test_threads_default_to_every_available_cpu(runner, monkeypatch): + """As RAVEN's getBlast/getDiamond, which pass feature('numcores')'s assigned count.""" + import raven_toolbox.reconstruction.homology.blast as blast + + calls = [] + monkeypatch.setattr(blast, "_available_cpus", lambda: 7) + monkeypatch.setattr(blast, "resolve_binary", lambda name, binary=None: name) + monkeypatch.setattr(blast, "_run", lambda cmd: calls.append(cmd) or "") + runner("org", "org.faa", ["m"], ["m.faa"]) + searches = [c for c in calls if c[0] in ("blastp", "diamond") and "makedb" not in c] + assert searches + for cmd in searches: + flag = "-num_threads" if cmd[0] == "blastp" else "--threads" + assert cmd[cmd.index(flag) + 1] == "7" From 202f4e82af96b8387da59e0b019584f58d1cbed8 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:51:36 +0100 Subject: [PATCH 11/35] merge_compartments: keep the annotation of the kept copy; drop_single_metabolite_reactions defaults to False The merged metabolite kept the id, name, formula, charge and notes of the chosen copy and dropped its annotation. It now keeps the annotation too, as RAVEN's mergeCompartments keeps every field of the kept row. drop_single_metabolite_reactions defaults to False, the default of RAVEN's deleteRxnsWithOneMet. --- CHANGELOG.md | 4 ++- .../manipulation/compartments.py | 26 +++++++++++-------- tests/test_manipulation_compartments.py | 18 +++++++++++++ 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eebb605..feafe3d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,7 +12,9 @@ Milestones in the raven-toolbox port. For function-level status see * `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`). + 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`). * `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 diff --git a/src/raven_toolbox/manipulation/compartments.py b/src/raven_toolbox/manipulation/compartments.py index 6cd405e..9cfe333 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -9,6 +9,7 @@ """ from __future__ import annotations +import copy import warnings from collections.abc import Callable, Iterable @@ -24,7 +25,7 @@ def merge_compartments( merged_id: str = _MERGED_COMPARTMENT, merged_name: str = "system", base_metabolite: Callable[[cobra.Metabolite], str] | None = None, - drop_single_metabolite_reactions: bool = True, + drop_single_metabolite_reactions: bool = False, deduplicate_reactions: bool = True, ) -> tuple[cobra.Model, list[str], list[str]]: """Merge every metabolite of ``model`` into one ``merged_id`` compartment. @@ -48,19 +49,20 @@ def merge_compartments( yeast-GEM's ``s_####`` ids — need ``base_metabolite=lambda m: m.name`` (or an annotation key) so the compartments actually unify; pass the same callable to :func:`~raven_toolbox.localization.assign_compartments`. The key **only decides grouping**. The - merged metabolite keeps the id, name, formula, charge and notes of one of its copies, never an id - minted from the key: the copy already in ``merged_id`` if there is one, otherwise the copy whose id - sorts first (``ACCOA_c`` before ``ACCOA_m``); among several copies in ``merged_id`` the one whose id - sorts first. RAVEN's ``mergeCompartments`` applies the same rule. The id can therefore name a + merged metabolite keeps the id, name, formula, charge, annotation and notes of one of its copies, + never an id minted from the key: the copy already in ``merged_id`` if there is one, otherwise + the copy whose id sorts first (``ACCOA_c`` before ``ACCOA_m``); among several copies in + ``merged_id`` the one whose id sorts first. RAVEN's ``mergeCompartments`` applies the same rule. The id can therefore name a compartment that no longer exists. A name-based key like ``"1-acyl-sn-glycerol 3-phosphate (16:0)"`` unifies the copies while the merged id stays the model's own, already solver/SBML-safe identifier. Metabolites that share a key collapse into one entity in the merged compartment; their - stoichiometric contributions are summed per reaction. Reactions that *become* trivial - through the merge (e.g. ``A[c] → A[m]`` becomes ``A → A`` = nothing) are deleted by - default (RAVEN's ``deleteRxnsWithOneMet``). A reaction that already had a single - metabolite before merging is an exchange (``=> glucose``) and is always kept — - dropping those would silently remove the model's medium and its biomass reaction. + stoichiometric contributions are summed per reaction. A reaction whose metabolites all cancel + (``A[c] → A[m]`` becomes ``A → A``) is always deleted. With + ``drop_single_metabolite_reactions=True`` (RAVEN's ``deleteRxnsWithOneMet``; default ``False`` + on both sides) a reaction reduced to one metabolite by the merge is deleted too. A reaction that + already had a single metabolite before merging is an exchange (``=> glucose``) and is always + kept, as in RAVEN; dropping those would remove the model's medium and its biomass reaction. Reactions that become identical after merging are deduplicated (one survives). The objective is carried over to the returned model. """ @@ -84,7 +86,9 @@ def merge_compartments( mid = proto.id new_met = cobra.Metabolite(mid, name=proto.name, compartment=merged_id, formula=proto.formula, charge=proto.charge) - new_met.notes = dict(proto.notes or {}) + new_met.notes = copy.deepcopy(proto.notes or {}) + new_met.annotation = copy.deepcopy(proto.annotation or {}) + new_met._bound = proto._bound canonical[base] = new_met # Rewrite all reactions: replace each metabolite with its canonical, summing diff --git a/tests/test_manipulation_compartments.py b/tests/test_manipulation_compartments.py index 113d7bc..0c71c9f 100644 --- a/tests/test_manipulation_compartments.py +++ b/tests/test_manipulation_compartments.py @@ -111,6 +111,10 @@ def test_merge_compartments_keeps_single_met_reactions_when_asked(): # With keep, sym survives as a one-met reaction (consumes H). assert "sym" not in deleted_keep assert "sym" in {r.id for r in merged_keep.reactions} + merged_default, deleted_default, _ = merge_compartments(m) + # The default is RAVEN's deleteRxnsWithOneMet default, false. + assert "sym" in {r.id for r in merged_default.reactions} + assert deleted_default == deleted_keep def test_merge_compartments_keeps_pre_existing_exchange_reactions(): @@ -246,6 +250,20 @@ def test_merge_compartments_prefers_the_copy_in_the_merged_compartment(): assert {x.id for x in merged.metabolites} == {"ACCOA_z", "X_c"} +def test_merge_compartments_keeps_the_annotation_and_notes_of_the_chosen_copy(): + """RAVEN keeps every field of the kept row, metMiriams included.""" + m = _acetyl_coa_model(("ACCOA_m", "m"), ("ACCOA_c", "c")) + m.metabolites.ACCOA_m.annotation = {"chebi": "CHEBI:1"} + m.metabolites.ACCOA_c.annotation = {"chebi": "CHEBI:2", "kegg.compound": "C00024"} + m.metabolites.ACCOA_c.notes = {"note": "cytosolic copy"} + merged, _, _ = merge_compartments(m, base_metabolite=lambda mm: mm.name) + kept = merged.metabolites.ACCOA_c + assert dict(kept.annotation) == {"chebi": "CHEBI:2", "kegg.compound": "C00024"} + assert kept.notes == {"note": "cytosolic copy"} + kept.annotation["chebi"] = "changed" + assert m.metabolites.ACCOA_c.annotation["chebi"] == "CHEBI:2" + + # ----------------------------------------------------------------- copy_to_compartment def test_copy_to_compartment_basic(): From 2fb60d926b7f1f9fa2b6490e0fc17cb49aeb24ef Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:51:36 +0100 Subject: [PATCH 12/35] localization: the default metabolite key of a merged model matches the universal copy merge_compartments keeps the id of one copy, so a merged ACCOA_c sits in s and its own-compartment suffix rule left the key ACCOA_c. A universal ACCOA_m has the key ACCOA, so apply_assignment added a second copy of the metabolite in m. assign_compartments, apply_assignment, relocate_reactions and curation_priority now also strip a suffix that names a compartment of the model, the universal model, the scores or the placements. --- CHANGELOG.md | 4 ++ src/raven_toolbox/localization/assign.py | 48 ++++++++++++++++++++-- src/raven_toolbox/localization/certify.py | 7 +++- src/raven_toolbox/localization/curation.py | 14 ++++++- src/raven_toolbox/localization/relocate.py | 10 ++++- tests/test_localization_certify.py | 38 +++++++++++++++++ 6 files changed, 112 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index feafe3d..db3422b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,10 @@ Milestones in the raven-toolbox port. For function-level status see 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`). +* `localization.assign_compartments`, `apply_assignment`, `relocate_reactions` and + `curation_priority`: the default metabolite key also strips a trailing `_` suffix + that names a compartment of the model, the universal model, the scores or the placements, so + the merged `ACCOA_c` in `s` and a universal `ACCOA_m` share the key `ACCOA`. * `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 diff --git a/src/raven_toolbox/localization/assign.py b/src/raven_toolbox/localization/assign.py index 38cf74e..be8629d 100644 --- a/src/raven_toolbox/localization/assign.py +++ b/src/raven_toolbox/localization/assign.py @@ -20,7 +20,7 @@ """ from __future__ import annotations -from collections.abc import Callable, Mapping +from collections.abc import Callable, Iterable, Mapping from dataclasses import dataclass, field import cobra @@ -35,6 +35,40 @@ def _base_met(m: cobra.Metabolite) -> str: return m.id +def _default_base(compartments: Iterable[str]) -> Callable[[cobra.Metabolite], str]: + """The default compartment-agnostic key: the id without a trailing ``_`` suffix, + for the metabolite's own compartment or for any of ``compartments``. + + A metabolite of a :func:`~raven_toolbox.manipulation.merge_compartments` model keeps the id of + one of its copies, so ``ACCOA_c`` in ``s`` has the key ``ACCOA``, the key of ``ACCOA_m`` of a + universal model. Callers pass the compartments of the model, the universal model, the scores and + the placements. + """ + suffixes = sorted({f"_{c}" for c in compartments if c}, key=len, reverse=True) + + def base(m: cobra.Metabolite) -> str: + own = _base_met(m) + if own != m.id: + return own + for sfx in suffixes: + if len(m.id) > len(sfx) and m.id.endswith(sfx): + return m.id[: -len(sfx)] + return m.id + + return base + + +def _known_compartments(model: cobra.Model, universal: cobra.Model | None, + *others: Iterable[str]) -> set[str]: + """Compartment ids of ``model``, of ``universal`` and of ``others``.""" + comps = set(model.compartments) + if universal is not None: + comps |= set(universal.compartments) + for extra in others: + comps |= set(extra) + return comps + + def _rxn_compartment(rxn: cobra.Reaction) -> str | None: comps = {m.compartment for m in rxn.metabolites if m.compartment} return next(iter(comps)) if len(comps) == 1 else None @@ -122,10 +156,18 @@ def apply_assignment( ``default_compartment`` is the hub compartment star-topology transports route through. Pass the same ``base_metabolite`` and ``universal`` used in :func:`assign_compartments` so the same compartment-agnostic keying is used (existing per-compartment metabolites are reused rather than - duplicated) and the chosen gap-fill reactions are added. + duplicated) and the chosen gap-fill reactions are added. The default key is the id without a + trailing ``_`` suffix, for the metabolite's own compartment or for a compartment of + the model, ``universal`` or the proposal, so a merged ``ACCOA_c`` in ``s`` and a universal + ``ACCOA_m`` share the key ``ACCOA``. """ out = model.copy() - base = base_metabolite if base_metabolite is not None else _base_met + if base_metabolite is not None: + base = base_metabolite + else: + base = _default_base(_known_compartments( + model, universal, [default_compartment], *proposal.placements.values(), + [c for _, c in proposal.added_transports])) # (base key, compartment) -> existing metabolite, so a reaction moved into a compartment # reuses the species already there instead of creating a parallel duplicate. diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index 6e7f642..94d2ae8 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -41,7 +41,8 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, GrowthCondition, - _base_met, + _default_base, + _known_compartments, _rxn_compartment, apply_assignment, ) @@ -102,7 +103,7 @@ def _prepare_scope(model, scores, reactions_to_relocate, *, transportable, base_ pinned = [r for r in model.reactions if r.id not in movable_ids] pinned_comp = {r.id: _rxn_compartment(r) for r in pinned} - base = base_metabolite if base_metabolite is not None else _base_met + base = base_metabolite movable_base = {base(m) for r in movable for m in r.metabolites} transp = set(movable_base) if transportable is None else set(transportable) & movable_base @@ -649,6 +650,8 @@ def assign_compartments( if default_compartment not in compartments: raise ValueError(f"default_compartment={default_compartment!r} not in {compartments}") + if base_metabolite is None: + base_metabolite = _default_base(_known_compartments(model, universal, compartments)) sc = _prepare_scope(model, scores, reactions_to_relocate, transportable=transportable, base_metabolite=base_metabolite, biomass_reaction=biomass_reaction, min_growth=min_growth) diff --git a/src/raven_toolbox/localization/curation.py b/src/raven_toolbox/localization/curation.py index b03dbaa..2a690d8 100644 --- a/src/raven_toolbox/localization/curation.py +++ b/src/raven_toolbox/localization/curation.py @@ -51,7 +51,12 @@ import cobra -from raven_toolbox.localization.assign import AssignmentProposal, _base_met, _rxn_compartment +from raven_toolbox.localization.assign import ( + AssignmentProposal, + _default_base, + _known_compartments, + _rxn_compartment, +) from raven_toolbox.localization.scores import LocalizationScores __all__ = ["curation_priority"] @@ -170,7 +175,12 @@ def curation_priority( """ import pandas as pd - base = base_metabolite if base_metabolite is not None else _base_met + if base_metabolite is not None: + base = base_metabolite + else: + base = _default_base(_known_compartments( + model, universal, [default_compartment], scores.compartments, + *proposal.placements.values())) df = scores.df extra_currency = {_norm_key(x) for x in (currency_metabolites or ())} unplaced = set(proposal.unplaced_reactions) diff --git a/src/raven_toolbox/localization/relocate.py b/src/raven_toolbox/localization/relocate.py index 1cad005..950006c 100644 --- a/src/raven_toolbox/localization/relocate.py +++ b/src/raven_toolbox/localization/relocate.py @@ -27,7 +27,8 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, - _base_met, + _default_base, + _known_compartments, apply_assignment, ) from raven_toolbox.localization.curation import _is_impermeant @@ -116,7 +117,12 @@ def relocate_reactions( removes a transport of species ``b`` at a reaction's old compartment when nothing else uses ``b`` there any more. The transport repair only touches species of the moved reactions. """ - base = base_metabolite if base_metabolite is not None else _base_met + if base_metabolite is not None: + base = base_metabolite + else: + base = _default_base(_known_compartments( + model, universal, [default_compartment], decisions.values(), + *proposal.placements.values())) floor = float(min_growth if min_growth is not None else proposal.min_growth) warnings: list[str] = [] diff --git a/tests/test_localization_certify.py b/tests/test_localization_certify.py index c31883e..dc43bab 100644 --- a/tests/test_localization_certify.py +++ b/tests/test_localization_certify.py @@ -431,3 +431,41 @@ def test_growth_conditions_are_certified_per_medium(): assert res.certified assert res.growths["__primary__"] >= 1.0 assert res.growths["onB"] >= 1.0 + + +def _merged_two_compartment_draft(): + """A draft merged by merge_compartments: A has copies A_c and A_m, so the merged A keeps the id + A_c in compartment s.""" + from raven_toolbox.manipulation import merge_compartments + + m = cobra.Model("draft") + a_c = cobra.Metabolite("A_c", name="A", compartment="c") + a_m = cobra.Metabolite("A_m", name="A", compartment="m") + x_c = cobra.Metabolite("X_c", name="X", compartment="c") + r1 = cobra.Reaction("R1", lower_bound=0, upper_bound=1000) + r2 = cobra.Reaction("R2", lower_bound=0, upper_bound=1000) + ex = cobra.Reaction("EX_X", lower_bound=-10, upper_bound=1000) + m.add_reactions([r1, r2, ex]) + r1.add_metabolites({x_c: -1, a_c: 1}) + r2.add_metabolites({a_m: -1, x_c: 1}) + ex.add_metabolites({x_c: -1}) + merged, _, _ = merge_compartments(m) + assert merged.metabolites.A_c.compartment == "s" + return merged + + +def test_apply_assignment_matches_a_merged_metabolite_to_the_universal_copy(): + """The merged A_c (in s) and the universal A_m share the key A, so the gap-fill reaction lands + on the A that R1 makes in m instead of on a second copy.""" + merged = _merged_two_compartment_draft() + universal = cobra.Model("u") + u1 = cobra.Reaction("U1", lower_bound=0, upper_bound=1000) + universal.add_reactions([u1]) + u1.add_metabolites({cobra.Metabolite("A_m", name="A", compartment="m"): -1, + cobra.Metabolite("Y_m", name="Y", compartment="m"): 1}) + proposal = AssignmentProposal(placements={"R1": ["m"]}, added_reactions=["U1"]) + out = apply_assignment(merged, proposal, universal=universal) + a_in_m = [x for x in out.metabolites if x.name == "A" and x.compartment == "m"] + assert len(a_in_m) == 1 + assert a_in_m[0] in out.reactions.R1.metabolites + assert a_in_m[0] in out.reactions.U1.metabolites From 422304219f93895daa6166b4a09e3c7277b51f30 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:51:49 +0100 Subject: [PATCH 13/35] parse_gr_rule: whitespace is what MATLAB's isspace counts The tokeniser split words at Python's Unicode whitespace, so a no-break space (U+00A0), U+0085, U+2007 and U+202F separated words, while RAVEN's parseGrRule, which uses isspace, keeps them in the word and refuses G1and G2. The whitespace class is now the list of code points for which isspace(char(0:65535)) is true in MATLAB R2024b. --- CHANGELOG.md | 3 ++- src/raven_toolbox/utils/gpr.py | 16 +++++++++++----- tests/test_utils_gpr.py | 29 +++++++++++++++++++++++++++++ 3 files changed, 42 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 01b6892..c10b656 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,8 @@ Milestones in the raven-toolbox port. For function-level status see 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. + 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. * `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 diff --git a/src/raven_toolbox/utils/gpr.py b/src/raven_toolbox/utils/gpr.py index c84f763..5b88671 100644 --- a/src/raven_toolbox/utils/gpr.py +++ b/src/raven_toolbox/utils/gpr.py @@ -108,7 +108,11 @@ def is_dnf(gpr: GPR | str | None) -> bool: return _is_dnf_node(parsed.body) -_RULE_TOKEN = re.compile(r"\(|\)|[^\s()]+") +# The characters for which MATLAB's isspace is true, the whitespace of RAVEN's parseGrRule. Python's +# str.isspace also counts U+0085, U+00A0, U+2007 and U+202F; here they are part of a word. +_RULE_SPACE = r"\t\n\x0b\x0c\r\x1c-\x20\u1680\u2000-\u2006\u2008-\u200a\u2028\u2029\u205f\u3000" +_RULE_TOKEN = re.compile(rf"\(|\)|[^{_RULE_SPACE}()]+") +_RULE_BLANK = re.compile(rf"[{_RULE_SPACE}]*") _AND_WORDS = frozenset({"and", "&", "&&"}) _OR_WORDS = frozenset({"or", "|", "||"}) @@ -116,9 +120,11 @@ def is_dnf(gpr: GPR | str | None) -> bool: def rule_tokens(rule: str) -> list[tuple[str, str]]: """Split a gene rule into ``(kind, text)`` pairs; kind is ``(``, ``)``, ``and``, ``or`` or ``gene``. - A word is everything up to whitespace or a bracket. The operators are ``and``/``or`` in any - case and the symbols ``&``, ``&&``, ``|`` and ``||``, each as a whole word, as RAVEN's - ``parseGrRule`` reads them; every other word is a gene id, so ``G1&G2`` is one gene. + A word is everything up to whitespace or a bracket, whitespace being the characters for which + MATLAB's ``isspace`` is true (not U+00A0, U+0085, U+2007 or U+202F). The operators are + ``and``/``or`` in any case and the symbols ``&``, ``&&``, ``|`` and ``||``, each as a whole + word, as RAVEN's ``parseGrRule`` reads them; every other word is a gene id, so ``G1&G2`` is + one gene. """ tokens = [] for match in _RULE_TOKEN.finditer(rule): @@ -167,7 +173,7 @@ def parse_gr_rule(rule: str | None) -> GPR | None: parentheses, or two genes without an operator between them. The messages are those of RAVEN's ``RAVEN:badGrRule`` error. """ - if rule is None or not rule.strip(): + if rule is None or _RULE_BLANK.fullmatch(rule): return None tokens = rule_tokens(rule) diff --git a/tests/test_utils_gpr.py b/tests/test_utils_gpr.py index e29e3b7..0c3d2ab 100644 --- a/tests/test_utils_gpr.py +++ b/tests/test_utils_gpr.py @@ -138,3 +138,32 @@ def test_parse_gr_rule_blank_is_none(rule): def test_parse_gr_rule_raises_ravens_message(rule, message): with pytest.raises(ValueError, match=re.escape(message)): parse_gr_rule(rule) + + +# The code points for which MATLAB's isspace is true, enumerated over char(0:65535). +_MATLAB_ISSPACE = [9, 10, 11, 12, 13, 28, 29, 30, 31, 32, 5760, *range(8192, 8199), 8200, 8201, + 8202, 8232, 8233, 8287, 12288] + + +def test_parse_gr_rule_whitespace_is_matlab_isspace(): + for code in range(65536): + rule = f"G1{chr(code)}and G2" + if code in _MATLAB_ISSPACE: + assert sorted(parse_gr_rule(rule).genes) == ["G1", "G2"], code + elif chr(code) not in "()": + with pytest.raises(ValueError): + parse_gr_rule(rule) + + +@pytest.mark.parametrize("char", [" ", "\u0085", " ", " "]) +def test_parse_gr_rule_a_space_that_isspace_refuses_is_part_of_a_word(char): + with pytest.raises(ValueError, match=re.escape(f'Unexpected "G2" in grRule: G1{char}and G2')): + parse_gr_rule(f"G1{char}and G2") + assert sorted(parse_gr_rule(f"G1{char}G2").genes) == [f"G1{char}G2"] + assert parse_gr_rule(char) is not None + + +@pytest.mark.parametrize("char", ["\x1c", " ", "　"]) +def test_parse_gr_rule_isspace_characters_separate_words(char): + assert parse_gr_rule(f"G1{char}or{char}G2").to_string() == "G1 or G2" + assert parse_gr_rule(char * 2) is None From 35417c0dfe23aa34341875f1bfff2e830bd24fab Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:51:49 +0100 Subject: [PATCH 14/35] batch_curate: store the gene rule as parse_gr_rule reads it The rule was checked with parse_gr_rule, for which G1&G2 is one gene id, and then stored as a string that cobra reads as G1 and G2. The reaction now gets the parsed rule, so its genes are the ones that were checked, as in RAVEN's curateModelFromTables. --- CHANGELOG.md | 3 +++ src/raven_toolbox/curation/batch.py | 34 +++++++++-------------------- tests/test_curation.py | 20 +++++++++++++++++ 3 files changed, 33 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c10b656..5a8b509 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,9 @@ Milestones in the raven-toolbox port. For function-level status see 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`). * `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 diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index d42565b..0545e61 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -12,8 +12,9 @@ import cobra import numpy as np import pandas as pd +from cobra.core.gene import GPR -from raven_toolbox.utils.gpr import _RULE_TOKEN, parse_gr_rule, rule_tokens +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 @@ -435,10 +436,6 @@ def _row_has_content(row, columns=None) -> bool: return any(_has_value(v) for v in values) -#: Spellings of the operators that cobra reads as they are. -_COBRA_OPERATORS = frozenset({"and", "AND", "&", "or", "OR", "|"}) - - def _rule_genes(rule: str) -> list[str]: """Gene ids of a gene rule, in order; :class:`ValueError` when the rule is not well formed. @@ -462,23 +459,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() - kind = rule_tokens(word)[0][0] - if word in _COBRA_OPERATORS or kind not in ("and", "or"): - continue - replacement = kind - 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: @@ -848,7 +834,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. @@ -904,7 +890,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")): diff --git a/tests/test_curation.py b/tests/test_curation.py index d4d5c5b..eca5aaa 100644 --- a/tests/test_curation.py +++ b/tests/test_curation.py @@ -1285,6 +1285,26 @@ def test_a_valid_gene_rule_in_any_spelling_is_stored(rule, genes): assert sorted(g.id for g in _by_name(m, "fresh").genes) == genes +@pytest.mark.parametrize("name, ctable", [("fresh", ("-2", "2")), ("ATP hydrolysis", ("-1", "1"))]) +@pytest.mark.parametrize("rule, genes", [ + ("G1&G2 or G3", ["G1&G2", "G3"]), + ("G1|G2", ["G1|G2"]), +]) +def test_an_operator_glued_to_a_gene_id_stays_part_of_the_stored_gene(name, ctable, rule, genes): + """RAVEN's parseGrRule reads G1&G2 as one gene id; the stored rule has the same genes.""" + m = _base_model() + rxns_df = pd.DataFrame([{"rxnNames": name, "grRules": rule}]) + coeffs_df = pd.DataFrame([ + {"rxnNames": name, "metNames": "ATP", "comps": "c", "coefficient": ctable[0]}, + {"rxnNames": name, "metNames": "ADP", "comps": "c", "coefficient": ctable[1]}]) + batch_curate(m, genes_df=_genes_df(*genes), rxns_df=rxns_df, + rxns_coeffs_df=coeffs_df, rxn_id_prefix="r_") + rxn = _by_name(m, name) + assert sorted(g.id for g in rxn.genes) == genes + assert rxn.gene_reaction_rule == rule + assert not {"G1", "G2"} & {g.id for g in m.genes} + + # --- header-only tables, atomicity of several rows, genes of new reactions, headers --------- _RXN_HEADER = "rxnIdx\trxnNames\tlb\tub\n" From 2201873efbb2f48af9ab536b9e8e62a801857b95 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:51:49 +0100 Subject: [PATCH 15/35] run_hmmsearch, get_kegg_model_from_sequences: threads default to every CPU the process may use The default was the machine's CPU count minus one. It is now the number of CPUs the process may run on, as in run_blast and run_diamond and in RAVEN's getKEGGModelForOrganism, which passes feature('numcores')'s assigned count to hmmsearch. --- CHANGELOG.md | 7 +++--- .../reconstruction/kegg/query.py | 13 ++++++----- tests/test_reconstruction_kegg_query.py | 22 +++++++++++++++++++ 3 files changed, 34 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a8b509..76d43f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,10 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased -* `reconstruction.homology.run_blast` and `run_diamond` use every CPU the process may run on - when `threads` is not given, as RAVEN's `getBlast` and `getDiamond` do (was: the machine's - CPU count minus one). +* `reconstruction.homology.run_blast` and `run_diamond`, and `reconstruction.kegg.run_hmmsearch` + and `get_kegg_model_from_sequences`, use every CPU the process may run on when `threads` is + not given, as RAVEN's `getBlast`, `getDiamond` 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 diff --git a/src/raven_toolbox/reconstruction/kegg/query.py b/src/raven_toolbox/reconstruction/kegg/query.py index f949325..b49e017 100644 --- a/src/raven_toolbox/reconstruction/kegg/query.py +++ b/src/raven_toolbox/reconstruction/kegg/query.py @@ -19,7 +19,6 @@ from __future__ import annotations import math -import os import subprocess import tempfile from pathlib import Path @@ -29,6 +28,7 @@ from raven_toolbox.binaries import resolve_binary from raven_toolbox.io.yaml import read_yaml_model +from raven_toolbox.reconstruction.homology.blast import _available_cpus from raven_toolbox.reconstruction.kegg.assemble import assemble_model_from_ko_genes from raven_toolbox.reconstruction.kegg.parse import _resolve_artefact, read_kegg_table @@ -50,7 +50,7 @@ def run_hmmsearch( fasta: str | Path, library: str | Path, *, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, hmmsearch: str | Path | None = None, ) -> str: """Search the profile ``library`` against proteome ``fasta``; return tblout text. @@ -59,8 +59,10 @@ def run_hmmsearch( against ``fasta`` (target) — the fast search direction, parallelised with ``--cpu``, and no ``hmmpress`` needed. ``-Z`` is fixed to the profile count so the per-hit E-values match the convention :func:`assign_kos` is calibrated against - (identical to a ``hmmscan`` against the same library). + (identical to a ``hmmscan`` against the same library). ``threads`` defaults to every + CPU the process may use, as in RAVEN's ``getKEGGModelForOrganism``. """ + threads = _available_cpus() if threads is None else threads exe = resolve_binary("hmmsearch", binary=hmmsearch) nprofiles = _count_profiles(library) z_opt = ["-Z", str(nprofiles)] if nprofiles > 0 else [] @@ -190,14 +192,15 @@ def get_kegg_model_from_sequences( keep_undefined_stoich: bool = True, keep_incomplete: bool = True, keep_general: bool = False, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, hmmsearch: str | Path | None = None, ) -> cobra.Model: """Reconstruct a draft model for a proteome by HMM-searching the KO library. Searches the ``library`` against ``fasta`` (3b.3), assigns KOs (:func:`assign_kos`), and assembles the model against ``reference_model`` / - ``ko_reaction``. Genes are the query proteome's identifiers. + ``ko_reaction``. Genes are the query proteome's identifiers. ``threads`` is passed to + :func:`run_hmmsearch`. """ if model_id is None: # RAVEN always sets model.id (to organismID); default it here so the diff --git a/tests/test_reconstruction_kegg_query.py b/tests/test_reconstruction_kegg_query.py index bce7003..23a7d85 100644 --- a/tests/test_reconstruction_kegg_query.py +++ b/tests/test_reconstruction_kegg_query.py @@ -174,3 +174,25 @@ def test_prune_orthology_keeps_only_matched_kos(): # Default (organism path) keeps the full reference KO list. full, _ = assemble_model_from_ko_genes(ref, ko_reaction, ko_to_genes) assert full.reactions.get_by_id("R1").annotation["kegg.orthology"] == ["K1", "K2", "K3"] + + +def test_hmmsearch_threads_default_to_every_available_cpu(tmp_path, monkeypatch): + """As RAVEN's getKEGGModelForOrganism, which passes feature('numcores')'s assigned count.""" + import subprocess + + import raven_toolbox.reconstruction.kegg.query as query + + calls = [] + monkeypatch.setattr(query, "_available_cpus", lambda: 7) + monkeypatch.setattr(query, "resolve_binary", lambda name, binary=None: name) + + def fake_run(cmd, **kwargs): + calls.append(cmd) + open(cmd[cmd.index("--tblout") + 1], "w").close() + return subprocess.CompletedProcess(cmd, 0, "", "") + + monkeypatch.setattr(query.subprocess, "run", fake_run) + library = tmp_path / "lib.hmm" + library.write_text("HMMER3/f\n") + query.run_hmmsearch("org.faa", library) + assert calls[0][calls[0].index("--cpu") + 1] == "7" From 134020ba08c1e8615aeec47f128a0dca50ef8b14 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:26:18 +0100 Subject: [PATCH 16/35] remove_genes: constrain blocked reactions by default and remove unused genes, as removeGenes --- CHANGELOG.md | 6 ++++++ src/raven_toolbox/manipulation/remove.py | 27 ++++++++++++++++++------ tests/test_manipulation_remove.py | 23 ++++++++++++++++++++ 3 files changed, 50 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index db3422b..91fb58d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ 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. 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. diff --git a/src/raven_toolbox/manipulation/remove.py b/src/raven_toolbox/manipulation/remove.py index 43f966f..e1df9f2 100644 --- a/src/raven_toolbox/manipulation/remove.py +++ b/src/raven_toolbox/manipulation/remove.py @@ -11,8 +11,9 @@ * ``remove_genes`` — cobra's ``cobra.manipulation.remove_genes`` already rewrites GPRs through the boolean AST (removing one gene of ``A and B`` empties the rule, of ``A or B`` keeps the other) — exactly RAVEN's intent, without its - ``eval``. The gap is RAVEN's default of **constraining** flux-blocked reactions - to zero instead of deleting them; exposed as ``blocked_reactions``. + ``eval``. What it adds is RAVEN's default of **constraining** flux-blocked + reactions to zero instead of deleting them (``blocked_reactions``) and RAVEN's + removal of genes that no remaining reaction uses. """ from __future__ import annotations @@ -66,7 +67,7 @@ def remove_genes( model: cobra.Model, genes: str | Gene | Iterable, *, - blocked_reactions: str = "remove", + blocked_reactions: str = "constrain", remove_orphans: bool = False, ) -> list[str]: """Remove genes and handle reactions left unable to carry flux. @@ -74,12 +75,20 @@ def remove_genes( GPR rewriting (with correct AND/OR semantics) and gene deletion are done by cobra; this adds a policy for reactions whose GPR becomes empty (no enzyme left): - * ``"remove"`` — delete them (cobra's default). - * ``"constrain"`` — keep them but set bounds to ``(0, 0)``. + * ``"constrain"`` (default, as RAVEN's ``removeGenes``) — keep them with an empty + GPR and bounds ``(0, 0)``. + * ``"remove"`` — delete them. * ``"keep"`` — leave them with an empty GPR and unchanged bounds. + Besides the requested genes, every gene that no remaining reaction named in its + GPR before the removal is deleted too, as RAVEN's ``removeGenes`` does: genes + that no reaction uses (orphan genes, also when not requested) and, with + ``"remove"``, genes used only by the deleted reactions. A gene whose only + reaction is blocked and kept (``"constrain"`` or ``"keep"``), such as the other + subunit of a complex that lost one, stays in ``model.genes``. + ``remove_orphans`` (only meaningful with ``blocked_reactions="remove"``) - passes through to cobra: drop metabolites *and* genes orphaned by the removal. + passes through to cobra: also drop metabolites orphaned by the removal. Returns ------- @@ -97,6 +106,7 @@ def remove_genes( present = [gid for gid in requested if gid in model.genes] if not present: return [] + named_before = {r.id: {g.id for g in r.genes} for r in model.reactions} # Reactions touched by these genes that currently have a GPR. affected: set[str] = set() @@ -117,4 +127,9 @@ def remove_genes( for rid in blocked: model.reactions.get_by_id(rid).bounds = (0, 0) + used = set().union(*(named_before[r.id] for r in model.reactions)) + unused = [g.id for g in model.genes if g.id not in used] + if unused: + _cobra_remove_genes(model, unused, remove_reactions=False) + return sorted(blocked) diff --git a/tests/test_manipulation_remove.py b/tests/test_manipulation_remove.py index 3bfe0be..3162c5c 100644 --- a/tests/test_manipulation_remove.py +++ b/tests/test_manipulation_remove.py @@ -95,3 +95,26 @@ def test_remove_genes_absent_gene_is_noop(model): def test_remove_genes_bad_policy(model): with pytest.raises(ValueError, match="blocked_reactions"): remove_genes(model, ["G1"], blocked_reactions="explode") + + +def test_remove_genes_default_constrains_as_raven(model): + n_reactions = len(model.reactions) + blocked = remove_genes(model, ["G1"]) + assert blocked == ["R1"] + assert len(model.reactions) == n_reactions + r1 = model.reactions.get_by_id("R1") + assert r1.bounds == (0, 0) + assert r1.gene_reaction_rule == "" + + +def test_remove_genes_drops_unused_genes(model): + model.genes.append(cobra.Gene("G9")) # orphan gene, not requested + remove_genes(model, ["G1"]) + # G2 was named by R1, which is kept (constrained), so it stays; G9 goes. + assert sorted(g.id for g in model.genes) == ["G2", "G3", "G4"] + + +def test_remove_genes_remove_mode_drops_genes_of_deleted_reactions(model): + model.genes.append(cobra.Gene("G9")) + remove_genes(model, ["G1"], blocked_reactions="remove") + assert sorted(g.id for g in model.genes) == ["G3", "G4"] From 9aadd82e46a69082c65d8805432559d654989d7c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:26:35 +0100 Subject: [PATCH 17/35] tests: merge_linear marks the reaction flipped in a later pass as reversed --- tests/test_init_merge.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_init_merge.py b/tests/test_init_merge.py index 29254f0..0604c65 100644 --- a/tests/test_init_merge.py +++ b/tests/test_init_merge.py @@ -119,6 +119,24 @@ def test_groups_are_numbered_by_their_first_reaction(): assert sorted(r.id for r in reduced.reactions) == ["R1", "R5"] +def test_reaction_flipped_in_later_pass_is_marked_reversed(): + """R1 and R2 merge through A in the first pass; the second pass merges R1 and R4 through + Y and flips R4, so R4 (not R3, the third reaction) is the reversed one.""" + m = cobra.Model("lateflip") + mets = {i: cobra.Metabolite(i, name=i, compartment="c") for i in "AYQW"} + spec = [("R1", {"A": 1, "Y": 1}, 0), ("R2", {"A": -1, "Y": 1}, 0), ("R3", {"Q": -1}, 0), + ("R4", {"W": -1, "Y": 1}, -1000), ("R5", {"W": 1}, 0), ("R6", {"W": -1}, 0)] + for rid, stoichiometry, lb in spec: + r = cobra.Reaction(rid, lower_bound=lb, upper_bound=1000) + r.add_metabolites({mets[k]: v for k, v in stoichiometry.items()}) + m.add_reactions([r]) + + reduced, orig_ids, group_ids, reversed_rxns = merge_linear(m) + assert group_ids == [1, 1, 0, 1, 0, 0] + assert [rid for rid, flag in zip(orig_ids, reversed_rxns, strict=True) if flag] == ["R4"] + assert sorted(r.id for r in reduced.reactions) == ["R1", "R3", "R5", "R6"] + + def test_group_scores_zero_handling(): """Genuine-zero score → 0.01; a group cancelling to zero with nonzero members → 0.01.""" reduced, orig_ids, group_ids, _ = merge_linear(make_test_model()) From 635de5cbb39d89ae6e5c5a531136c7a7a6c6f4c1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:31:35 +0100 Subject: [PATCH 18/35] read_yaml_model: notes that write_yaml_model writes as a notes mapping go back into notes write_yaml_model writes the notes without a RAVEN field as a notes mapping, which read_yaml_model nested under notes["note"]. They are now spread back into notes, with the free-text note under "note". --- src/raven_toolbox/io/yaml.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 325d69a..3b2dddc 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -44,6 +44,7 @@ import re import warnings from collections import OrderedDict +from collections.abc import Mapping from datetime import date from pathlib import Path @@ -141,6 +142,9 @@ def _capture_entry_fields(entries, fields): Returns a list of ``{notes_key: value}`` dicts aligned with ``entries`` (so cobra never sees these keys), to be attached to the built objects afterwards. + A ``notes`` mapping, which :func:`write_yaml_model` writes for notes that have no + RAVEN field, is spread back into the notes dict, with its ``note`` entry as the + free-text note. """ captured = [] for entry in entries: @@ -148,6 +152,8 @@ def _capture_entry_fields(entries, fields): for yaml_key, notes_key in fields: if yaml_key in entry: notes[notes_key] = entry.pop(yaml_key) + if isinstance(notes.get("note"), Mapping): + notes.update(notes.pop("note")) captured.append(notes) return captured From 652e8d347b68000fdbca5de4be07a36c8a85a37c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:31:35 +0100 Subject: [PATCH 19/35] localization: the default metabolite key is the base id merge_compartments stores The default key guessed a compartment suffix from the compartments each function saw, so assign_compartments and apply_assignment could key one metabolite differently (I and I_g), and distinct species such as X_c and X_m in s, or vitamin_a and vitamin_c, got one key. merge_compartments now stores the base id of each merged metabolite in notes["merged_base_id"] (the own id when two would coincide), and the default key is that stored id, otherwise the id without the suffix of the metabolite's own compartment. assign_compartments, apply_assignment, relocate_reactions and curation_priority raise ValueError when two metabolites in one compartment have the same key. --- CHANGELOG.md | 14 ++- src/raven_toolbox/localization/assign.py | 86 +++++++-------- src/raven_toolbox/localization/certify.py | 9 +- src/raven_toolbox/localization/curation.py | 12 +-- src/raven_toolbox/localization/relocate.py | 12 +-- .../manipulation/compartments.py | 33 +++++- tests/test_localization_certify.py | 101 ++++++++++++++++++ tests/test_manipulation_compartments.py | 43 +++++++- 8 files changed, 235 insertions(+), 75 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index db3422b..8186386 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,10 +15,16 @@ Milestones in the raven-toolbox port. For function-level status see 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`). -* `localization.assign_compartments`, `apply_assignment`, `relocate_reactions` and - `curation_priority`: the default metabolite key also strips a trailing `_` suffix - that names a compartment of the model, the universal model, the scores or the placements, so - the merged `ACCOA_c` in `s` and a universal `ACCOA_m` share the key `ACCOA`. +* `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). +* `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"]`). * `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 diff --git a/src/raven_toolbox/localization/assign.py b/src/raven_toolbox/localization/assign.py index be8629d..acf879f 100644 --- a/src/raven_toolbox/localization/assign.py +++ b/src/raven_toolbox/localization/assign.py @@ -25,48 +25,43 @@ import cobra +from raven_toolbox.manipulation.compartments import MERGED_BASE_ID_NOTE + __all__ = ["AssignmentProposal", "GrowthCondition", "apply_assignment"] # --------------------------------------------------------------------------- helpers def _base_met(m: cobra.Metabolite) -> str: + """The default compartment-agnostic key of a metabolite. + + The base id that :func:`~raven_toolbox.manipulation.merge_compartments` stored in + ``notes["merged_base_id"]``, if there is one; otherwise the id without the trailing + ``_`` suffix of the metabolite's own compartment. + """ + stored = (m.notes or {}).get(MERGED_BASE_ID_NOTE) + if isinstance(stored, str) and stored: + return stored if m.compartment and m.id.endswith(f"_{m.compartment}"): return m.id[: -(len(m.compartment) + 1)] return m.id -def _default_base(compartments: Iterable[str]) -> Callable[[cobra.Metabolite], str]: - """The default compartment-agnostic key: the id without a trailing ``_`` suffix, - for the metabolite's own compartment or for any of ``compartments``. - - A metabolite of a :func:`~raven_toolbox.manipulation.merge_compartments` model keeps the id of - one of its copies, so ``ACCOA_c`` in ``s`` has the key ``ACCOA``, the key of ``ACCOA_m`` of a - universal model. Callers pass the compartments of the model, the universal model, the scores and - the placements. - """ - suffixes = sorted({f"_{c}" for c in compartments if c}, key=len, reverse=True) - - def base(m: cobra.Metabolite) -> str: - own = _base_met(m) - if own != m.id: - return own - for sfx in suffixes: - if len(m.id) > len(sfx) and m.id.endswith(sfx): - return m.id[: -len(sfx)] - return m.id - - return base - - -def _known_compartments(model: cobra.Model, universal: cobra.Model | None, - *others: Iterable[str]) -> set[str]: - """Compartment ids of ``model``, of ``universal`` and of ``others``.""" - comps = set(model.compartments) - if universal is not None: - comps |= set(universal.compartments) - for extra in others: - comps |= set(extra) - return comps +def _check_unique_keys(metabolites: Iterable[cobra.Metabolite], + base: Callable[[cobra.Metabolite], str]) -> None: + """Raise :class:`ValueError` when two metabolites in one compartment have the same key, which + would make them one species.""" + seen: dict[tuple[str, str], str] = {} + clashes: list[str] = [] + for m in metabolites: + key = (base(m), m.compartment) + if key in seen and seen[key] != m.id: + clashes.append(f"{seen[key]} and {m.id} (key {key[0]!r} in {key[1]!r})") + else: + seen.setdefault(key, m.id) + if clashes: + raise ValueError("metabolites with the same key in one compartment: " + + "; ".join(clashes[:5]) + + "; pass a base_metabolite that tells them apart") def _rxn_compartment(rxn: cobra.Reaction) -> str | None: @@ -156,18 +151,14 @@ def apply_assignment( ``default_compartment`` is the hub compartment star-topology transports route through. Pass the same ``base_metabolite`` and ``universal`` used in :func:`assign_compartments` so the same compartment-agnostic keying is used (existing per-compartment metabolites are reused rather than - duplicated) and the chosen gap-fill reactions are added. The default key is the id without a - trailing ``_`` suffix, for the metabolite's own compartment or for a compartment of - the model, ``universal`` or the proposal, so a merged ``ACCOA_c`` in ``s`` and a universal - ``ACCOA_m`` share the key ``ACCOA``. + duplicated) and the chosen gap-fill reactions are added. The default key is the base id that + :func:`~raven_toolbox.manipulation.merge_compartments` stored on a merged metabolite, otherwise + the id without the suffix of the metabolite's own compartment. :class:`ValueError` when two + metabolites in one compartment of ``model`` have the same key. """ out = model.copy() - if base_metabolite is not None: - base = base_metabolite - else: - base = _default_base(_known_compartments( - model, universal, [default_compartment], *proposal.placements.values(), - [c for _, c in proposal.added_transports])) + base = base_metabolite if base_metabolite is not None else _base_met + _check_unique_keys(out.metabolites, base) # (base key, compartment) -> existing metabolite, so a reaction moved into a compartment # reuses the species already there instead of creating a parallel duplicate. @@ -180,10 +171,13 @@ def resolve(template: cobra.Metabolite, compartment: str) -> cobra.Metabolite: met = index.get(key) if met is None: new_id = f"{template.id}__{compartment}" - met = (out.metabolites.get_by_id(new_id) if new_id in out.metabolites - else cobra.Metabolite(new_id, name=template.name, compartment=compartment, - formula=template.formula, charge=template.charge)) - if new_id not in out.metabolites: + if new_id in out.metabolites: + met = out.metabolites.get_by_id(new_id) + else: + met = cobra.Metabolite(new_id, name=template.name, compartment=compartment, + formula=template.formula, charge=template.charge) + if MERGED_BASE_ID_NOTE in (template.notes or {}): + met.notes[MERGED_BASE_ID_NOTE] = template.notes[MERGED_BASE_ID_NOTE] out.add_metabolites([met]) index[key] = met return met diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index 94d2ae8..fd43cd1 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -41,8 +41,8 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, GrowthCondition, - _default_base, - _known_compartments, + _base_met, + _check_unique_keys, _rxn_compartment, apply_assignment, ) @@ -103,7 +103,7 @@ def _prepare_scope(model, scores, reactions_to_relocate, *, transportable, base_ pinned = [r for r in model.reactions if r.id not in movable_ids] pinned_comp = {r.id: _rxn_compartment(r) for r in pinned} - base = base_metabolite + base = base_metabolite if base_metabolite is not None else _base_met movable_base = {base(m) for r in movable for m in r.metabolites} transp = set(movable_base) if transportable is None else set(transportable) & movable_base @@ -650,8 +650,7 @@ def assign_compartments( if default_compartment not in compartments: raise ValueError(f"default_compartment={default_compartment!r} not in {compartments}") - if base_metabolite is None: - base_metabolite = _default_base(_known_compartments(model, universal, compartments)) + _check_unique_keys(model.metabolites, base_metabolite or _base_met) sc = _prepare_scope(model, scores, reactions_to_relocate, transportable=transportable, base_metabolite=base_metabolite, biomass_reaction=biomass_reaction, min_growth=min_growth) diff --git a/src/raven_toolbox/localization/curation.py b/src/raven_toolbox/localization/curation.py index 2a690d8..0ec0e46 100644 --- a/src/raven_toolbox/localization/curation.py +++ b/src/raven_toolbox/localization/curation.py @@ -53,8 +53,8 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, - _default_base, - _known_compartments, + _base_met, + _check_unique_keys, _rxn_compartment, ) from raven_toolbox.localization.scores import LocalizationScores @@ -175,12 +175,8 @@ def curation_priority( """ import pandas as pd - if base_metabolite is not None: - base = base_metabolite - else: - base = _default_base(_known_compartments( - model, universal, [default_compartment], scores.compartments, - *proposal.placements.values())) + base = base_metabolite if base_metabolite is not None else _base_met + _check_unique_keys(model.metabolites, base) df = scores.df extra_currency = {_norm_key(x) for x in (currency_metabolites or ())} unplaced = set(proposal.unplaced_reactions) diff --git a/src/raven_toolbox/localization/relocate.py b/src/raven_toolbox/localization/relocate.py index 950006c..3e254e1 100644 --- a/src/raven_toolbox/localization/relocate.py +++ b/src/raven_toolbox/localization/relocate.py @@ -27,8 +27,8 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, - _default_base, - _known_compartments, + _base_met, + _check_unique_keys, apply_assignment, ) from raven_toolbox.localization.curation import _is_impermeant @@ -117,12 +117,8 @@ def relocate_reactions( removes a transport of species ``b`` at a reaction's old compartment when nothing else uses ``b`` there any more. The transport repair only touches species of the moved reactions. """ - if base_metabolite is not None: - base = base_metabolite - else: - base = _default_base(_known_compartments( - model, universal, [default_compartment], decisions.values(), - *proposal.placements.values())) + base = base_metabolite if base_metabolite is not None else _base_met + _check_unique_keys(model.metabolites, base) floor = float(min_growth if min_growth is not None else proposal.min_growth) warnings: list[str] = [] diff --git a/src/raven_toolbox/manipulation/compartments.py b/src/raven_toolbox/manipulation/compartments.py index 9cfe333..a6466e8 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -18,6 +18,10 @@ # Compartments produced by merge_compartments (RAVEN uses 's' for "system"). _MERGED_COMPARTMENT = "s" +#: Key in ``Metabolite.notes`` where :func:`merge_compartments` stores the base id of a merged +#: metabolite, read by the default metabolite key of :mod:`raven_toolbox.localization`. +MERGED_BASE_ID_NOTE = "merged_base_id" + def merge_compartments( model: cobra.Model, @@ -64,7 +68,23 @@ def merge_compartments( already had a single metabolite before merging is an exchange (``=> glucose``) and is always kept, as in RAVEN; dropping those would remove the model's medium and its biomass reaction. Reactions that become identical after merging are deduplicated (one survives). The - objective is carried over to the returned model. + objective is carried over to the returned model. ``deleted_single_met_reactions`` lists the + reactions whose metabolites all cancelled and, with ``drop_single_metabolite_reactions=True``, + those reduced to one metabolite. RAVEN's ``deletedRxns`` lists only the latter, so it is empty + with ``deleteRxnsWithOneMet`` off. + + Each merged metabolite stores its base id in ``notes["merged_base_id"]`` + (:data:`MERGED_BASE_ID_NOTE`): the id of the kept copy without the suffix of the compartment it + came from (``ACCOA`` for ``ACCOA_c``). The default metabolite key of + :func:`~raven_toolbox.localization.assign_compartments`, + :func:`~raven_toolbox.localization.apply_assignment`, + :func:`~raven_toolbox.localization.relocate_reactions` and + :func:`~raven_toolbox.localization.curation_priority` reads it, so the merged ``ACCOA_c`` in + ``s`` and ``ACCOA_m`` of a universal model are the same species there. When two merged + metabolites would get the same base id (a ``base_metabolite`` that keeps ``X_c`` and ``X_m`` + apart), both store their own id instead, so they stay distinct species. The note survives + ``model.copy()``, cobra's JSON and SBML files and :func:`~raven_toolbox.io.write_yaml_model` + with :func:`~raven_toolbox.io.read_yaml_model`. """ out = model.copy() key = base_metabolite if base_metabolite is not None else _base_id @@ -79,15 +99,22 @@ def merge_compartments( # 2. Build the merged metabolites and rewrite reactions. Each merged metabolite is a copy of one # member, chosen as the docstring states. Groups are disjoint and model ids are unique, so the # chosen ids are unique too. - canonical: dict[str, cobra.Metabolite] = {} + protos = {} for base, mets in new_to_old.items(): in_target = [m for m in mets if m.compartment == merged_id] - proto = min(in_target or mets, key=lambda m: m.id) + protos[base] = min(in_target or mets, key=lambda m: m.id) + base_id_count: dict[str, int] = {} + for proto in protos.values(): + base_id_count[_base_id(proto)] = base_id_count.get(_base_id(proto), 0) + 1 + canonical: dict[str, cobra.Metabolite] = {} + for base, proto in protos.items(): mid = proto.id new_met = cobra.Metabolite(mid, name=proto.name, compartment=merged_id, formula=proto.formula, charge=proto.charge) new_met.notes = copy.deepcopy(proto.notes or {}) new_met.annotation = copy.deepcopy(proto.annotation or {}) + stored = _base_id(proto) + new_met.notes[MERGED_BASE_ID_NOTE] = stored if base_id_count[stored] == 1 else mid new_met._bound = proto._bound canonical[base] = new_met diff --git a/tests/test_localization_certify.py b/tests/test_localization_certify.py index dc43bab..680f253 100644 --- a/tests/test_localization_certify.py +++ b/tests/test_localization_certify.py @@ -469,3 +469,104 @@ def test_apply_assignment_matches_a_merged_metabolite_to_the_universal_copy(): assert len(a_in_m) == 1 assert a_in_m[0] in out.reactions.R1.metabolites assert a_in_m[0] in out.reactions.U1.metabolites + + +def _draft_merged_into_c(): + """A draft merged into c by merge_compartments: the intermediate I had copies I_g and I_x only, + so the merged I keeps the id I_g in c.""" + from raven_toolbox.manipulation import merge_compartments + + o = cobra.Model("orig") + a = cobra.Metabolite("A_c", name="A", compartment="c") + b = cobra.Metabolite("B_c", name="B", compartment="c") + ig = cobra.Metabolite("I_g", name="I", compartment="g") + ix = cobra.Metabolite("I_x", name="I", compartment="x") + ex = cobra.Reaction("EX_A", lower_bound=-10, upper_bound=1000) + r1 = cobra.Reaction("r1", lower_bound=0, upper_bound=1000) + r2 = cobra.Reaction("r2", lower_bound=0, upper_bound=1000) + t = cobra.Reaction("tI", lower_bound=-1000, upper_bound=1000) + bio = cobra.Reaction("bio", lower_bound=0, upper_bound=1000) + o.add_reactions([ex, r1, r2, t, bio]) + ex.add_metabolites({a: -1}) + r1.add_metabolites({a: -1, ig: 1}) + r2.add_metabolites({ix: -1, b: 1}) + t.add_metabolites({ig: -1, ix: 1}) + bio.add_metabolites({b: -1}) + r1.gene_reaction_rule = "g1" + r2.gene_reaction_rule = "g2" + o.objective = "bio" + draft, _, _ = merge_compartments(o, merged_id="c", merged_name="cytosol") + assert draft.metabolites.I_g.compartment == "c" + scores = LocalizationScores(pd.DataFrame( + {"c": [0.1, 0.1], "m": [0.9, 0.0], "p": [0.0, 0.9], "g": [0.0, 0.0], "x": [0.0, 0.0]}, + index=["g1", "g2"])) + return draft, scores + + +def test_default_key_of_a_merged_metabolite_is_the_same_in_every_localization_function(): + """assign_compartments and apply_assignment both key the merged I_g as I, so the model built + from the proposal has every transport the proposal lists and grows.""" + from raven_toolbox.localization import curation_priority, relocate_reactions + + draft, scores = _draft_merged_into_c() + res = assign_compartments(draft, scores, ["r1", "r2"]) + assert res.certified + assert {b for b, _ in res.added_transports} == {"A", "I", "B"} + out = apply_assignment(draft, res) + assert len([r for r in out.reactions if r.id.startswith("tr_")]) == len(res.added_transports) + assert out.slim_optimize(error_value=0.0) > 1e-6 + moved = relocate_reactions(draft, res, {"r2": "m"}) + assert {b for b, _ in moved.transports_added + moved.transports_removed} <= {"A", "I", "B"} + curation_priority(draft, res, scores, check_essential=False) + + +def test_name_keyed_merge_keeps_distinct_species_distinct_in_apply_assignment(): + """X_c (name X) and X_m (name X mito) are two species; after a name-keyed merge both are in s + and apply_assignment makes a copy of each in m.""" + from raven_toolbox.manipulation import merge_compartments + + m = cobra.Model("m") + xc = cobra.Metabolite("X_c", name="X", compartment="c") + xm = cobra.Metabolite("X_m", name="X mito", compartment="m") + yc = cobra.Metabolite("Y_c", name="Y", compartment="c") + ra = cobra.Reaction("Ra", lower_bound=0, upper_bound=1000) + rb = cobra.Reaction("Rb", lower_bound=0, upper_bound=1000) + m.add_reactions([ra, rb]) + ra.add_metabolites({yc: -1, xc: 1}) + rb.add_metabolites({xm: -1, yc: 1}) + merged, _, _ = merge_compartments(m, base_metabolite=lambda x: x.name) + out = apply_assignment(merged, AssignmentProposal(placements={"Ra": ["m"], "Rb": ["m"]})) + in_m = sorted(x.name for x in out.metabolites if x.compartment == "m") + assert in_m == ["X", "X mito", "Y"] + + +def test_id_that_ends_in_another_compartment_id_keeps_its_own_key(): + """vitamin_a in c and vitamin_c in e stay two species when the universal model has a + compartment a: only the suffix of a metabolite's own compartment is stripped.""" + m = cobra.Model("m") + va = cobra.Metabolite("vitamin_a", name="retinol", compartment="c") + vc = cobra.Metabolite("vitamin_c", name="ascorbate", compartment="e") + z = cobra.Metabolite("Z_c", name="Z", compartment="c") + rv = cobra.Reaction("Rv", lower_bound=0, upper_bound=1000) + rw = cobra.Reaction("Rw", lower_bound=0, upper_bound=1000) + m.add_reactions([rv, rw]) + rv.add_metabolites({va: -1, z: 1}) + rw.add_metabolites({vc: -1, z: 1}) + u = cobra.Model("u") + u.add_metabolites([cobra.Metabolite("Q_a", name="Q", compartment="a")]) + out = apply_assignment(m, AssignmentProposal(placements={"Rv": ["m"], "Rw": ["m"]}), universal=u) + assert sorted(x.name for x in out.metabolites if x.compartment == "m") == ["Z", "ascorbate", "retinol"] + + +def test_two_metabolites_with_one_key_in_a_compartment_raise(): + m = cobra.Model("m") + h1 = cobra.Metabolite("h_c", name="H+", compartment="c") + h2 = cobra.Metabolite("h", name="proton", compartment="c") + r = cobra.Reaction("R", lower_bound=0, upper_bound=1000) + m.add_reactions([r]) + r.add_metabolites({h1: -1, h2: 1}) + with pytest.raises(ValueError, match="same key"): + apply_assignment(m, AssignmentProposal(placements={"R": ["m"]})) + scores = LocalizationScores(pd.DataFrame({"c": [0.1], "m": [0.9]}, index=["g1"])) + with pytest.raises(ValueError, match="same key"): + assign_compartments(m, scores, ["R"], biomass_reaction="R") diff --git a/tests/test_manipulation_compartments.py b/tests/test_manipulation_compartments.py index 0c71c9f..72bc8a6 100644 --- a/tests/test_manipulation_compartments.py +++ b/tests/test_manipulation_compartments.py @@ -259,11 +259,52 @@ def test_merge_compartments_keeps_the_annotation_and_notes_of_the_chosen_copy(): merged, _, _ = merge_compartments(m, base_metabolite=lambda mm: mm.name) kept = merged.metabolites.ACCOA_c assert dict(kept.annotation) == {"chebi": "CHEBI:2", "kegg.compound": "C00024"} - assert kept.notes == {"note": "cytosolic copy"} + assert kept.notes == {"note": "cytosolic copy", "merged_base_id": "ACCOA"} kept.annotation["chebi"] = "changed" assert m.metabolites.ACCOA_c.annotation["chebi"] == "CHEBI:2" +def test_merge_compartments_stores_the_base_id_of_each_merged_metabolite(): + """The base id is the kept copy's id without the suffix of the compartment it came from.""" + from raven_toolbox.manipulation.compartments import MERGED_BASE_ID_NOTE + + merged, _, _ = merge_compartments(_two_compartment_model()) + assert {x.id: x.notes[MERGED_BASE_ID_NOTE] for x in merged.metabolites} == {"A_c": "A", "B_c": "B"} + named, _, _ = merge_compartments(_name_keyed_model(), base_metabolite=lambda x: x.name) + assert {x.id: x.notes[MERGED_BASE_ID_NOTE] for x in named.metabolites} == {"s_1": "s_1", "s_3": "s_3"} + + +def test_merge_compartments_stores_the_own_id_when_two_base_ids_coincide(): + """X_c (name X) and X_m (name X mito) are kept apart by a name key; their base ids are both X, + so each stores its own id and the two stay distinct species.""" + from raven_toolbox.manipulation.compartments import MERGED_BASE_ID_NOTE + + m = cobra.Model("m") + m.add_metabolites([cobra.Metabolite("X_c", name="X", compartment="c"), + cobra.Metabolite("X_m", name="X mito", compartment="m"), + cobra.Metabolite("Y_c", name="Y", compartment="c")]) + merged, _, _ = merge_compartments(m, base_metabolite=lambda x: x.name) + assert {x.id: x.notes[MERGED_BASE_ID_NOTE] for x in merged.metabolites} == { + "X_c": "X_c", "X_m": "X_m", "Y_c": "Y"} + + +def test_merged_base_id_survives_copy_and_model_files(tmp_path): + from raven_toolbox.io import read_yaml_model, write_yaml_model + from raven_toolbox.manipulation.compartments import MERGED_BASE_ID_NOTE + + m = _two_compartment_model() + m.metabolites.A_c.notes = {"note": "free text"} + merged, _, _ = merge_compartments(m) + write_yaml_model(merged, tmp_path / "m.yml") + cobra.io.write_sbml_model(merged, str(tmp_path / "m.xml")) + for model in (merged.copy(), read_yaml_model(tmp_path / "m.yml"), + cobra.io.from_json(cobra.io.to_json(merged)), + cobra.io.read_sbml_model(str(tmp_path / "m.xml"))): + assert model.metabolites.A_c.notes[MERGED_BASE_ID_NOTE] == "A" + assert model.metabolites.B_c.notes[MERGED_BASE_ID_NOTE] == "B" + assert read_yaml_model(tmp_path / "m.yml").metabolites.A_c.notes["note"] == "free text" + + # ----------------------------------------------------------------- copy_to_compartment def test_copy_to_compartment_basic(): From 66f48046f4489e36e03268c6805dbde0705b14c9 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:31:52 +0100 Subject: [PATCH 20/35] read_yaml_model: read gene rules with parse_gr_rule cobra's model_from_dict read the rule G1&G2 or G3, written by write_yaml_model after batch_curate, as the genes G1, G2 and G3. The reader now gives each reaction the rule as parse_gr_rule reads it, so G1&G2 stays one gene; a rule that parse_gr_rule refuses is read by cobra as before. The batch_curate help states that cobra's JSON reader still splits such an id. --- CHANGELOG.md | 3 +++ src/raven_toolbox/curation/batch.py | 4 +++ src/raven_toolbox/io/yaml.py | 28 +++++++++++++++++++++ tests/test_curation.py | 18 +++++++++++++ tests/test_io_yaml.py | 39 +++++++++++++++++++++++++++++ 5 files changed, 92 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 76d43f3..9b02fd2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,9 @@ Milestones in the raven-toolbox port. For function-level status see * `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. cobra's JSON reader still splits `G1&G2`. * `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 diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index 0545e61..6163cca 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -179,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 diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 325d69a..87012b4 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -58,6 +58,7 @@ ec_data_from_yaml_sections, ec_data_to_yaml_sections, ) +from raven_toolbox.utils.gpr import parse_gr_rule # A dedicated instance for *writing* only (reading still goes through # _cobra_yaml, whose parsing is unaffected by any of this). Mutating @@ -152,6 +153,26 @@ def _capture_entry_fields(entries, fields): return captured +def _pop_gene_rule(entry) -> str | None: + """Remove and return the gene rule of a reaction entry (``None`` when it has none).""" + if not isinstance(entry, dict): + return None + rule = entry.pop("gene_reaction_rule", None) + return None if rule is None or str(rule).strip() == "" else str(rule) + + +def _set_gene_rule(rxn: cobra.Reaction, rule: str) -> None: + """Give ``rxn`` the rule as :func:`~raven_toolbox.utils.parse_gr_rule` reads it. A rule that + grammar refuses is handed to cobra as it is.""" + try: + gpr = parse_gr_rule(rule) + except ValueError: + rxn.gene_reaction_rule = rule + return + if gpr is not None: + rxn.gpr = gpr + + def read_yaml_model(path: str | Path) -> cobra.Model: """Read a RAVEN/cobrapy YAML model into a ``cobra.Model``. @@ -271,8 +292,15 @@ def model_from_yaml_data(raw: dict) -> cobra.Model: rxn_notes = _capture_entry_fields(raw.get("reactions", []), _RXN_FIELDS) gene_notes = _capture_entry_fields(raw.get("genes", []), _GENE_FIELDS) + # Gene rules are read with RAVEN's grammar (parse_gr_rule), in which G1&G2 is one gene id; + # cobra's model_from_dict would read it as G1 and G2. + rules = [_pop_gene_rule(entry) for entry in raw.get("reactions", []) or ()] + model = model_from_dict(raw) + for rxn, rule in zip(model.reactions, rules, strict=False): + if rule is not None: + _set_gene_rule(rxn, rule) for met, notes in zip(model.metabolites, met_notes, strict=False): met.notes = notes for rxn, notes in zip(model.reactions, rxn_notes, strict=False): diff --git a/tests/test_curation.py b/tests/test_curation.py index eca5aaa..654d0a0 100644 --- a/tests/test_curation.py +++ b/tests/test_curation.py @@ -1305,6 +1305,24 @@ def test_an_operator_glued_to_a_gene_id_stays_part_of_the_stored_gene(name, ctab assert not {"G1", "G2"} & {g.id for g in m.genes} +def test_a_glued_gene_id_stays_one_gene_through_the_yaml_files(tmp_path): + from raven_toolbox.io import read_yaml_model, write_yaml_model + + m = _base_model() + rxns_df = pd.DataFrame([{"rxnNames": "fresh", "grRules": "G1&G2 or G3"}]) + coeffs_df = pd.DataFrame([ + {"rxnNames": "fresh", "metNames": "ATP", "comps": "c", "coefficient": "-1"}, + {"rxnNames": "fresh", "metNames": "ADP", "comps": "c", "coefficient": "1"}]) + batch_curate(m, genes_df=_genes_df("G1&G2", "G3"), rxns_df=rxns_df, + rxns_coeffs_df=coeffs_df, rxn_id_prefix="r_") + write_yaml_model(m, tmp_path / "m.yml") + back = read_yaml_model(tmp_path / "m.yml") + rxn = _by_name(back, "fresh") + assert sorted(g.id for g in rxn.genes) == ["G1&G2", "G3"] + assert rxn.gene_reaction_rule == "G1&G2 or G3" + assert not {"G1", "G2"} & {g.id for g in back.genes} + + # --- header-only tables, atomicity of several rows, genes of new reactions, headers --------- _RXN_HEADER = "rxnIdx\trxnNames\tlb\tub\n" diff --git a/tests/test_io_yaml.py b/tests/test_io_yaml.py index d232567..c35a382 100644 --- a/tests/test_io_yaml.py +++ b/tests/test_io_yaml.py @@ -506,3 +506,42 @@ def test_only_the_stored_default_bound_is_written(tmp_path): text = out.read_text() assert "defaultLB" in text assert "defaultUB" not in text + + +@pytest.mark.parametrize("rule, genes, text", [ + ("G1&G2 or G3", ["G1&G2", "G3"], "G1&G2 or G3"), + ("(G1 and G2) or G3", ["G1", "G2", "G3"], "(G1 and G2) or G3"), + ("G1 && G2", ["G1", "G2"], "G1 and G2"), +]) +def test_read_yaml_model_reads_gene_rules_as_parse_gr_rule(tmp_path, rule, genes, text): + """RAVEN's parseGrRule reads G1&G2 as one gene id and && as and.""" + m = cobra.Model("t") + a = cobra.Metabolite("a_c", compartment="c") + r = cobra.Reaction("R1") + m.add_reactions([r]) + r.add_metabolites({a: -1}) + write_yaml_model(m, tmp_path / "m.yml") + data = (tmp_path / "m.yml").read_text().replace( + "- id: R1\n", f"- id: R1\n - gene_reaction_rule: {rule}\n") + (tmp_path / "m.yml").write_text(data) + back = read_yaml_model(tmp_path / "m.yml") + assert sorted(g.id for g in back.reactions.R1.genes) == genes + assert back.reactions.R1.gene_reaction_rule == text + + +def test_read_yaml_model_hands_a_rule_parse_gr_rule_refuses_to_cobra(tmp_path): + """Two genes without an operator are a malformed rule for parse_gr_rule; the file is still + read, with cobra's reading of that rule.""" + m = cobra.Model("t") + a = cobra.Metabolite("a_c", compartment="c") + r = cobra.Reaction("R1") + m.add_reactions([r]) + r.add_metabolites({a: -1}) + r.gene_reaction_rule = "G1 or G2" + write_yaml_model(m, tmp_path / "m.yml") + data = (tmp_path / "m.yml").read_text().replace("G1 or G2", "G1 G2 or G3") + (tmp_path / "m.yml").write_text(data) + expected = cobra.Reaction("x") + expected.gene_reaction_rule = "G1 G2 or G3" + back = read_yaml_model(tmp_path / "m.yml") + assert back.reactions.R1.gene_reaction_rule == expected.gene_reaction_rule From ea8fbd0c87dae946d11aefd8e3ffa77c535845df Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:31:52 +0100 Subject: [PATCH 21/35] run_supporth, build_ko_hmm, build_hmm_library: threads default to every CPU the process may use As run_blast and run_hmmsearch, and as RAVEN's runSuppOrth, which passes the cores assigned to MATLAB (was: the machine's CPU count minus one). --- CHANGELOG.md | 7 +++-- .../reconstruction/homology/supporth.py | 7 +++-- src/raven_toolbox/reconstruction/kegg/hmm.py | 12 +++++--- tests/test_reconstruction_kegg_hmm.py | 30 +++++++++++++++++++ tests/test_reconstruction_supporth.py | 11 +++++++ 5 files changed, 58 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9b02fd2..81c4a33 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,10 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased -* `reconstruction.homology.run_blast` and `run_diamond`, and `reconstruction.kegg.run_hmmsearch` - and `get_kegg_model_from_sequences`, use every CPU the process may run on when `threads` is - not given, as RAVEN's `getBlast`, `getDiamond` and `getKEGGModelForOrganism` do (was: the +* `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 diff --git a/src/raven_toolbox/reconstruction/homology/supporth.py b/src/raven_toolbox/reconstruction/homology/supporth.py index ecd1ac0..56d7bd3 100644 --- a/src/raven_toolbox/reconstruction/homology/supporth.py +++ b/src/raven_toolbox/reconstruction/homology/supporth.py @@ -27,6 +27,7 @@ import pandas as pd +from raven_toolbox.reconstruction.homology.blast import _available_cpus from raven_toolbox.reconstruction.homology.hits import HIT_COLUMNS SOURCE = "SuppOrth" @@ -331,7 +332,7 @@ def run_supporth( support_lv: int | None = None, identity_filter: float = 80.0, alignment: str = "diamond", - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, extra_args: Sequence[str] = (), supporth: str | os.PathLike | None = None, **reader_kwargs, @@ -348,7 +349,8 @@ def run_supporth( it, or if its identity is above ``identity_filter``. ``support_lv`` defaults to 3, or to the number of ``tools`` when fewer are selected. ``extra_args`` are appended to the command, for example ``["--resume"]`` or - ``["--oma", "path/to/Output"]``. Output is written to ``outdir``. + ``["--oma", "path/to/Output"]``. Output is written to ``outdir``. ``threads`` defaults to + every CPU the process may use, as in RAVEN's ``runSuppOrth``. SuppOrth requires protein ids to be unique over all proteomes. When two FASTA files share an id (for example two templates of the same species), every id is @@ -369,6 +371,7 @@ def run_supporth( "SuppOrth's predictors run on Linux and macOS only. Run SuppOrth under WSL or on " "a Linux host, and read suppOrth.tsv with hits_from_supporth()." ) + threads = _available_cpus() if threads is None else threads organism_id = str(organism_id) model_ids = [model_ids] if isinstance(model_ids, str) else [str(m) for m in model_ids] ref_fastas = [ref_fastas] if isinstance(ref_fastas, (str, Path)) else list(ref_fastas) diff --git a/src/raven_toolbox/reconstruction/kegg/hmm.py b/src/raven_toolbox/reconstruction/kegg/hmm.py index 0fe0ed1..df16624 100644 --- a/src/raven_toolbox/reconstruction/kegg/hmm.py +++ b/src/raven_toolbox/reconstruction/kegg/hmm.py @@ -38,6 +38,7 @@ from tqdm.contrib.logging import logging_redirect_tqdm from raven_toolbox.binaries import resolve_binary +from raven_toolbox.reconstruction.homology.blast import _available_cpus from raven_toolbox.reconstruction.kegg.taxonomy import organisms_in_domain logger = logging.getLogger(__name__) @@ -326,7 +327,7 @@ def build_ko_hmm( *, seq_identity: float = 0.9, parttree_residues: int | None = None, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, fast: bool = True, verbose: bool = False, cdhit: str | Path | None = None, @@ -347,8 +348,10 @@ def build_ko_hmm( simple residue-count cutoff. ``verbose`` logs (via the ``logging`` module, INFO/DEBUG) which tool is running for this KO, sequence counts at each stage, timings, and the tools' own - output. Returns ``out_hmm``. + output. ``threads`` (CD-HIT, MAFFT and hmmbuild) defaults to every CPU the process may + use. Returns ``out_hmm``. """ + threads = _available_cpus() if threads is None else threads ko_fasta = Path(ko_fasta) out_hmm = Path(out_hmm) label = out_hmm.stem @@ -431,7 +434,7 @@ def build_hmm_library( domain: str, seq_identity: float = 0.9, parttree_residues: int | None = None, - threads: int = max(1, (os.cpu_count() or 2) - 1), + threads: int | None = None, fast: bool = True, verbose: bool = False, progress: bool = False, @@ -456,7 +459,8 @@ def build_hmm_library( ``progress`` shows tqdm bars for the multi-FASTA and per-KO HMM passes (the headline ``N of M`` counter for the long build); it is independent of ``verbose`` (per-KO INFO logging). With both on, per-KO log lines are routed - through ``tqdm.write`` so they don't corrupt the bar. + through ``tqdm.write`` so they don't corrupt the bar. ``threads`` goes to + :func:`build_ko_hmm` and defaults to every CPU the process may use. """ out_dir = Path(out_dir) fasta_dir = out_dir / "fasta" diff --git a/tests/test_reconstruction_kegg_hmm.py b/tests/test_reconstruction_kegg_hmm.py index 48bba63..b1107c4 100644 --- a/tests/test_reconstruction_kegg_hmm.py +++ b/tests/test_reconstruction_kegg_hmm.py @@ -207,6 +207,36 @@ def fake_run(cmd, *, stdout_path=None): assert out.read_text() == "HMM\n" +def test_build_ko_hmm_threads_default_to_every_available_cpu(tmp_path, monkeypatch): + import os + + fasta = tmp_path / "K90001.fa" + fasta.write_text(">a\nMKV\n>b\nMRV\n") + cmds = [] + monkeypatch.setattr(os, "sched_getaffinity", lambda pid: set(range(7)), raising=False) + monkeypatch.setattr(os, "cpu_count", lambda: 16) + monkeypatch.setattr( + "raven_toolbox.reconstruction.kegg.hmm.resolve_binary", + lambda exe, binary=None: binary or exe, + ) + + def fake_run(cmd, *, stdout_path=None): + cmds.append(cmd) + if stdout_path is not None: + Path(stdout_path).write_text(">a\nMKV\n>b\nMRV\n") + if Path(cmd[0]).name == "cd-hit": + Path(cmd[cmd.index("-o") + 1]).write_text(">a\nMKV\n>b\nMRV\n") + if Path(cmd[0]).name == "hmmbuild": + Path(cmd[-2]).write_text("HMM\n") + return "" + + monkeypatch.setattr("raven_toolbox.reconstruction.kegg.hmm._run", fake_run) + build_ko_hmm(fasta, tmp_path / "K90001.hmm") + flags = {"cd-hit": "-T", "mafft": "--thread", "hmmbuild": "--cpu"} + assert {Path(c[0]).name: c[c.index(flags[Path(c[0]).name]) + 1] for c in cmds} == { + "cd-hit": "7", "mafft": "7", "hmmbuild": "7"} + + def test_build_ko_hmm_single_sequence_skips_align(tmp_path, monkeypatch): fasta = tmp_path / "K9.fa" fasta.write_text(">only\nMKV\n") diff --git a/tests/test_reconstruction_supporth.py b/tests/test_reconstruction_supporth.py index 7e2e379..6ff1a06 100644 --- a/tests/test_reconstruction_supporth.py +++ b/tests/test_reconstruction_supporth.py @@ -305,6 +305,17 @@ def test_run_supporth_command(fake_run, fastas, tmp_path): assert not (tmp_path / "out" / "prefixed_inputs").exists() +def test_run_supporth_threads_default_to_every_available_cpu(fake_run, fastas, tmp_path, monkeypatch): + """As RAVEN's runSuppOrth, which passes feature('numcores')'s assigned count.""" + import os + + monkeypatch.setattr(os, "sched_getaffinity", lambda pid: set(range(7)), raising=False) + monkeypatch.setattr(os, "cpu_count", lambda: 16) + run_supporth("bug", fastas["bug"], "templateGEM", fastas["templateGEM"], tmp_path / "out", + supporth="suppOrth") + assert fake_run["cmd"][fake_run["cmd"].index("-t") + 1] == "7" + + def test_run_supporth_default_support(fake_run, fastas, tmp_path): run_supporth("bug", fastas["bug"], "templateGEM", fastas["templateGEM"], tmp_path / "out", supporth="suppOrth") From a3b9c1294ef108ef8dbe9abc104417b8a379727f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:31:52 +0100 Subject: [PATCH 22/35] tests: hmmsearch thread default compares the --cpu value The test patched a helper that only exists once the default is in place, so without it the test stopped with an AttributeError. It now fixes the CPUs the process may use and checks the --cpu value, also for get_kegg_model_from_sequences without threads. --- tests/test_reconstruction_kegg_query.py | 29 +++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/tests/test_reconstruction_kegg_query.py b/tests/test_reconstruction_kegg_query.py index 23a7d85..6992c3c 100644 --- a/tests/test_reconstruction_kegg_query.py +++ b/tests/test_reconstruction_kegg_query.py @@ -176,14 +176,15 @@ def test_prune_orthology_keeps_only_matched_kos(): assert full.reactions.get_by_id("R1").annotation["kegg.orthology"] == ["K1", "K2", "K3"] -def test_hmmsearch_threads_default_to_every_available_cpu(tmp_path, monkeypatch): - """As RAVEN's getKEGGModelForOrganism, which passes feature('numcores')'s assigned count.""" +def _fake_hmmsearch(monkeypatch, calls): + """Make hmmsearch a stub that records its command; the process may run on 7 of 16 CPUs.""" + import os import subprocess import raven_toolbox.reconstruction.kegg.query as query - calls = [] - monkeypatch.setattr(query, "_available_cpus", lambda: 7) + monkeypatch.setattr(os, "sched_getaffinity", lambda pid: set(range(7)), raising=False) + monkeypatch.setattr(os, "cpu_count", lambda: 16) monkeypatch.setattr(query, "resolve_binary", lambda name, binary=None: name) def fake_run(cmd, **kwargs): @@ -192,7 +193,27 @@ def fake_run(cmd, **kwargs): return subprocess.CompletedProcess(cmd, 0, "", "") monkeypatch.setattr(query.subprocess, "run", fake_run) + return query + + +def test_hmmsearch_threads_default_to_every_available_cpu(tmp_path, monkeypatch): + """As RAVEN's getKEGGModelForOrganism, which passes feature('numcores')'s assigned count.""" + calls: list = [] + query = _fake_hmmsearch(monkeypatch, calls) library = tmp_path / "lib.hmm" library.write_text("HMMER3/f\n") query.run_hmmsearch("org.faa", library) assert calls[0][calls[0].index("--cpu") + 1] == "7" + query.run_hmmsearch("org.faa", library, threads=3) + assert calls[1][calls[1].index("--cpu") + 1] == "3" + + +def test_kegg_model_from_sequences_passes_the_default_threads_on(tmp_path, monkeypatch, reference_and_tables): + calls: list = [] + query = _fake_hmmsearch(monkeypatch, calls) + model_ref, tables = reference_and_tables + library = tmp_path / "lib.hmm" + library.write_text("HMMER3/f\n") + query.get_kegg_model_from_sequences("org.faa", model_ref, tables["ko_reaction"], library, + rxn_flags=tables["rxn_flags"]) + assert calls[0][calls[0].index("--cpu") + 1] == "7" From 0535635848380d2962b14fd8560d89b64e05db1e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:11:59 +0100 Subject: [PATCH 23/35] Keep RAVEN's reading of a gene rule when a function copies or builds it utils.gpr.set_gene_rule copies a GPR object, or reads a rule string with parse_gr_rule, so G1&G2 stays one gene id. merge_compartments, copy_to_compartment, apply_assignment, convert_to_irreversible, merge_models, add_reactions_from_model, expand_model, change_gene_reaction_rules, add_reactions_from_equations, fill_tasks, remove_low_score_genes, the homology and KEGG model builders and read_matlab_model set rules through it instead of assigning the string to gene_reaction_rule, which cobra read with its own grammar. --- CHANGELOG.md | 11 +- scripts/benchmark_certified_multiorg.py | 3 +- src/raven_toolbox/init/genes.py | 4 +- src/raven_toolbox/init/taskfill.py | 3 +- src/raven_toolbox/io/mat.py | 3 +- src/raven_toolbox/io/yaml.py | 16 +- src/raven_toolbox/localization/assign.py | 5 +- src/raven_toolbox/manipulation/add.py | 12 +- src/raven_toolbox/manipulation/change.py | 12 +- .../manipulation/compartments.py | 6 +- src/raven_toolbox/manipulation/expand.py | 4 +- .../manipulation/irreversible.py | 4 +- src/raven_toolbox/manipulation/merge.py | 4 +- src/raven_toolbox/manipulation/transfer.py | 5 +- .../reconstruction/homology/homology.py | 3 +- .../reconstruction/kegg/assemble.py | 4 +- src/raven_toolbox/utils/gpr.py | 24 +++ tests/test_gene_rule_copy.py | 163 ++++++++++++++++++ 18 files changed, 244 insertions(+), 42 deletions(-) create mode 100644 tests/test_gene_rule_copy.py diff --git a/CHANGELOG.md b/CHANGELOG.md index dd7fc6b..b291b00 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -54,7 +54,16 @@ Milestones in the raven-toolbox port. For function-level status see `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. cobra's JSON reader still splits `G1&G2`. + `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. * `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 diff --git a/scripts/benchmark_certified_multiorg.py b/scripts/benchmark_certified_multiorg.py index 2e8cf14..f472fb8 100644 --- a/scripts/benchmark_certified_multiorg.py +++ b/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) diff --git a/src/raven_toolbox/init/genes.py b/src/raven_toolbox/init/genes.py index fa78a08..a964d2c 100644 --- a/src/raven_toolbox/init/genes.py +++ b/src/raven_toolbox/init/genes.py @@ -16,7 +16,7 @@ import cobra from cobra.manipulation import remove_genes -from raven_toolbox.utils.gpr import resolve_aggregators +from raven_toolbox.utils.gpr import resolve_aggregators, set_gene_rule def _prune(node, scores, iso, cplx, ties: list) -> tuple[str | None, float | None]: @@ -100,7 +100,7 @@ def remove_low_score_genes( if ties and tied_reactions is not None: tied_reactions.append(rxn.id) if pruned is not None: - rxn.gene_reaction_rule = pruned + set_gene_rule(rxn, pruned) used = {g.id for rxn in out.reactions for g in rxn.genes} removed = sorted(g.id for g in out.genes if g.id not in used) diff --git a/src/raven_toolbox/init/taskfill.py b/src/raven_toolbox/init/taskfill.py index 0d0bd7f..b502df3 100644 --- a/src/raven_toolbox/init/taskfill.py +++ b/src/raven_toolbox/init/taskfill.py @@ -31,6 +31,7 @@ apply_task_constraints, task_name_maps, ) +from raven_toolbox.utils.gpr import set_gene_rule _FILL_TIME_LIMIT = 300.0 # RAVEN ftINITFillGapsMILP TimeLimit (per task) _FILL_SEED = 26 # RAVEN ftINITFillGapsMILP Seed @@ -115,7 +116,7 @@ def _add_reference_reactions( target.add_reactions([new for new, _ in built]) for new, src in built: new.add_metabolites({target.metabolites.get_by_id(m.id): c for m, c in src.metabolites.items()}) - new.gene_reaction_rule = src.gene_reaction_rule + set_gene_rule(new, src.gpr) def _set_fill_solver(model: cobra.Model, time_limit: float | None, seed: int) -> None: diff --git a/src/raven_toolbox/io/mat.py b/src/raven_toolbox/io/mat.py index ddfd046..8f8678a 100644 --- a/src/raven_toolbox/io/mat.py +++ b/src/raven_toolbox/io/mat.py @@ -34,6 +34,7 @@ from raven_toolbox.io.metadata import get_model_metadata from raven_toolbox.io.yaml import _default_bounds +from raven_toolbox.utils.gpr import set_gene_rule # Field order RAVEN's exporter emits. Writing in this order means a model # round-tripped through here keeps the same field order as one written by @@ -621,7 +622,7 @@ def read_matlab_model(path: str | Path, *, varname: str | None = None) -> cobra. {metabolites[i]: float(v) for i, v in zip(column.row, column.data, strict=False)} ) if gr_rules[j]: - rxn.gene_reaction_rule = gr_rules[j] + set_gene_rule(rxn, gr_rules[j]) if subsystems is not None: entry = np.asarray(subsystems).reshape(-1, order="F") if j < entry.size: diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index b2c92bf..22179cf 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -59,7 +59,7 @@ ec_data_from_yaml_sections, ec_data_to_yaml_sections, ) -from raven_toolbox.utils.gpr import parse_gr_rule +from raven_toolbox.utils.gpr import set_gene_rule # A dedicated instance for *writing* only (reading still goes through # _cobra_yaml, whose parsing is unaffected by any of this). Mutating @@ -167,18 +167,6 @@ def _pop_gene_rule(entry) -> str | None: return None if rule is None or str(rule).strip() == "" else str(rule) -def _set_gene_rule(rxn: cobra.Reaction, rule: str) -> None: - """Give ``rxn`` the rule as :func:`~raven_toolbox.utils.parse_gr_rule` reads it. A rule that - grammar refuses is handed to cobra as it is.""" - try: - gpr = parse_gr_rule(rule) - except ValueError: - rxn.gene_reaction_rule = rule - return - if gpr is not None: - rxn.gpr = gpr - - def read_yaml_model(path: str | Path) -> cobra.Model: """Read a RAVEN/cobrapy YAML model into a ``cobra.Model``. @@ -306,7 +294,7 @@ def model_from_yaml_data(raw: dict) -> cobra.Model: for rxn, rule in zip(model.reactions, rules, strict=False): if rule is not None: - _set_gene_rule(rxn, rule) + set_gene_rule(rxn, rule) for met, notes in zip(model.metabolites, met_notes, strict=False): met.notes = notes for rxn, notes in zip(model.reactions, rxn_notes, strict=False): diff --git a/src/raven_toolbox/localization/assign.py b/src/raven_toolbox/localization/assign.py index acf879f..29fb56f 100644 --- a/src/raven_toolbox/localization/assign.py +++ b/src/raven_toolbox/localization/assign.py @@ -26,6 +26,7 @@ import cobra from raven_toolbox.manipulation.compartments import MERGED_BASE_ID_NOTE +from raven_toolbox.utils.gpr import set_gene_rule __all__ = ["AssignmentProposal", "GrowthCondition", "apply_assignment"] @@ -210,7 +211,7 @@ def _add_universal_reaction(model, urxn, resolve): # copy, leaving the gap-fill reaction a disconnected island. Keying on (base, compartment) lands on # the species already there. new.add_metabolites({resolve(m, m.compartment): coeff for m, coeff in urxn.metabolites.items()}) - new.gene_reaction_rule = urxn.gene_reaction_rule + set_gene_rule(new, urxn.gpr) def _move_reaction(rxn, compartment, resolve): @@ -229,7 +230,7 @@ def _duplicate_reaction(model, rxn, compartment, resolve): lower_bound=rxn.lower_bound, upper_bound=rxn.upper_bound) model.add_reactions([dup]) dup.add_metabolites({resolve(m, compartment): coeff for m, coeff in rxn.metabolites.items()}) - dup.gene_reaction_rule = rxn.gene_reaction_rule + set_gene_rule(dup, rxn.gpr) def _add_transport(model, base_key, compartment, default_compartment, index, resolve, n): diff --git a/src/raven_toolbox/manipulation/add.py b/src/raven_toolbox/manipulation/add.py index a76bc7b..b546b62 100644 --- a/src/raven_toolbox/manipulation/add.py +++ b/src/raven_toolbox/manipulation/add.py @@ -9,7 +9,7 @@ and reversibility arrows (``<=>``, ``-->``, ``=>``) and creates unknown metabolites — but only matching metabolites **by ID**, and it leaves new metabolites with ``compartment=None``. -* assigning ``reaction.gene_reaction_rule`` auto-creates ``Gene`` objects. +* setting a reaction's GPR auto-creates ``Gene`` objects. So this port keeps only the parts cobra lacks: @@ -33,8 +33,9 @@ import cobra from cobra import Metabolite, Reaction -from cobra.core.gene import GPR, Gene +from cobra.core.gene import Gene +from raven_toolbox.utils.gpr import rule_to_gpr, set_gene_rule from raven_toolbox.utils.parse import parse_name_comp, subsystem_to_str # Reversibility arrows. ``<=>`` must be tried before ``=>`` (it contains it). @@ -284,6 +285,7 @@ def add_reactions_from_equations( * ``name`` — reaction name. * ``bounds`` — ``(lower, upper)`` tuple; overrides the arrow. * ``gene_reaction_rule`` — GPR string. + Read with :func:`~raven_toolbox.utils.parse_gr_rule`, RAVEN's grammar. * ``subsystem`` — subsystem name. mets_by How bare equation tokens (without ``[comp]``) are matched: @@ -361,7 +363,7 @@ def add_reactions_from_equations( # iteration order of the sets cobra builds from each rule. new_genes = sorted( {g for spec in reactions if spec.get("gene_reaction_rule") - for g in GPR.from_string(spec["gene_reaction_rule"]).genes} - known_genes + for g in rule_to_gpr(spec["gene_reaction_rule"]).genes} - known_genes ) for gene_id in new_genes: gene = Gene(gene_id) @@ -374,13 +376,13 @@ def add_reactions_from_equations( rule = spec.get("gene_reaction_rule", "") if rule: if not allow_new_genes: - missing = sorted(set(GPR.from_string(rule).genes) - known_genes) + missing = sorted(set(rule_to_gpr(rule).genes) - known_genes) if missing: raise ValueError( f"Reaction {rxn.id!r} references genes not in the model: " f"{missing}. Set allow_new_genes=True or add them first." ) - rxn.gene_reaction_rule = rule + set_gene_rule(rxn, rule) known_genes.update(gene.id for gene in rxn.genes) return added diff --git a/src/raven_toolbox/manipulation/change.py b/src/raven_toolbox/manipulation/change.py index 01ab3f3..6718fbc 100644 --- a/src/raven_toolbox/manipulation/change.py +++ b/src/raven_toolbox/manipulation/change.py @@ -19,6 +19,7 @@ from cobra import Reaction from raven_toolbox.manipulation.add import _build_met_index, _stoichiometry +from raven_toolbox.utils.gpr import set_gene_rule __all__ = ["change_reaction_equations", "change_gene_reaction_rules"] @@ -90,10 +91,11 @@ def change_gene_reaction_rules( replace: bool = True, ) -> list[Reaction]: """Set or append gene-reaction rules on existing reactions. - cobra already does the heavy lifting on assignment to - ``reaction.gene_reaction_rule``: it auto-creates any new ``Gene`` objects and - normalises the rule. So the value here is batching plus RAVEN's ``replace`` - option to **append** rather than overwrite. + + Rules are read with :func:`~raven_toolbox.utils.parse_gr_rule`, the grammar of RAVEN's + ``parseGrRule`` (``G1&G2`` is one gene id); a rule that grammar refuses is read with + cobra's parser. New genes are created in the model. ``replace=False`` appends the rule + as an isozyme, as RAVEN's ``changeGrRules`` with ``replace`` false does. Parameters ---------- @@ -122,7 +124,7 @@ def change_gene_reaction_rules( else: new_rule = f"({rxn.gene_reaction_rule}) or ({rule})" - rxn.gene_reaction_rule = new_rule # cobra creates genes + normalises + set_gene_rule(rxn, new_rule) changed.append(rxn) return changed diff --git a/src/raven_toolbox/manipulation/compartments.py b/src/raven_toolbox/manipulation/compartments.py index a6466e8..d67e078 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -15,6 +15,8 @@ import cobra +from raven_toolbox.utils.gpr import set_gene_rule + # Compartments produced by merge_compartments (RAVEN uses 's' for "system"). _MERGED_COMPARTMENT = "s" @@ -162,7 +164,7 @@ def merge_compartments( new_r = cobra.Reaction(r.id, name=r.name, lower_bound=r.lower_bound, upper_bound=r.upper_bound) new_r.add_metabolites({merged.metabolites.get_by_id(mid): c for mid, c in stoich.items()}) - new_r.gene_reaction_rule = r.gene_reaction_rule + set_gene_rule(new_r, r.gpr) if r.subsystem: new_r.subsystem = r.subsystem new_r.notes = dict(r.notes or {}) @@ -238,7 +240,7 @@ def copy_to_compartment( new_r = cobra.Reaction(new_id, name=src.name, lower_bound=src.lower_bound, upper_bound=src.upper_bound) new_r.add_metabolites(new_stoich) - new_r.gene_reaction_rule = src.gene_reaction_rule + set_gene_rule(new_r, src.gpr) if src.subsystem: new_r.subsystem = src.subsystem new_r.notes = dict(src.notes or {}) diff --git a/src/raven_toolbox/manipulation/expand.py b/src/raven_toolbox/manipulation/expand.py index 2bd98c1..fe77d56 100644 --- a/src/raven_toolbox/manipulation/expand.py +++ b/src/raven_toolbox/manipulation/expand.py @@ -15,7 +15,7 @@ import cobra from cobra.core.gene import GPR -from raven_toolbox.utils.gpr import parse_gr_rule +from raven_toolbox.utils.gpr import parse_gr_rule, set_gene_rule __all__ = ["expand_model", "gpr_to_dnf"] @@ -123,7 +123,7 @@ def expand_model(model: cobra.Model) -> list[str]: new_rxn.upper_bound = original_rxn.upper_bound new_rxn.add_metabolites(dict(original_rxn.metabolites.items())) new_rxn.subsystem = original_rxn.subsystem - new_rxn.gene_reaction_rule = " and ".join(clause) + set_gene_rule(new_rxn, " and ".join(clause)) # Propagate per-reaction metadata (notably ec-code / annotations) # so downstream functions see the same annotations on expanded # reactions as on the original. Deep-copy so siblings are independent. diff --git a/src/raven_toolbox/manipulation/irreversible.py b/src/raven_toolbox/manipulation/irreversible.py index 8e6f8b3..3283e97 100644 --- a/src/raven_toolbox/manipulation/irreversible.py +++ b/src/raven_toolbox/manipulation/irreversible.py @@ -9,6 +9,8 @@ import cobra +from raven_toolbox.utils.gpr import set_gene_rule + def convert_to_irreversible( model: cobra.Model, rxns: list[str] | None = None @@ -76,7 +78,7 @@ def convert_to_irreversible( rev_rxn.lower_bound = 0.0 rev_rxn.upper_bound = -original_lb rev_rxn.add_metabolites({m: -c for m, c in rxn.metabolites.items()}) - rev_rxn.gene_reaction_rule = rxn.gene_reaction_rule + set_gene_rule(rev_rxn, rxn.gpr) rev_rxn.subsystem = rxn.subsystem rev_rxn.annotation = copy.deepcopy(rxn.annotation) rev_rxn.notes = copy.deepcopy(rxn.notes) diff --git a/src/raven_toolbox/manipulation/merge.py b/src/raven_toolbox/manipulation/merge.py index 8366562..504e014 100644 --- a/src/raven_toolbox/manipulation/merge.py +++ b/src/raven_toolbox/manipulation/merge.py @@ -23,6 +23,8 @@ import cobra from cobra import Gene, Metabolite, Model, Reaction +from raven_toolbox.utils.gpr import set_gene_rule + def _unique_id(existing, base: str, suffix: str) -> str: """Return base, or base_suffix (then base_suffix_2, ...) if it collides.""" @@ -154,7 +156,7 @@ def ensure_metabolite(src: Metabolite, origin: str) -> Metabolite: {ensure_metabolite(m, origin): coef for m, coef in rxn.metabolites.items()} ) if rxn.gene_reaction_rule: - new_rxn.gene_reaction_rule = rxn.gene_reaction_rule + set_gene_rule(new_rxn, rxn.gpr) new_rxn.annotation = copy.deepcopy(rxn.annotation) new_rxn.notes = copy.deepcopy(rxn.notes) if track_origin: diff --git a/src/raven_toolbox/manipulation/transfer.py b/src/raven_toolbox/manipulation/transfer.py index 4f744c0..c2570e4 100644 --- a/src/raven_toolbox/manipulation/transfer.py +++ b/src/raven_toolbox/manipulation/transfer.py @@ -18,6 +18,7 @@ from cobra import Metabolite, Reaction from raven_toolbox.manipulation.add import _new_met_id +from raven_toolbox.utils.gpr import set_gene_rule def _name_comp(met: Metabolite) -> str: @@ -72,7 +73,7 @@ def add_reactions_from_model( if genes is False: rules = [""] * len(source_rxns) elif genes is True: - rules = [r.gene_reaction_rule for r in source_rxns] + rules = [r.gpr for r in source_rxns] elif isinstance(genes, str): rules = [genes] * len(source_rxns) else: @@ -135,7 +136,7 @@ def add_reactions_from_model( {draft_by_name[_name_comp(met)]: coef for met, coef in srx.metabolites.items()} ) if rule: - rxn.gene_reaction_rule = rule + set_gene_rule(rxn, rule) rxn.annotation = copy.deepcopy(srx.annotation) notes = copy.deepcopy(srx.notes) if note is not None: diff --git a/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index e8aea82..7ce6abb 100644 --- a/src/raven_toolbox/reconstruction/homology/homology.py +++ b/src/raven_toolbox/reconstruction/homology/homology.py @@ -22,6 +22,7 @@ from raven_toolbox.manipulation.merge import merge_models from raven_toolbox.reconstruction.homology.hits import validate_hits +from raven_toolbox.utils.gpr import set_gene_rule # cobra's GPR parser splits a gene id at whitespace, parentheses, "|" and "&". _GPR_SPLIT = re.compile(r"[\s()|&]") @@ -498,7 +499,7 @@ def _transfer(model_by_id, order, ortho, model_for, model_ids, complex_policy, m.remove_reactions([r for r in m.reactions if r.id not in keep], remove_orphans=True) for rid, gpr in keep.items(): r = m.reactions.get_by_id(rid) - r.gene_reaction_rule = gpr + set_gene_rule(r, gpr) if preserve_notes: r.notes = {**r.notes, "homology_source": mid} else: diff --git a/src/raven_toolbox/reconstruction/kegg/assemble.py b/src/raven_toolbox/reconstruction/kegg/assemble.py index 5c86775..44d0c84 100644 --- a/src/raven_toolbox/reconstruction/kegg/assemble.py +++ b/src/raven_toolbox/reconstruction/kegg/assemble.py @@ -11,6 +11,8 @@ import cobra import pandas as pd +from raven_toolbox.utils.gpr import set_gene_rule + _DOMAINS = {"eukaryotes", "prokaryotes"} @@ -93,7 +95,7 @@ def assemble_model_from_ko_genes( [r for r in model.reactions if r.id not in keep], remove_orphans=True ) for rid, genes in gpr_map.items(): - model.reactions.get_by_id(rid).gene_reaction_rule = " or ".join(genes) + set_gene_rule(model.reactions.get_by_id(rid), " or ".join(genes)) if prune_orthology: # Restrict each gene-backed reaction's kegg.orthology to the KOs that # matched a gene, preserving the reference annotation's order. (The diff --git a/src/raven_toolbox/utils/gpr.py b/src/raven_toolbox/utils/gpr.py index 5b88671..a2f3cb3 100644 --- a/src/raven_toolbox/utils/gpr.py +++ b/src/raven_toolbox/utils/gpr.py @@ -212,6 +212,30 @@ def expr(pos: int) -> tuple[ast.expr, int]: return GPR(ast.Expression(body=body)) +def rule_to_gpr(rule: str | None) -> GPR: + """The :class:`~cobra.core.gene.GPR` of a rule string as :func:`parse_gr_rule` reads it. + + An empty or ``None`` rule gives an empty GPR. A rule that grammar refuses is read with + cobra's ``GPR.from_string`` instead. + """ + try: + gpr = parse_gr_rule(rule) + except ValueError: + return GPR.from_string(rule) + return gpr if gpr is not None else GPR() + + +def set_gene_rule(rxn: cobra.Reaction, rule: GPR | str | None) -> None: + """Give ``rxn`` a gene rule in RAVEN's grammar. + + A :class:`~cobra.core.gene.GPR` is copied as it is; a string is read with + :func:`rule_to_gpr`. Assigning a string to ``rxn.gene_reaction_rule`` would read it + with cobra's grammar, which splits a gene id such as ``G1&G2`` into ``G1`` and ``G2``; + every function in this package that copies or builds a rule goes through here. + """ + rxn.gpr = rule.copy() if isinstance(rule, GPR) else rule_to_gpr(rule) + + @dataclass(frozen=True) class GPRIssue: """A reaction whose GPR is flagged by the linter. diff --git a/tests/test_gene_rule_copy.py b/tests/test_gene_rule_copy.py new file mode 100644 index 0000000..de64d06 --- /dev/null +++ b/tests/test_gene_rule_copy.py @@ -0,0 +1,163 @@ +"""Every function that copies or builds a gene rule keeps RAVEN's reading of it. + +In RAVEN's grammar ``G1&G2`` is one gene id (``&`` is an operator only as a whole word); +cobra's string parser would split it into ``G1`` and ``G2``. +""" +import cobra +import pytest + +from raven_toolbox.io.yaml import read_yaml_model, write_yaml_model +from raven_toolbox.utils import parse_gr_rule + +RULE = "G1&G2 or G3" +GENES = ["G1&G2", "G3"] + + +def _model(tmp_path): + """A model whose R1 has RAVEN's rule ``G1&G2 or G3``, read back from YAML.""" + m = cobra.Model("y") + a = cobra.Metabolite("A_c", name="A", compartment="c") + b = cobra.Metabolite("B_c", name="B", compartment="c") + m.add_metabolites([a, b]) + r = cobra.Reaction("R1", lower_bound=-10, upper_bound=10) + m.add_reactions([r]) + r.add_metabolites({a: -1, b: 1}) + r.gpr = parse_gr_rule(RULE) + m.compartments = {"c": "cytosol", "m": "mito"} + path = tmp_path / "m.yml" + write_yaml_model(m, path) + return read_yaml_model(path) + + +def _merge_compartments(m, tmp_path): + from raven_toolbox.manipulation import merge_compartments + return merge_compartments(m)[0], "R1" + + +def _copy_to_compartment(m, tmp_path): + from raven_toolbox.manipulation import copy_to_compartment + out, new, _ = copy_to_compartment(m, ["R1"], "m") + return out, new[0] + + +def _apply_assignment_duplicate(m, tmp_path): + from raven_toolbox.localization.assign import AssignmentProposal, apply_assignment + return apply_assignment(m, AssignmentProposal(placements={"R1": ["c", "m"]})), "R1_m" + + +def _apply_assignment_universal(m, tmp_path): + from raven_toolbox.localization.assign import AssignmentProposal, apply_assignment + universal = m.copy() + universal.reactions.R1.id = "U1" + universal.repair() + m.remove_reactions(["R1"], remove_orphans=True) + return apply_assignment(m, AssignmentProposal(added_reactions=["U1"]), + universal=universal), "U1" + + +def _convert_to_irreversible(m, tmp_path): + from raven_toolbox.manipulation import convert_to_irreversible + convert_to_irreversible(m) + return m, "R1_REV" + + +def _merge_models(m, tmp_path): + from raven_toolbox.manipulation import merge_models + return merge_models([m, cobra.Model("z")]), "R1" + + +def _task_fill_reference(m, tmp_path): + from raven_toolbox.init.taskfill import _add_reference_reactions + target = cobra.Model("t") + _add_reference_reactions(target, m, ["R1"]) + return target, "R1" + + +def _add_reactions_from_model(m, tmp_path): + from raven_toolbox.manipulation import add_reactions_from_model + target = cobra.Model("t") + target.add_metabolites([cobra.Metabolite("A_c", name="A", compartment="c"), + cobra.Metabolite("B_c", name="B", compartment="c")]) + add_reactions_from_model(target, m, ["R1"], genes=True) + return target, "R1" + + +def _matlab_round_trip(m, tmp_path): + from raven_toolbox.io.mat import read_matlab_model, write_matlab_model + path = tmp_path / "m.mat" + write_matlab_model(m, path) + return read_matlab_model(path), "R1" + + +def _change_rules(m, tmp_path): + from raven_toolbox.manipulation import change_gene_reaction_rules + change_gene_reaction_rules(m, {"R1": RULE}) + return m, "R1" + + +def _change_rules_append(m, tmp_path): + from raven_toolbox.manipulation import change_gene_reaction_rules + m.reactions.R1.gpr = parse_gr_rule("G1&G2") + change_gene_reaction_rules(m, {"R1": "G3"}, replace=False) + return m, "R1" + + +def _add_from_equations(m, tmp_path): + from raven_toolbox.manipulation import add_reactions_from_equations + add_reactions_from_equations(m, [{"id": "R9", "equation": "A_c => B_c", + "gene_reaction_rule": RULE}]) + return m, "R9" + + +def _remove_low_score_genes(m, tmp_path): + from raven_toolbox.init.genes import remove_low_score_genes + m.reactions.R1.gpr = parse_gr_rule("G1&G2 or G3 or G4") + out, _ = remove_low_score_genes(m, {"G1&G2": 1.0, "G3": 1.0, "G4": -5.0}) + return out, "R1" + + +@pytest.mark.parametrize("operation", [ + _merge_compartments, _copy_to_compartment, _apply_assignment_duplicate, + _apply_assignment_universal, _convert_to_irreversible, _merge_models, + _task_fill_reference, _add_reactions_from_model, _matlab_round_trip, _change_rules, + _change_rules_append, _add_from_equations, _remove_low_score_genes, +], ids=lambda f: f.__name__.lstrip("_")) +def test_operation_keeps_raven_gene_ids(operation, tmp_path): + out, rid = operation(_model(tmp_path), tmp_path) + assert sorted(g.id for g in out.reactions.get_by_id(rid).genes) == GENES + assert {"G1", "G2"}.isdisjoint(g.id for g in out.genes) + + +def test_expand_model_keeps_raven_gene_ids(tmp_path): + from raven_toolbox.manipulation import expand_model + m = _model(tmp_path) + m.reactions.R1.gpr = parse_gr_rule("(G1&G2 and G3) or G4") + expand_model(m) + expanded = [r for r in m.reactions if r.id.startswith("R1_EXP")] + assert sorted({g.id for r in expanded for g in r.genes}) == ["G1&G2", "G3", "G4"] + + +def test_set_gene_rule_copies_a_gpr_and_reads_strings_in_raven_grammar(): + from raven_toolbox.utils.gpr import set_gene_rule + + source = cobra.Reaction("S") + source.gpr = parse_gr_rule(RULE) + target = cobra.Reaction("T") + set_gene_rule(target, source.gpr) + assert sorted(g.id for g in target.genes) == GENES + assert target.gpr is not source.gpr + set_gene_rule(target, "") + assert target.gene_reaction_rule == "" and not target.genes + set_gene_rule(target, RULE) + assert sorted(g.id for g in target.genes) == GENES + + +@pytest.mark.filterwarnings("ignore:Malformed gene_reaction_rule") +def test_rule_to_gpr_falls_back_to_cobra_for_a_rule_raven_refuses(): + from cobra.core.gene import GPR + + from raven_toolbox.utils.gpr import rule_to_gpr + + rule = "G1 G2 or G3" + assert rule_to_gpr(rule).to_string() == GPR.from_string(rule).to_string() + assert rule_to_gpr(None).body is None From c2e1b7cd7b5596b7dae193871edec0923df0482e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:11:59 +0100 Subject: [PATCH 24/35] remove_genes: state that nothing changes when no requested gene is in the model --- CHANGELOG.md | 3 ++- src/raven_toolbox/manipulation/remove.py | 14 ++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b291b00..2ac696a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,8 @@ Milestones in the raven-toolbox port. For function-level status see 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. A gene kept only by a constrained reaction stays. + `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. diff --git a/src/raven_toolbox/manipulation/remove.py b/src/raven_toolbox/manipulation/remove.py index e1df9f2..074cc18 100644 --- a/src/raven_toolbox/manipulation/remove.py +++ b/src/raven_toolbox/manipulation/remove.py @@ -80,12 +80,14 @@ def remove_genes( * ``"remove"`` — delete them. * ``"keep"`` — leave them with an empty GPR and unchanged bounds. - Besides the requested genes, every gene that no remaining reaction named in its - GPR before the removal is deleted too, as RAVEN's ``removeGenes`` does: genes - that no reaction uses (orphan genes, also when not requested) and, with - ``"remove"``, genes used only by the deleted reactions. A gene whose only - reaction is blocked and kept (``"constrain"`` or ``"keep"``), such as the other - subunit of a complex that lost one, stays in ``model.genes``. + Requested genes that are not in the model are ignored. When none of the requested + genes is in the model, the model is not changed and ``[]`` is returned. Otherwise, + besides the requested genes, every gene that no remaining reaction named in its GPR + before the removal is deleted too, as RAVEN's ``removeGenes`` does: genes that no + reaction uses (orphan genes, also when not requested) and, with ``"remove"``, genes + used only by the deleted reactions. A gene whose only reaction is blocked and kept + (``"constrain"`` or ``"keep"``), such as the other subunit of a complex that lost + one, stays in ``model.genes``. ``remove_orphans`` (only meaningful with ``blocked_reactions="remove"``) passes through to cobra: also drop metabolites orphaned by the removal. From 1d5e55b40820dbf21966cddd8358771b8383cab8 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:12:00 +0100 Subject: [PATCH 25/35] benchmark_replicate_yeast_gem: drop_single_metabolite_reactions defaults to False --- scripts/benchmark_replicate_yeast_gem.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/benchmark_replicate_yeast_gem.py b/scripts/benchmark_replicate_yeast_gem.py index a72549e..091e0c4 100644 --- a/scripts/benchmark_replicate_yeast_gem.py +++ b/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.""" From 91ff9d062c2d168ac27337946bdd9039b16682e2 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:12:00 +0100 Subject: [PATCH 26/35] tests: ftinit ranks a merged reaction by the smallest original id of its group --- tests/test_init_pipeline.py | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 8d6b89d..a11264d 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -10,6 +10,7 @@ TEST_MODEL_FTINIT_WITH_TASK, TEST_MODEL_SCORES, TEST_MODEL_TASK_ESSENTIAL_MERGED, + _build, expr_for_rxn_score, make_test_model, make_test_task, @@ -251,6 +252,47 @@ def test_ftinit_resolve_ties_matches_oracle_and_is_stable(): assert {r.id for r in out1.reactions} == {r.id for r in out2.reactions} +def _two_merged_paths_model(): + """a[c] reaches m[c] by path P (Rm, Rc) or path Q (Rn, Rb); each path is one linear merge group. + + Model order and group numbers put P first ({R1, R2} is group 1, P group 2, Q group 3), and + the merged reactions are named by their first members Rm and Rn, but the smallest original + id of Q (Rb) sorts before that of P (Rc). + """ + mets = {"as": ("a", "s"), "ac": ("a", "c"), "xc": ("x", "c"), "yc": ("y", "c"), + "mc": ("m", "c"), "ec": ("e", "c"), "es": ("e", "s")} + rxns = { + "R1": ({"as": 1}, 0, 1000, ""), + "R2": ({"as": -1, "ac": 1}, -1000, 1000, ""), + "Rm": ({"ac": -1, "xc": 1}, 0, 1000, "Gm"), + "Rc": ({"xc": -1, "mc": 1}, 0, 1000, "Gc"), + "Rn": ({"ac": -1, "yc": 1}, 0, 1000, "Gn"), + "Rb": ({"yc": -1, "mc": 1}, 0, 1000, "Gb"), + "E": ({"mc": -1, "ec": 1}, 0, 1000, "GE"), + "R7": ({"ec": -1, "es": 1}, 0, 1000, "G7"), + "R8": ({"es": -1}, 0, 1000, ""), + } + return _build("twoPaths", mets, rxns, "R8") + + +def test_ftinit_resolve_ties_ranks_a_merged_reaction_by_its_smallest_original_id(): + """The tie between P and Q goes to Q, whose smallest member id Rb sorts first. + + Ranking by the merged reaction's id (Rm, Rn), by model order or by group number would keep P. + """ + model = _two_merged_paths_model() + prep = prep_init_model(model, [make_test_task()], ext_comp="s") + groups = prep.group_of + assert groups["Rm"] == groups["Rc"] < groups["Rn"] == groups["Rb"] + assert {"Rm", "Rn"} <= {r.id for r in prep.min_model.reactions} + scores = score_reactions_from_genes( + model, {"Gm": -1.0, "Gc": -1.0, "Gn": -1.0, "Gb": -1.0, "GE": 5.0, "G7": 5.0}) + out = ftinit(prep, rxn_scores=scores, resolve_ties=True) + kept = {r.id for r in out.reactions} + assert {"Rn", "Rb"} <= kept + assert kept.isdisjoint({"Rm", "Rc"}) + + # --------------------------------------------------------------------------- # # metabolomics (production-bonus for detected metabolites) — full pipeline. # Exercises what run_ftinit's own metabolomics tests cannot: metabolite *name* From d71652b6797aa8a0cad75ade3d25f2e541fbc880 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:25:01 +0100 Subject: [PATCH 27/35] tests: the paths of the rank-key test model are reversible, as in RAVEN's test --- tests/test_init_pipeline.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index a11264d..027ae54 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -257,17 +257,18 @@ def _two_merged_paths_model(): Model order and group numbers put P first ({R1, R2} is group 1, P group 2, Q group 3), and the merged reactions are named by their first members Rm and Rn, but the smallest original - id of Q (Rb) sorts before that of P (Rc). + id of Q (Rb) sorts before that of P (Rc). The model is that of RAVEN's tINIT test + ftINITResolveTiesRanksMergedRxnBySmallestOriginalId. """ mets = {"as": ("a", "s"), "ac": ("a", "c"), "xc": ("x", "c"), "yc": ("y", "c"), "mc": ("m", "c"), "ec": ("e", "c"), "es": ("e", "s")} rxns = { "R1": ({"as": 1}, 0, 1000, ""), "R2": ({"as": -1, "ac": 1}, -1000, 1000, ""), - "Rm": ({"ac": -1, "xc": 1}, 0, 1000, "Gm"), - "Rc": ({"xc": -1, "mc": 1}, 0, 1000, "Gc"), - "Rn": ({"ac": -1, "yc": 1}, 0, 1000, "Gn"), - "Rb": ({"yc": -1, "mc": 1}, 0, 1000, "Gb"), + "Rm": ({"ac": -1, "xc": 1}, -1000, 1000, "Gm"), + "Rc": ({"xc": -1, "mc": 1}, -1000, 1000, "Gc"), + "Rn": ({"ac": -1, "yc": 1}, -1000, 1000, "Gn"), + "Rb": ({"yc": -1, "mc": 1}, -1000, 1000, "Gb"), "E": ({"mc": -1, "ec": 1}, 0, 1000, "GE"), "R7": ({"ec": -1, "es": 1}, 0, 1000, "G7"), "R8": ({"es": -1}, 0, 1000, ""), From 481ee987f737e94086bba0a14c7fab6877cc564b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:34:36 +0100 Subject: [PATCH 28/35] run_ftinit: resolve_ties warns when the minimal-flux or the direction phase finds no solution --- CHANGELOG.md | 3 +++ src/raven_toolbox/init/ftinit.py | 36 ++++++++++++++++++++++++-------- 2 files changed, 30 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2ac696a..e2aded0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -85,6 +85,9 @@ Milestones in the raven-toolbox port. For function-level status see 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` warn (`phase2c-minflux`, + `phase2d-direction`) when the minimal-flux or the direction phase finds no solution, as + RAVEN does; the fluxes of the phase before it are kept (was: kept without a warning). ## 3.0.0b1 — 2026-09-07 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 6c3b497..21357cb 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -188,9 +188,12 @@ def run_ftinit( exhaust ``time_limit``; when that happens the (still-adopted) incumbent is an unproven tie-break and a warning is raised. The tie-break is skipped when the primary solve itself stops at ``time_limit``, since there is then no proven optimum to hold. - After it, the reported ``fluxes`` are a least-total-flux solution for the selected on + Phases 2a and 2b get half of ``time_limit`` each (rounded down, at least 1 s). + After them, the reported ``fluxes`` are a least-total-flux solution for the selected on set (:func:`_min_abs_flux`) and, among those, the one with the least id-rank-weighted - backward flux (:func:`_min_backward_flux`). See :func:`_resolve_ties`. + backward flux (:func:`_min_backward_flux`), each with the full ``time_limit``. When one + of these two flux phases finds no solution, the fluxes of the phase before it are kept + and a warning is raised. See :func:`_resolve_ties`. ``rank_keys`` maps a reaction id to the key the id-rank phases of ``resolve_ties`` sort it by (default: the id itself). :func:`ftinit` passes the smallest original @@ -668,9 +671,9 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): Returns ``(abs_expr, total, free_sign)`` if it produced a usable solution (the caller then reads its fluxes): the ``Σ|v|`` expression, its value and the ids of the - reactions whose flux can take either sign. Returns ``None`` otherwise (the caller keeps - the phase-2b fluxes). It never fails the extraction; a phase ending at the time limit - is adopted with a warning. + reactions whose flux can take either sign. Returns ``None`` otherwise, with a + ``phase2c-minflux`` warning (the caller keeps the phase-2b fluxes). It never fails the + extraction; a phase ending at the time limit is adopted with a warning. """ # Read every value before changing a bound: a bound change discards the solution. values = [(ind, ind.primal or 0.0) for ind, _ in indicators.values()] @@ -713,6 +716,12 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): ) if opt.status in ("optimal", "feasible", "suboptimal", "time_limit") and _has_solution(opt): return abs_expr, float(opt.objective.value), free_sign + warnings.warn( + "ftINIT minimal-flux phase (resolve_ties) found no solution (phase2c-minflux): " + "the fluxes, and hence the next step's forced directions, stay those of the " + "tie-broken on/off state.", + stacklevel=3, + ) return None @@ -732,7 +741,8 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k the full ``time_limit`` as its budget. Returns ``True`` if it produced a usable solution (the caller then reads its fluxes), - ``False`` otherwise (the caller keeps the phase-2c fluxes). It never fails the + ``False`` otherwise (the caller keeps the phase-2c fluxes), with a + ``phase2d-direction`` warning when the solve found nothing. It never fails the extraction; a phase ending at the time limit is adopted with a warning. """ if not free_sign: @@ -773,8 +783,15 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k "for a proven solution.", stacklevel=3, ) - return (opt.status in ("optimal", "feasible", "suboptimal", "time_limit") - and _has_solution(opt)) + found = (opt.status in ("optimal", "feasible", "suboptimal", "time_limit") + and _has_solution(opt)) + if not found: + warnings.warn( + "ftINIT direction phase (resolve_ties) found no solution (phase2d-direction): " + "the fluxes stay those of the minimal-flux phase.", + stacklevel=3, + ) + return found def _nudge_scores(rxn_scores: Mapping[str, float]) -> dict[str, float]: @@ -1011,7 +1028,8 @@ def ftinit( solver — applied both to each extraction step and to the task gap-fill. In each extraction step it then takes the fluxes from the least-total-flux solution for that on set, so the directions and forced magnitudes the next step reads from - them do not depend on the solver's vertex either. + them do not depend on the solver's vertex either. A step whose primary solve + stops at ``time_limit`` keeps its primary solution without these phases. * ``prove_abs_gap`` replaces the loose relative-gap escalation with a single solve per step proven to this fixed *absolute* gap, so the kept-set objective is optimal at any scale — not an arbitrary within-gap incumbent. **Try 1.0**: on Human-GEM/DLD1 From df789921f9f97d2f71532366e9035861da966340 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:34:36 +0100 Subject: [PATCH 29/35] tests: resolve_ties warns when phase 2c or 2d finds no solution --- tests/test_init_ftinit.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 4585e25..3dfa157 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -295,6 +295,43 @@ def solve(rank_keys=None): assert solve({"R1": "R9"}) == pytest.approx((-0.1, 0.0), abs=1e-6) +def _no_solution_in(monkeypatch, constraint_prefix): + """Make ``_has_solution`` report no solution while ``opt`` holds a constraint whose name + starts with ``constraint_prefix``, so the phase that adds it finds nothing.""" + import importlib + + ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") + real = ftinit_mod._has_solution + + def has_solution(opt): + if any(c.name.startswith(constraint_prefix) for c in opt.constraints): + return False + return real(opt) + + monkeypatch.setattr(ftinit_mod, "_has_solution", has_solution) + + +def test_resolve_ties_warns_when_the_minimal_flux_phase_finds_nothing(monkeypatch): + """When phase 2c finds no solution, a phase2c-minflux warning is raised and the + tie-broken on set is kept.""" + _no_solution_in(monkeypatch, "_absp_") + with pytest.warns(UserWarning, match="phase2c-minflux"): + res = run_ftinit(_direction_model(), {"R2": -1.0, "R3": 1.0}, essential_rxns=["E"], + resolve_ties=True) + assert "R3" in res.on_reactions and "R2" not in res.kept_reactions + assert abs(res.fluxes["E"]) >= 0.1 - 1e-6 + + +def test_resolve_ties_warns_when_the_direction_phase_finds_nothing(monkeypatch): + """When phase 2d finds no solution, a phase2d-direction warning is raised and the + fluxes stay a least-total-flux solution of phase 2c.""" + _no_solution_in(monkeypatch, "_minflux_cap") + with pytest.warns(UserWarning, match="phase2d-direction"): + res = run_ftinit(_two_route_model(), {"N": -1.0}, essential_rxns=["E"], + resolve_ties=True) + assert abs(res.fluxes["R1"]) + abs(res.fluxes["R2"]) == pytest.approx(0.1, abs=1e-6) + + # --------------------------------------------------------------------------- # # seed as an explicit parameter, and the unproven-incumbent warning. # --------------------------------------------------------------------------- # From b23b9bdf1f3471ef5b58e75ed4a35feddb6c351c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:43:24 +0100 Subject: [PATCH 30/35] ftinit: resolve_ties orients a merged reaction as its smallest member and skips the flux phases in the last step --- CHANGELOG.md | 7 ++++ src/raven_toolbox/init/ftinit.py | 55 +++++++++++++++++++++++++------- 2 files changed, 50 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e2aded0..711f341 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -88,6 +88,13 @@ Milestones in the raven-toolbox port. For function-level status see * `init.ftinit` and `init.run_ftinit` with `resolve_ties` warn (`phase2c-minflux`, `phase2d-direction`) when the minimal-flux or the direction phase finds no solution, as RAVEN does; the fluxes of the phase before it are kept (was: kept without a warning). +* `init.ftinit` with `resolve_ties`: the direction phase minimises the backward flux of a + merged reaction in the orientation of its member with the smallest original id, and a + turned-on reaction with zero flux is forced in that orientation, so neither depends on which + member names the merged reaction (new `run_ftinit` argument `orientation`). +* `init.ftinit` with `resolve_ties`: the last step skips the minimal-flux and direction phases, + as no step reads its fluxes (new `run_ftinit` argument `pin_fluxes`). The reaction set does + not change. ## 3.0.0b1 — 2026-09-07 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 21357cb..e53c36a 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -143,6 +143,8 @@ def run_ftinit( prove_abs_gap: float | None = None, resolve_ties: bool = False, rank_keys: Mapping[str, str] | None = None, + orientation: Mapping[str, int] | None = None, + pin_fluxes: bool = True, metabolomics: Mapping[str, Iterable[str]] | None = None, prod_weight: float = 0.5, seed: int = _EXTRACT_SEED, @@ -200,6 +202,17 @@ def run_ftinit( reaction id of each merge group, so the order does not depend on which member names the merged reaction. + ``orientation`` maps a reaction id to ``1`` or ``-1`` (default ``1``): the sign that + turns its flux into the direction of its reference orientation. The direction phase of + ``resolve_ties`` minimises the backward flux in that orientation. :func:`ftinit` passes + ``-1`` for a merged reaction in which the member with the smallest original id runs + reversed, so the preferred direction does not depend on which member names the merged + reaction. + + ``pin_fluxes`` (default true): if false, ``resolve_ties`` runs phases 2a and 2b only and + the reported ``fluxes`` are those of phase 2b. :func:`ftinit` passes false in its last + step, whose fluxes no later step reads. + ``metabolomics`` (already in **this model's** id space — :func:`ftinit` resolves a caller-facing metabolite-name list into this once) maps a detected-metabolite label to the set of reaction ids that produce it. A metabolite's ``mon`` indicator @@ -524,13 +537,15 @@ def rank_key(rid: str) -> str: primary_obj, time_limit, rank_key): on, fluxes = _read_solution() # Like phases 2a/2b, only when there is a removable (negative-score) reaction. + # Phases 2c and 2d set the fluxes only, which the next ftINIT step reads. held = (_min_abs_flux(opt, prob, indicators, flux_terms, time_limit) - if any(score < 0 for _, score in indicators.values()) else None) + if pin_fluxes and any(score < 0 for _, score in indicators.values()) + else None) if held is not None: fluxes = _read_solution()[1] abs_expr, total, free_sign = held if _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, - rank_key, time_limit): + rank_key, time_limit, orientation): fluxes = _read_solution()[1] kept = free_or_essential | on @@ -726,14 +741,15 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_key, - time_limit) -> bool: + time_limit, orientation=None) -> bool: """Phase 2d of ``resolve_ties``: fix the directions left open by the minimal total flux. Phase 2c can leave the direction of a flux free when two directions carry the same total flux, for example a pair of reversible reactions that close a loop in either direction. This phase keeps the on/off state of phase 2c, holds ``Σ|v|`` at most at its phase-2c - value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -v)`` over the reactions - whose flux can take either sign. The weight is the reaction's id rank from the end: + value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -o·v)`` over the reactions + whose flux can take either sign, with ``o`` the reaction's ``orientation`` (default 1). + The weight is the reaction's id rank from the end: sorted by ``rank_key``, the first of ``n`` reactions gets ``n`` and the last ``1``, so a reaction with a lower id prefers its forward direction. On Gurobi the phase-2c solution is its MIP start, as the tight cap on ``Σ|v|`` leaves the solver's heuristics little @@ -753,8 +769,9 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k n = len(free_sign) weight = {rid: float(n - i) for i, rid in enumerate(sorted(free_sign, key=rank_key))} extra, bwd_terms = [prob.Constraint(abs_expr, ub=total + tol, name="_minflux_cap")], [] + orient = {rid: (orientation or {}).get(rid, 1) for rid in free_sign} for rid in free_sign: - net = add([mul([Real(sign), var]) for var, sign in flux_terms[rid]]) + net = add([mul([Real(orient[rid] * sign), var]) for var, sign in flux_terms[rid]]) b = prob.Variable(f"_bwd_{rid}", lb=0.0) extra += [b, prob.Constraint(b + net, lb=0.0, name=f"_bwdc_{rid}")] bwd_terms.append(mul([Real(weight[rid]), b])) @@ -770,7 +787,7 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k grb.getVarByName(name).Start = val for rid in free_sign: net = sum(sign * start[var.name] for var, sign in flux_terms[rid]) - grb.getVarByName(f"_bwd_{rid}").Start = max(0.0, -net) + grb.getVarByName(f"_bwd_{rid}").Start = max(0.0, -orient[rid] * net) if time_limit is not None: opt.configuration.timeout = max(1, int(time_limit)) opt.optimize() @@ -872,7 +889,7 @@ def _metabolomics_producers(prep, metabolomics: Iterable[str]) -> dict[str, set[ def _solve_step( min_model, scores, step, *, essential, directions, ess_force, force_on, big_m, mip_gap, mip_gap_abs, time_limit, prove_abs_gap=None, resolve_ties=False, - rank_keys=None, metabolomics=None, prod_weight=0.5, + rank_keys=None, orientation=None, pin_fluxes=True, metabolomics=None, prod_weight=0.5, seed=_EXTRACT_SEED, threads=_EXTRACT_THREADS, feas_tol=_EXTRACT_TOL, opt_tol=_EXTRACT_TOL, int_feas_tol=_EXTRACT_TOL, presolve=_EXTRACT_PRESOLVE, @@ -894,6 +911,7 @@ def _run(mg, mga, tl, prove=None): force_on=force_on, force_on_ess=force_on, big_m=big_m, mip_gap=mg, mip_gap_abs=mga, time_limit=tl, prove_abs_gap=prove, resolve_ties=resolve_ties, rank_keys=rank_keys, + orientation=orientation, pin_fluxes=pin_fluxes, metabolomics=metabolomics, prod_weight=prod_weight, seed=seed, threads=threads, feas_tol=feas_tol, opt_tol=opt_tol, int_feas_tol=int_feas_tol, presolve=presolve, @@ -1028,8 +1046,11 @@ def ftinit( solver — applied both to each extraction step and to the task gap-fill. In each extraction step it then takes the fluxes from the least-total-flux solution for that on set, so the directions and forced magnitudes the next step reads from - them do not depend on the solver's vertex either. A step whose primary solve - stops at ``time_limit`` keeps its primary solution without these phases. + them do not depend on the solver's vertex either; among those, it prefers the + forward direction of the lower-ranked reaction, a merged reaction ranked and + oriented by the member with the smallest original id. The last step, whose + fluxes no step reads, skips these flux phases. A step whose primary solve stops + at ``time_limit`` keeps its primary solution without any of these phases. * ``prove_abs_gap`` replaces the loose relative-gap escalation with a single solve per step proven to this fixed *absolute* gap, so the kept-set objective is optimal at any scale — not an arbitrary within-gap incumbent. **Try 1.0**: on Human-GEM/DLD1 @@ -1095,13 +1116,17 @@ def ftinit( # the permanent (prep) essentials too, not only reactions turned on by a prior step. flux_of: dict[str, float] = {r.id: force_on for r in min_model.reactions} # resolve_ties ranks a merged reaction by its group's smallest original id, which does - # not depend on which member names the merged reaction. + # not depend on which member names the merged reaction. The reference orientation of a + # merged reaction is that of the same member: -1 when that member runs reversed in the + # merged reaction. smallest: dict[int, str] = {} for rid, gid in zip(prep.orig_rxn_ids, prep.group_ids, strict=True): if gid and (gid not in smallest or rid < smallest[gid]): smallest[gid] = rid rank_keys = {r.id: smallest[group_of[r.id]] for r in min_model.reactions if group_of.get(r.id, 0)} + is_reversed = dict(zip(prep.orig_rxn_ids, prep.reversed_rxns, strict=True)) + orientation = {rid: -1 for rid, key in rank_keys.items() if is_reversed[key]} n_steps = len(steps) for i, step in enumerate(steps): if verbose: @@ -1114,7 +1139,9 @@ def ftinit( if step.how_to_use_prev == "essential": for rid, flux in turned_on.items(): essential.add(rid) - directions[rid] = 1 if flux >= 0 else -1 + # A zero flux takes the reference orientation of the reaction. + directions[rid] = (1 if flux > 0 else -1 if flux < 0 + else orientation.get(rid, 1)) ess_force = {rid: min(abs(flux_of.get(rid, force_on)) * 0.99, force_on) for rid in essential} res = _solve_step( @@ -1122,6 +1149,10 @@ def ftinit( ess_force=ess_force, force_on=force_on, big_m=big_m, mip_gap=mip_gap, mip_gap_abs=mip_gap_abs, time_limit=time_limit, prove_abs_gap=prove_abs_gap, resolve_ties=resolve_ties, rank_keys=rank_keys, + orientation=orientation, + # The fluxes of a step set the forced directions of the next one, so the last + # step skips the flux phases of resolve_ties. + pin_fluxes=i < n_steps - 1, metabolomics=met_producers, prod_weight=prod_weight, seed=seed, threads=threads, feas_tol=feas_tol, opt_tol=opt_tol, int_feas_tol=int_feas_tol, presolve=presolve, From 83ee5699c1a364ccfaef61fed64d0aad378d4d85 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:43:24 +0100 Subject: [PATCH 31/35] tests: resolve_ties follows orientation, pin_fluxes skips the flux phases, ftinit passes both --- tests/test_init_ftinit.py | 24 +++++++++++++++++++ tests/test_init_pipeline.py | 46 +++++++++++++++++++++++++++++++++++++ 2 files changed, 70 insertions(+) diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 3dfa157..0927897 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -295,6 +295,30 @@ def solve(rank_keys=None): assert solve({"R1": "R9"}) == pytest.approx((-0.1, 0.0), abs=1e-6) +def test_resolve_ties_follows_the_orientation_of_a_reaction(): + """With orientation -1 for R1, its backward flux in that orientation is a forward flux, + so R1, which ranks first, runs backward instead of R2.""" + res = run_ftinit(_two_route_model(), {"N": -1.0}, essential_rxns=["E"], + resolve_ties=True, orientation={"R1": -1}) + assert (res.fluxes["R1"], res.fluxes["R2"]) == pytest.approx((-0.1, 0.0), abs=1e-6) + + +def test_resolve_ties_without_pin_fluxes_skips_the_flux_phases(monkeypatch): + """With pin_fluxes=False, phases 2c and 2d do not run and R1 is still the kept one.""" + import importlib + + ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") + calls = [] + monkeypatch.setattr(ftinit_mod, "_min_abs_flux", lambda *a, **k: calls.append(a)) + res = run_ftinit(_degenerate_model(), {"R1": -1.0, "R2": -1.0}, essential_rxns=["E"], + resolve_ties=True, pin_fluxes=False) + assert calls == [] + assert "R1" in res.kept_reactions and "R2" not in res.kept_reactions + run_ftinit(_degenerate_model(), {"R1": -1.0, "R2": -1.0}, essential_rxns=["E"], + resolve_ties=True) + assert len(calls) == 1 + + def _no_solution_in(monkeypatch, constraint_prefix): """Make ``_has_solution`` report no solution while ``opt`` holds a constraint whose name starts with ``constraint_prefix``, so the phase that adds it finds nothing.""" diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 027ae54..3a10734 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -294,6 +294,52 @@ def test_ftinit_resolve_ties_ranks_a_merged_reaction_by_its_smallest_original_id assert kept.isdisjoint({"Rm", "Rc"}) +def _flipped_member_model(): + """h[c] reaches g[c] through y[c] (Rz, Ra) or through y2[c] (Rp, Rq). + + Rz and Ra are both written as producers of y[c], so the merge runs Ra reversed in the + merged reaction Rz. The model is that of RAVEN's tINIT test + prepINITModelReturnsReversedRxnsOfTheMerge. + """ + mets = {"hs": ("h", "s"), "hc": ("h", "c"), "gc": ("g", "c"), "y1": ("y", "c"), + "y2": ("y2", "c"), "gs": ("g", "s")} + rxns = { + "EX_h": ({"hs": 1}, -1000, 1000, ""), + "Th": ({"hs": -1, "hc": 1}, -1000, 1000, ""), + "Rz": ({"hc": -1, "y1": 1}, -1000, 1000, "Gz"), + "Ra": ({"gc": -1, "y1": 1}, -1000, 1000, "Ga"), + "Rp": ({"hc": -1, "y2": 1}, -1000, 1000, "Gp"), + "Rq": ({"y2": -1, "gc": 1}, -1000, 1000, "Gq"), + "Tg": ({"gc": -1, "gs": 1}, -1000, 1000, ""), + "EX_g": ({"gs": -1}, -1000, 1000, ""), + } + return _build("flippedMember", mets, rxns, "EX_g") + + +def test_ftinit_orients_a_merged_reaction_as_its_smallest_member_and_pins_fluxes_until_the_last_step( + monkeypatch): + """ftinit passes orientation -1 for the merged Rz, in which its smallest member Ra runs + reversed, and skips the flux phases of resolve_ties in its last step only.""" + import importlib + + ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") + model = _flipped_member_model() + prep = prep_init_model(model, ext_comp="s") + assert [r for r, rev in zip(prep.orig_rxn_ids, prep.reversed_rxns, strict=True) if rev] == ["Ra"] + assert "Rz" in {r.id for r in prep.min_model.reactions} + seen = [] + real = ftinit_mod._solve_step + + def spy(*args, **kwargs): + seen.append((kwargs["orientation"], kwargs["pin_fluxes"])) + return real(*args, **kwargs) + + monkeypatch.setattr(ftinit_mod, "_solve_step", spy) + scores = score_reactions_from_genes(model, {"Gz": -1.0, "Ga": -1.0, "Gp": 1.0, "Gq": 1.0}) + ftinit(prep, rxn_scores=scores, resolve_ties=True, fill_gaps=False) + assert seen == [({"Rz": -1}, True), ({"Rz": -1}, False)] + + # --------------------------------------------------------------------------- # # metabolomics (production-bonus for detected metabolites) — full pipeline. # Exercises what run_ftinit's own metabolomics tests cannot: metabolite *name* From c812536b6ede090a7945023ceed7f0535313dc09 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:50:02 +0100 Subject: [PATCH 32/35] merge_compartments: docstring states where notes["merged_base_id"] survives and what happens when it is lost --- CHANGELOG.md | 5 +++++ src/raven_toolbox/manipulation/compartments.py | 11 ++++++++++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 711f341..18aed0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -30,6 +30,11 @@ Milestones in the raven-toolbox port. For function-level status see 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 diff --git a/src/raven_toolbox/manipulation/compartments.py b/src/raven_toolbox/manipulation/compartments.py index d67e078..a051c4c 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -86,7 +86,16 @@ def merge_compartments( metabolites would get the same base id (a ``base_metabolite`` that keeps ``X_c`` and ``X_m`` apart), both store their own id instead, so they stay distinct species. The note survives ``model.copy()``, cobra's JSON and SBML files and :func:`~raven_toolbox.io.write_yaml_model` - with :func:`~raven_toolbox.io.read_yaml_model`. + with :func:`~raven_toolbox.io.read_yaml_model`. RAVEN's ``readYAMLmodel`` shows it in + ``metNotes``, for example ``merged_base_id: ACCOA note: free text note`` for a metabolite + that also has the free-text note ``free text note``; no RAVEN function uses it. A round trip + through RAVEN YAML (``readYAMLmodel`` and ``writeYAMLmodel``) or through a .mat file loses + it: RAVEN joins the key and the free text into one string, which + :func:`~raven_toolbox.io.read_yaml_model` reads back as the free-text note, and + :func:`~raven_toolbox.io.write_matlab_model` and :func:`~raven_toolbox.io.read_matlab_model` + carry no metabolite notes. The localization functions then key the metabolite by its id + without the suffix of its own compartment, so a merged ``ACCOA_c`` in ``s`` gets the key + ``ACCOA_c``, and no error is raised. """ out = model.copy() key = base_metabolite if base_metabolite is not None else _base_id From 7150361f40070a4ebed149c4b83c6b5dc478cbf5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 10:22:32 +0100 Subject: [PATCH 33/35] run_ftinit: docstring and comments of pin_fluxes name the phases it skips --- src/raven_toolbox/init/ftinit.py | 6 +++--- tests/test_init_pipeline.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index e53c36a..fb97626 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -209,8 +209,8 @@ def run_ftinit( reversed, so the preferred direction does not depend on which member names the merged reaction. - ``pin_fluxes`` (default true): if false, ``resolve_ties`` runs phases 2a and 2b only and - the reported ``fluxes`` are those of phase 2b. :func:`ftinit` passes false in its last + ``pin_fluxes`` (default true): if false, ``resolve_ties`` skips phases 2c and 2d and the + reported ``fluxes`` are those of phase 2b. :func:`ftinit` passes false in its last step, whose fluxes no later step reads. ``metabolomics`` (already in **this model's** id space — :func:`ftinit` resolves a @@ -537,7 +537,7 @@ def rank_key(rid: str) -> str: primary_obj, time_limit, rank_key): on, fluxes = _read_solution() # Like phases 2a/2b, only when there is a removable (negative-score) reaction. - # Phases 2c and 2d set the fluxes only, which the next ftINIT step reads. + # Phases 2c and 2d change the fluxes, not the on set; the next ftINIT step reads them. held = (_min_abs_flux(opt, prob, indicators, flux_terms, time_limit) if pin_fluxes and any(score < 0 for _, score in indicators.values()) else None) diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 3a10734..49fd009 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -319,7 +319,7 @@ def _flipped_member_model(): def test_ftinit_orients_a_merged_reaction_as_its_smallest_member_and_pins_fluxes_until_the_last_step( monkeypatch): """ftinit passes orientation -1 for the merged Rz, in which its smallest member Ra runs - reversed, and skips the flux phases of resolve_ties in its last step only.""" + reversed, and skips the flux phases of resolve_ties in its last step and in no other.""" import importlib ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") From e3dcb559a6a83818f3a0ada8e16b76001e598bfd Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 17:40:34 +0100 Subject: [PATCH 34/35] ftinit: resolve_ties keeps the orientation of each merged reaction as written _min_backward_flux weights the backward flux by the rank key only and a turned-on reaction with zero flux is forced forward in the next step. run_ftinit takes the private keyword _pin_fluxes, which ftinit sets false in its last step. The tests check that run_ftinit skips phases 2c and 2d with _pin_fluxes false and that ftinit sets it false in the last step only. --- CHANGELOG.md | 8 +--- src/raven_toolbox/init/ftinit.py | 55 ++++++++++----------------- tests/test_init_ftinit.py | 12 +----- tests/test_init_pipeline.py | 65 ++++++++++++-------------------- 4 files changed, 48 insertions(+), 92 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 18aed0a..6f66c3b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -93,13 +93,9 @@ Milestones in the raven-toolbox port. For function-level status see * `init.ftinit` and `init.run_ftinit` with `resolve_ties` warn (`phase2c-minflux`, `phase2d-direction`) when the minimal-flux or the direction phase finds no solution, as RAVEN does; the fluxes of the phase before it are kept (was: kept without a warning). -* `init.ftinit` with `resolve_ties`: the direction phase minimises the backward flux of a - merged reaction in the orientation of its member with the smallest original id, and a - turned-on reaction with zero flux is forced in that orientation, so neither depends on which - member names the merged reaction (new `run_ftinit` argument `orientation`). * `init.ftinit` with `resolve_ties`: the last step skips the minimal-flux and direction phases, - as no step reads its fluxes (new `run_ftinit` argument `pin_fluxes`). The reaction set does - not change. + 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 diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index fb97626..2a725cf 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -143,8 +143,7 @@ def run_ftinit( prove_abs_gap: float | None = None, resolve_ties: bool = False, rank_keys: Mapping[str, str] | None = None, - orientation: Mapping[str, int] | None = None, - pin_fluxes: bool = True, + _pin_fluxes: bool = True, metabolomics: Mapping[str, Iterable[str]] | None = None, prod_weight: float = 0.5, seed: int = _EXTRACT_SEED, @@ -202,16 +201,9 @@ def run_ftinit( reaction id of each merge group, so the order does not depend on which member names the merged reaction. - ``orientation`` maps a reaction id to ``1`` or ``-1`` (default ``1``): the sign that - turns its flux into the direction of its reference orientation. The direction phase of - ``resolve_ties`` minimises the backward flux in that orientation. :func:`ftinit` passes - ``-1`` for a merged reaction in which the member with the smallest original id runs - reversed, so the preferred direction does not depend on which member names the merged - reaction. - - ``pin_fluxes`` (default true): if false, ``resolve_ties`` skips phases 2c and 2d and the - reported ``fluxes`` are those of phase 2b. :func:`ftinit` passes false in its last - step, whose fluxes no later step reads. + ``_pin_fluxes`` (default true) is set by :func:`ftinit`, not by callers: false in its + last step, whose fluxes no later step reads. ``resolve_ties`` then skips phases 2c and + 2d and the reported ``fluxes`` are those of phase 2b. ``metabolomics`` (already in **this model's** id space — :func:`ftinit` resolves a caller-facing metabolite-name list into this once) maps a detected-metabolite label @@ -539,13 +531,13 @@ def rank_key(rid: str) -> str: # Like phases 2a/2b, only when there is a removable (negative-score) reaction. # Phases 2c and 2d change the fluxes, not the on set; the next ftINIT step reads them. held = (_min_abs_flux(opt, prob, indicators, flux_terms, time_limit) - if pin_fluxes and any(score < 0 for _, score in indicators.values()) + if _pin_fluxes and any(score < 0 for _, score in indicators.values()) else None) if held is not None: fluxes = _read_solution()[1] abs_expr, total, free_sign = held if _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, - rank_key, time_limit, orientation): + rank_key, time_limit): fluxes = _read_solution()[1] kept = free_or_essential | on @@ -741,15 +733,14 @@ def _min_abs_flux(opt, prob, indicators, flux_terms, time_limit): def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_key, - time_limit, orientation=None) -> bool: + time_limit) -> bool: """Phase 2d of ``resolve_ties``: fix the directions left open by the minimal total flux. Phase 2c can leave the direction of a flux free when two directions carry the same total flux, for example a pair of reversible reactions that close a loop in either direction. This phase keeps the on/off state of phase 2c, holds ``Σ|v|`` at most at its phase-2c - value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -o·v)`` over the reactions - whose flux can take either sign, with ``o`` the reaction's ``orientation`` (default 1). - The weight is the reaction's id rank from the end: + value plus ``max(1e-7·Σ|v|, 1e-6)``, and minimises ``Σ w·max(0, -v)`` over the reactions + whose flux can take either sign. The weight is the reaction's id rank from the end: sorted by ``rank_key``, the first of ``n`` reactions gets ``n`` and the last ``1``, so a reaction with a lower id prefers its forward direction. On Gurobi the phase-2c solution is its MIP start, as the tight cap on ``Σ|v|`` leaves the solver's heuristics little @@ -769,9 +760,8 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k n = len(free_sign) weight = {rid: float(n - i) for i, rid in enumerate(sorted(free_sign, key=rank_key))} extra, bwd_terms = [prob.Constraint(abs_expr, ub=total + tol, name="_minflux_cap")], [] - orient = {rid: (orientation or {}).get(rid, 1) for rid in free_sign} for rid in free_sign: - net = add([mul([Real(orient[rid] * sign), var]) for var, sign in flux_terms[rid]]) + net = add([mul([Real(sign), var]) for var, sign in flux_terms[rid]]) b = prob.Variable(f"_bwd_{rid}", lb=0.0) extra += [b, prob.Constraint(b + net, lb=0.0, name=f"_bwdc_{rid}")] bwd_terms.append(mul([Real(weight[rid]), b])) @@ -787,7 +777,7 @@ def _min_backward_flux(opt, prob, flux_terms, abs_expr, total, free_sign, rank_k grb.getVarByName(name).Start = val for rid in free_sign: net = sum(sign * start[var.name] for var, sign in flux_terms[rid]) - grb.getVarByName(f"_bwd_{rid}").Start = max(0.0, -orient[rid] * net) + grb.getVarByName(f"_bwd_{rid}").Start = max(0.0, -net) if time_limit is not None: opt.configuration.timeout = max(1, int(time_limit)) opt.optimize() @@ -889,7 +879,7 @@ def _metabolomics_producers(prep, metabolomics: Iterable[str]) -> dict[str, set[ def _solve_step( min_model, scores, step, *, essential, directions, ess_force, force_on, big_m, mip_gap, mip_gap_abs, time_limit, prove_abs_gap=None, resolve_ties=False, - rank_keys=None, orientation=None, pin_fluxes=True, metabolomics=None, prod_weight=0.5, + rank_keys=None, pin_fluxes=True, metabolomics=None, prod_weight=0.5, seed=_EXTRACT_SEED, threads=_EXTRACT_THREADS, feas_tol=_EXTRACT_TOL, opt_tol=_EXTRACT_TOL, int_feas_tol=_EXTRACT_TOL, presolve=_EXTRACT_PRESOLVE, @@ -911,7 +901,7 @@ def _run(mg, mga, tl, prove=None): force_on=force_on, force_on_ess=force_on, big_m=big_m, mip_gap=mg, mip_gap_abs=mga, time_limit=tl, prove_abs_gap=prove, resolve_ties=resolve_ties, rank_keys=rank_keys, - orientation=orientation, pin_fluxes=pin_fluxes, + _pin_fluxes=pin_fluxes, metabolomics=metabolomics, prod_weight=prod_weight, seed=seed, threads=threads, feas_tol=feas_tol, opt_tol=opt_tol, int_feas_tol=int_feas_tol, presolve=presolve, @@ -1047,10 +1037,10 @@ def ftinit( extraction step it then takes the fluxes from the least-total-flux solution for that on set, so the directions and forced magnitudes the next step reads from them do not depend on the solver's vertex either; among those, it prefers the - forward direction of the lower-ranked reaction, a merged reaction ranked and - oriented by the member with the smallest original id. The last step, whose - fluxes no step reads, skips these flux phases. A step whose primary solve stops - at ``time_limit`` keeps its primary solution without any of these phases. + forward direction of the lower-ranked reaction, a merged reaction ranked by the + member with the smallest original id. The last step, whose fluxes no step reads, + skips these flux phases. A step whose primary solve stops at ``time_limit`` keeps + its primary solution without any of these phases. * ``prove_abs_gap`` replaces the loose relative-gap escalation with a single solve per step proven to this fixed *absolute* gap, so the kept-set objective is optimal at any scale — not an arbitrary within-gap incumbent. **Try 1.0**: on Human-GEM/DLD1 @@ -1116,17 +1106,13 @@ def ftinit( # the permanent (prep) essentials too, not only reactions turned on by a prior step. flux_of: dict[str, float] = {r.id: force_on for r in min_model.reactions} # resolve_ties ranks a merged reaction by its group's smallest original id, which does - # not depend on which member names the merged reaction. The reference orientation of a - # merged reaction is that of the same member: -1 when that member runs reversed in the - # merged reaction. + # not depend on which member names the merged reaction. smallest: dict[int, str] = {} for rid, gid in zip(prep.orig_rxn_ids, prep.group_ids, strict=True): if gid and (gid not in smallest or rid < smallest[gid]): smallest[gid] = rid rank_keys = {r.id: smallest[group_of[r.id]] for r in min_model.reactions if group_of.get(r.id, 0)} - is_reversed = dict(zip(prep.orig_rxn_ids, prep.reversed_rxns, strict=True)) - orientation = {rid: -1 for rid, key in rank_keys.items() if is_reversed[key]} n_steps = len(steps) for i, step in enumerate(steps): if verbose: @@ -1139,9 +1125,7 @@ def ftinit( if step.how_to_use_prev == "essential": for rid, flux in turned_on.items(): essential.add(rid) - # A zero flux takes the reference orientation of the reaction. - directions[rid] = (1 if flux > 0 else -1 if flux < 0 - else orientation.get(rid, 1)) + directions[rid] = 1 if flux >= 0 else -1 ess_force = {rid: min(abs(flux_of.get(rid, force_on)) * 0.99, force_on) for rid in essential} res = _solve_step( @@ -1149,7 +1133,6 @@ def ftinit( ess_force=ess_force, force_on=force_on, big_m=big_m, mip_gap=mip_gap, mip_gap_abs=mip_gap_abs, time_limit=time_limit, prove_abs_gap=prove_abs_gap, resolve_ties=resolve_ties, rank_keys=rank_keys, - orientation=orientation, # The fluxes of a step set the forced directions of the next one, so the last # step skips the flux phases of resolve_ties. pin_fluxes=i < n_steps - 1, diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 0927897..89f80da 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -295,23 +295,15 @@ def solve(rank_keys=None): assert solve({"R1": "R9"}) == pytest.approx((-0.1, 0.0), abs=1e-6) -def test_resolve_ties_follows_the_orientation_of_a_reaction(): - """With orientation -1 for R1, its backward flux in that orientation is a forward flux, - so R1, which ranks first, runs backward instead of R2.""" - res = run_ftinit(_two_route_model(), {"N": -1.0}, essential_rxns=["E"], - resolve_ties=True, orientation={"R1": -1}) - assert (res.fluxes["R1"], res.fluxes["R2"]) == pytest.approx((-0.1, 0.0), abs=1e-6) - - def test_resolve_ties_without_pin_fluxes_skips_the_flux_phases(monkeypatch): - """With pin_fluxes=False, phases 2c and 2d do not run and R1 is still the kept one.""" + """With _pin_fluxes=False, phases 2c and 2d do not run and R1 is still the kept one.""" import importlib ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") calls = [] monkeypatch.setattr(ftinit_mod, "_min_abs_flux", lambda *a, **k: calls.append(a)) res = run_ftinit(_degenerate_model(), {"R1": -1.0, "R2": -1.0}, essential_rxns=["E"], - resolve_ties=True, pin_fluxes=False) + resolve_ties=True, _pin_fluxes=False) assert calls == [] assert "R1" in res.kept_reactions and "R2" not in res.kept_reactions run_ftinit(_degenerate_model(), {"R1": -1.0, "R2": -1.0}, essential_rxns=["E"], diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 49fd009..ab736e9 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -294,50 +294,35 @@ def test_ftinit_resolve_ties_ranks_a_merged_reaction_by_its_smallest_original_id assert kept.isdisjoint({"Rm", "Rc"}) -def _flipped_member_model(): - """h[c] reaches g[c] through y[c] (Rz, Ra) or through y2[c] (Rp, Rq). +def test_ftinit_pins_fluxes_in_every_step_but_the_last(monkeypatch): + """ftinit runs the flux phases of resolve_ties in each step whose fluxes a later step + reads and skips them in the last step, as RAVEN's tINIT test + ftINITPinsFluxesInEveryStepButTheLast checks.""" + import importlib - Rz and Ra are both written as producers of y[c], so the merge runs Ra reversed in the - merged reaction Rz. The model is that of RAVEN's tINIT test - prepINITModelReturnsReversedRxnsOfTheMerge. - """ - mets = {"hs": ("h", "s"), "hc": ("h", "c"), "gc": ("g", "c"), "y1": ("y", "c"), - "y2": ("y2", "c"), "gs": ("g", "s")} - rxns = { - "EX_h": ({"hs": 1}, -1000, 1000, ""), - "Th": ({"hs": -1, "hc": 1}, -1000, 1000, ""), - "Rz": ({"hc": -1, "y1": 1}, -1000, 1000, "Gz"), - "Ra": ({"gc": -1, "y1": 1}, -1000, 1000, "Ga"), - "Rp": ({"hc": -1, "y2": 1}, -1000, 1000, "Gp"), - "Rq": ({"y2": -1, "gc": 1}, -1000, 1000, "Gq"), - "Tg": ({"gc": -1, "gs": 1}, -1000, 1000, ""), - "EX_g": ({"gs": -1}, -1000, 1000, ""), - } - return _build("flippedMember", mets, rxns, "EX_g") + ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") + model = _two_merged_paths_model() + prep = prep_init_model(model, [make_test_task()], ext_comp="s") + step, seen = [-1], [] + real_step, real_run = ftinit_mod._solve_step, ftinit_mod.run_ftinit + def step_spy(*args, **kwargs): + step[0] += 1 + return real_step(*args, **kwargs) -def test_ftinit_orients_a_merged_reaction_as_its_smallest_member_and_pins_fluxes_until_the_last_step( - monkeypatch): - """ftinit passes orientation -1 for the merged Rz, in which its smallest member Ra runs - reversed, and skips the flux phases of resolve_ties in its last step and in no other.""" - import importlib + def run_spy(*args, **kwargs): + seen.append((step[0], kwargs["_pin_fluxes"])) + return real_run(*args, **kwargs) - ftinit_mod = importlib.import_module("raven_toolbox.init.ftinit") - model = _flipped_member_model() - prep = prep_init_model(model, ext_comp="s") - assert [r for r, rev in zip(prep.orig_rxn_ids, prep.reversed_rxns, strict=True) if rev] == ["Ra"] - assert "Rz" in {r.id for r in prep.min_model.reactions} - seen = [] - real = ftinit_mod._solve_step - - def spy(*args, **kwargs): - seen.append((kwargs["orientation"], kwargs["pin_fluxes"])) - return real(*args, **kwargs) - - monkeypatch.setattr(ftinit_mod, "_solve_step", spy) - scores = score_reactions_from_genes(model, {"Gz": -1.0, "Ga": -1.0, "Gp": 1.0, "Gq": 1.0}) - ftinit(prep, rxn_scores=scores, resolve_ties=True, fill_gaps=False) - assert seen == [({"Rz": -1}, True), ({"Rz": -1}, False)] + monkeypatch.setattr(ftinit_mod, "_solve_step", step_spy) + monkeypatch.setattr(ftinit_mod, "run_ftinit", run_spy) + scores = score_reactions_from_genes( + model, {"Gm": -1.0, "Gc": -1.0, "Gn": -1.0, "Gb": -1.0, "GE": 5.0, "G7": 5.0}) + out = ftinit(prep, rxn_scores=scores, resolve_ties=True, fill_gaps=False) + n_steps = len(get_init_steps()) + assert {s for s, _ in seen} == set(range(n_steps)) + assert all(pin == (s < n_steps - 1) for s, pin in seen) + assert {"Rn", "Rb"} <= {r.id for r in out.reactions} # --------------------------------------------------------------------------- # From ce9b2d32539c62d737cd4068d0a525dfa602d613 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 21:06:02 +0100 Subject: [PATCH 35/35] CHANGELOG: warnings of the resolve_ties minimal-flux and direction phases --- CHANGELOG.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6f66c3b..1185539 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -90,9 +90,11 @@ Milestones in the raven-toolbox port. For function-level status see 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` warn (`phase2c-minflux`, - `phase2d-direction`) when the minimal-flux or the direction phase finds no solution, as - RAVEN does; the fluxes of the phase before it are kept (was: kept without a warning). +* `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.