From 35569338fb71c8174139ac9edeb3c9fe1673802e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 10:40:16 +0200 Subject: [PATCH 01/71] Run FVA and blocked-reaction searches in process in the gap-filling code cobra's flux_variability_analysis and find_blocked_reactions start a worker pool by default. With Python 3.14 (forkserver start method) the pool did not return in gap_report on a cluster node: the main thread waited in Pool.join. The gap-filling modules and the localization certificate now pass processes=1, as the certificate's FVA already did. --- src/raven_toolbox/gapfilling/fast_lp.py | 4 ++-- src/raven_toolbox/gapfilling/fill.py | 6 ++++-- src/raven_toolbox/gapfilling/report.py | 6 +++++- src/raven_toolbox/localization/certify.py | 5 +++-- tests/test_gapfilling_report.py | 19 +++++++++++++++++++ 5 files changed, 33 insertions(+), 7 deletions(-) diff --git a/src/raven_toolbox/gapfilling/fast_lp.py b/src/raven_toolbox/gapfilling/fast_lp.py index 58822196..ed23f7df 100644 --- a/src/raven_toolbox/gapfilling/fast_lp.py +++ b/src/raven_toolbox/gapfilling/fast_lp.py @@ -267,7 +267,7 @@ def fill_gaps_fast_lp( templates = _as_models(templates) # ---- Find blocked reactions in draft ---- - blocked = set(find_blocked_reactions(model)) + blocked = set(find_blocked_reactions(model, processes=1)) if not blocked: if verbose: print("fill_gaps_fast_lp: no blocked reactions found.") @@ -287,7 +287,7 @@ def fill_gaps_fast_lp( blocked_in_model = [r for r in blocked if r in working.reactions] if blocked_in_model: fva = flux_variability_analysis( - working, reaction_list=blocked_in_model, fraction_of_optimum=0.0 + working, reaction_list=blocked_in_model, fraction_of_optimum=0.0, processes=1 ) rescuable = [ r diff --git a/src/raven_toolbox/gapfilling/fill.py b/src/raven_toolbox/gapfilling/fill.py index bdb849d5..cc5bc795 100644 --- a/src/raven_toolbox/gapfilling/fill.py +++ b/src/raven_toolbox/gapfilling/fill.py @@ -160,14 +160,16 @@ def connect_blocked_reactions( if penalty <= 0: raise ValueError(f"penalty must be > 0 (a non-positive cost malforms the MILP); got {penalty}.") templates = _as_models(templates) - blocked = set(find_blocked_reactions(model)) + blocked = set(find_blocked_reactions(model, processes=1)) candidates = [r for r in blocked if model.reactions.get_by_id(r).lower_bound >= 0] working, template_ids = _merge_templates(model, templates) target: list[str] = [] if candidates: - fva = flux_variability_analysis(working, reaction_list=candidates, fraction_of_optimum=0.0) + fva = flux_variability_analysis( + working, reaction_list=candidates, fraction_of_optimum=0.0, processes=1 + ) # A reaction can be missing from the FVA frame if the solver dropped it # (e.g. the reaction was eliminated upstream); treat that as "unreachable" # rather than letting the KeyError propagate. diff --git a/src/raven_toolbox/gapfilling/report.py b/src/raven_toolbox/gapfilling/report.py index 5b2b348a..c8c50885 100644 --- a/src/raven_toolbox/gapfilling/report.py +++ b/src/raven_toolbox/gapfilling/report.py @@ -103,7 +103,11 @@ class GapReportResult: def _blocked_ids(model: cobra.Model, reactions: list[cobra.Reaction] | None = None) -> set[str]: - return set(find_blocked_reactions(model, reaction_list=reactions, zero_cutoff=_FLUX_THRESHOLD)) + return set( + find_blocked_reactions( + model, reaction_list=reactions, zero_cutoff=_FLUX_THRESHOLD, processes=1 + ) + ) def _relaxed_copy(model: cobra.Model) -> cobra.Model: diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index 1074f6a0..6e7f6421 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -420,7 +420,8 @@ def _usable_transports(applied, proposal, biomass_id, min_growth, conditions, de with applied: _apply_medium(applied, medium) applied.objective = biomass_id - unblocked = set(applied.reactions.list_attr("id")) - set(find_blocked_reactions(applied)) + blocked_now = set(find_blocked_reactions(applied, processes=1)) + unblocked = set(applied.reactions.list_attr("id")) - blocked_now for trid, bc in tr_of.items(): if trid in unblocked: keep.add(bc) @@ -826,7 +827,7 @@ def _diagnose_growth_gap(applied, proposal, biomass_reaction, movable_ids, pinne return None bio = applied.reactions.get_by_id(biomass_reaction) bio_comp = _rxn_compartment(applied.reactions.get_by_id(biomass_reaction)) - blocked = set(find_blocked_reactions(applied)) + blocked = set(find_blocked_reactions(applied, processes=1)) precursors = {base(m) for m, coeff in bio.metabolites.items() if coeff < 0} for rid in proposal.placements: if rid not in movable_ids or rid in proposal.added_reactions: diff --git a/tests/test_gapfilling_report.py b/tests/test_gapfilling_report.py index c727c705..8e11d94a 100644 --- a/tests/test_gapfilling_report.py +++ b/tests/test_gapfilling_report.py @@ -163,3 +163,22 @@ def test_forced_flux_infeasible_once_closed_raises(gap_report_model): gap_report_model.reactions.get_by_id("R1").lower_bound = 1.0 with pytest.raises(ValueError, match="once its boundary reactions are closed.*R1"): gap_report(gap_report_model) + + +def test_blocked_reactions_are_found_in_process(gap_report_model, monkeypatch): + """cobra's default worker pool can hang (seen with Python 3.14 on a cluster). + + The report runs FVA in-process. + """ + import raven_toolbox.gapfilling.report as report + + seen = [] + original = report.find_blocked_reactions + + def recording(model, **kwargs): + seen.append(kwargs.get("processes")) + return original(model, **kwargs) + + monkeypatch.setattr(report, "find_blocked_reactions", recording) + gap_report(gap_report_model) + assert seen and set(seen) == {1} From 212bc7625e0a0a4edea6f8bc2d1791d5a3d2bd27 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 10:54:29 +0200 Subject: [PATCH 02/71] Align exportForGit with RAVEN: all five formats, exact txt, fuller dependencies The default formats are now yml, xml, mat, xlsx and txt, as in RAVEN. The txt table is headed 'Rxn id', writes bounds and objective coefficients with the digits needed to read them back (RAVEN's exportForGit does the same) and pads no field. dependencies.txt lists python, cobra, libSBML, ruamel.yaml, scipy, openpyxl and raven_toolbox in a fixed order, then the model's version and date; it holds no time of export. --- src/raven_toolbox/io/git.py | 65 +++++++++++++++++++++++++++++++------ tests/test_io_git.py | 40 ++++++++++++++++++++++- 2 files changed, 94 insertions(+), 11 deletions(-) diff --git a/src/raven_toolbox/io/git.py b/src/raven_toolbox/io/git.py index 15c83fbc..1f4f1919 100644 --- a/src/raven_toolbox/io/git.py +++ b/src/raven_toolbox/io/git.py @@ -13,6 +13,7 @@ import datetime import importlib.metadata as _md +import math import platform from collections.abc import Iterable from pathlib import Path @@ -145,23 +146,70 @@ def _write_sbml(model: cobra.Model, path: Path) -> None: libsbml.writeSBMLToFile(doc, str(path)) +def _shortest(value: float) -> str: + """The shortest text that reads back as exactly ``value`` (1/3 is 0.3333333333333333). + + The same rule as RAVEN's exportForGit, so a bound is never rounded in the file. + """ + if math.isinf(value): + return "Inf" if value > 0 else "-Inf" + if math.isnan(value): + return "NaN" + if value == math.floor(value) and abs(value) < 1e15: + return str(int(value)) + for digits in range(1, 18): + text = f"{value:.{digits}g}" + if float(text) == value: + return text + return repr(value) + + def _write_txt(model: cobra.Model, path: Path) -> None: - """Single-file, human-readable reaction table (RAVEN exportForGit txt).""" + """Single-file, human-readable reaction table (RAVEN exportForGit txt). + + Bounds and objective coefficients are written with the digits needed to read them back + unchanged, and no field is padded. + """ with open(path, "w", encoding="utf-8") as fh: - fh.write("Rxn name\tFormula\tGene-reaction association\tLB\tUB\tObjective\n") + fh.write("Rxn id\tFormula\tGene-reaction association\tLB\tUB\tObjective\n") for r in model.reactions: fh.write( - f"{r.id}\t{_equation(r)}\t{r.gene_reaction_rule}\t" - f"{r.lower_bound:g}\t{r.upper_bound:g}\t{r.objective_coefficient:g}\n" + f"{r.id}\t{_equation(r).strip()}\t{r.gene_reaction_rule}\t" + f"{_shortest(r.lower_bound)}\t{_shortest(r.upper_bound)}\t" + f"{_shortest(r.objective_coefficient)}\n" ) +def _write_dependencies(path: Path, model: cobra.Model) -> None: + """Versions of everything that writes the files, then the model's own version and date. + + The order is fixed and nothing depends on the time of the export, so the same model + exported twice gives the same file. + """ + lines = [ + ("python", platform.python_version()), + ("cobra", _version("cobra")), + ("libSBML", libsbml.getLibSBMLDottedVersion()), + ("ruamel.yaml", _version("ruamel.yaml")), + ("scipy", _version("scipy")), + ("openpyxl", _version("openpyxl")), + ("raven_toolbox", _version("raven_toolbox")), + ] + meta = get_model_metadata(model) + for key in ("version", "date"): + if meta.get(key) is not None: + lines.append((f"model_{key}", str(meta[key]))) + with open(path, "w", encoding="utf-8") as fh: + for name, value in lines: + fh.write(f"{name}\t{value}\n") + + def export_for_git( model: cobra.Model, path: str | Path = ".", *, prefix: str = "model", - formats: Iterable[str] = ("yml", "xml", "mat", "xlsx"), + formats: Iterable[str] = _ALL_FORMATS, sub_dirs: bool = True, varname: str | None = None, version: str | None = None, @@ -177,7 +225,7 @@ def export_for_git( Base filename for every format (default ``"model"``). formats Which formats to write; any of ``"yml"``, ``"xml"``, ``"mat"``, - ``"xlsx"``, ``"txt"`` (default ``yml``/``xml``/``mat``/``xlsx``). + ``"xlsx"``, ``"txt"`` (default: all five, as RAVEN's ``exportForGit``). sub_dirs If True (default), write ``model//.`` (standard-GEM layout); otherwise all files go directly in ``path``. @@ -224,9 +272,6 @@ def target(fmt: str) -> Path: if "txt" in formats: _write_txt(model, target("txt")) - with open(root / "dependencies.txt", "w", encoding="utf-8") as fh: - fh.write(f"python\t{platform.python_version()}\n") - fh.write(f"cobra\t{_version('cobra')}\n") - fh.write(f"raven_toolbox\t{_version('raven_toolbox')}\n") + _write_dependencies(root / "dependencies.txt", model) return root diff --git a/tests/test_io_git.py b/tests/test_io_git.py index 8e59cdb6..41602102 100644 --- a/tests/test_io_git.py +++ b/tests/test_io_git.py @@ -62,11 +62,49 @@ def test_does_not_mutate_model(model, tmp_path): def test_txt_table_content(model, tmp_path): root = export_for_git(model, tmp_path, formats=("txt",)) txt = (root / "txt" / "model.txt").read_text() - assert txt.splitlines()[0].startswith("Rxn name\t") + assert txt.splitlines()[0] == "Rxn id\tFormula\tGene-reaction association\tLB\tUB\tObjective" assert "R1" in txt assert "ATP[c]" in txt +def test_txt_bounds_are_written_exactly(model, tmp_path): + reaction = model.reactions.get_by_id("R1") + reaction.bounds = (-1 / 3, 0.1) + reaction.objective_coefficient = 1e-7 + root = export_for_git(model, tmp_path, formats=("txt",)) + fields = (root / "txt" / "model.txt").read_text().splitlines()[1].split("\t") + assert fields[3:] == ["-0.3333333333333333", "0.1", "1e-07"] + assert (float(fields[3]), float(fields[4]), float(fields[5])) == (-1 / 3, 0.1, 1e-7) + assert not fields[1].endswith(" ") + + +def test_txt_whole_and_infinite_bounds(model, tmp_path): + model.reactions.get_by_id("R1").bounds = (-1000.0, float("inf")) + root = export_for_git(model, tmp_path, formats=("txt",)) + fields = (root / "txt" / "model.txt").read_text().splitlines()[1].split("\t") + assert fields[3:5] == ["-1000", "Inf"] + + +def test_all_five_formats_are_written_by_default(model, tmp_path): + root = export_for_git(model, tmp_path) + for fmt in ("yml", "xml", "mat", "xlsx", "txt"): + assert (root / fmt / f"model.{fmt}").exists() + + +def test_dependencies_are_ordered_and_repeatable(model, tmp_path): + set_model_metadata(model, version="1.2.0", date="2026-10-05") + first = export_for_git(model, tmp_path / "a", formats=("yml",)) + second = export_for_git(model, tmp_path / "b", formats=("yml",)) + text = (first / "dependencies.txt").read_text() + assert text == (second / "dependencies.txt").read_text() + names = [line.split("\t")[0] for line in text.splitlines()] + assert names == [ + "python", "cobra", "libSBML", "ruamel.yaml", "scipy", "openpyxl", "raven_toolbox", + "model_version", "model_date", + ] + assert "model_version\t1.2.0" in text + + def test_bad_format(model, tmp_path): with pytest.raises(ValueError, match="Unknown format"): export_for_git(model, tmp_path, formats=("yml", "json")) From c7af46a95e781d9e4571e8904501ce642db8d6b4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 11:03:19 +0200 Subject: [PATCH 03/71] Write equations without a coefficient of 1, as RAVEN does The txt table and the Excel file wrote '1 ATP[c] => 1 ADP[c]' where RAVEN writes 'ATP[c] => ADP[c]'. A coefficient other than 1 is still written. --- src/raven_toolbox/io/excel.py | 4 +++- tests/test_io_git.py | 7 +++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/raven_toolbox/io/excel.py b/src/raven_toolbox/io/excel.py index 0439dfd5..39305bb7 100644 --- a/src/raven_toolbox/io/excel.py +++ b/src/raven_toolbox/io/excel.py @@ -35,8 +35,10 @@ def _equation(rxn: cobra.Reaction) -> str: """Human-readable equation in RAVEN ``name[comp]`` form.""" def side(items): + # A coefficient of 1 is left out, as RAVEN's constructEquations does. return " + ".join( - f"{abs(coef):g} {met.name}[{met.compartment}]" for met, coef in items + f"{'' if abs(coef) == 1 else f'{abs(coef):g} '}{met.name}[{met.compartment}]" + for met, coef in items ) reactants = [(m, c) for m, c in rxn.metabolites.items() if c < 0] diff --git a/tests/test_io_git.py b/tests/test_io_git.py index 41602102..8a314ec1 100644 --- a/tests/test_io_git.py +++ b/tests/test_io_git.py @@ -245,3 +245,10 @@ def test_sbml_gene_label_is_the_gene_name(model, tmp_path): (g.id, g.name) for g in plain_back.genes ) assert back.reactions.R1.gene_reaction_rule == "ENSG1 or ENSG2" + + +def test_equation_leaves_out_a_coefficient_of_one(model, tmp_path): + model.reactions.get_by_id("R1").add_metabolites({model.metabolites.atp_c: -1}) + root = export_for_git(model, tmp_path, formats=("txt",)) + equation = (root / "txt" / "model.txt").read_text().splitlines()[1].split("\t")[1] + assert equation == "2 ATP[c] <=> ADP[c]" From 537cc53966808f2071f7a3028d52d18000772769 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 14:23:41 +0200 Subject: [PATCH 04/71] Connect blocked reactions with RAVEN's threshold, forced flux, net production and score handling A reaction is blocked when it cannot carry a flux above 1, the check is repeated with the templates added, 0.1 is forced through each rescuable reaction, net production relaxes the steady state in the checks as well as in the MILP, and scores that are not negative take the least negative penalty. --- src/raven_toolbox/gapfilling/fill.py | 106 ++++++++++++++++++--------- tests/test_gapfilling.py | 62 ++++++++++++++++ 2 files changed, 135 insertions(+), 33 deletions(-) diff --git a/src/raven_toolbox/gapfilling/fill.py b/src/raven_toolbox/gapfilling/fill.py index cc5bc795..adf1a8ed 100644 --- a/src/raven_toolbox/gapfilling/fill.py +++ b/src/raven_toolbox/gapfilling/fill.py @@ -11,7 +11,9 @@ :func:`add_reactions_from_model`), so templates in a different identifier namespace than the model still work. Per-reaction ``scores`` (higher = prefer to include) map to RAVEN's ``rxnScores``; the MILP minimises the penalty ``-score`` (default penalty -``1.0``, i.e. minimise the number of reactions added). +``1.0``, i.e. minimise the number of reactions added). As in RAVEN, a score that is not +negative is replaced by the least negative score among the template reactions (or by 0 +when none is negative), so a reaction is never rewarded for being added. """ from __future__ import annotations @@ -20,7 +22,7 @@ from dataclasses import dataclass import cobra -from cobra.flux_analysis import find_blocked_reactions, flux_variability_analysis +from cobra.flux_analysis import flux_variability_analysis from optlang.interface import Variable # untyped (resolves to Any); used for container hints from raven_toolbox.manipulation.transfer import add_reactions_from_model @@ -61,6 +63,50 @@ def _merge_templates(model: cobra.Model, templates: list[cobra.Model]) -> tuple[ return working, template_ids +def _carrying_flux( + model: cobra.Model, reaction_ids: list[str], *, threshold: float, net_production: bool +) -> set[str]: + """Ids among ``reaction_ids`` that can carry a flux larger than ``threshold`` in either direction. + + Infinite bounds count as +-10000, and with ``net_production`` the steady state is + relaxed to ``S @ v >= 0``, as in RAVEN's ``haveFlux`` on a model whose ``b`` upper + bound is infinite. + """ + if not reaction_ids: + return set() + with model as m: + for rxn in m.reactions: + if rxn.lower_bound == -math.inf: + rxn.lower_bound = -10000.0 + if rxn.upper_bound == math.inf: + rxn.upper_bound = 10000.0 + if net_production: + for met in m.metabolites: + m.constraints[met.id].ub = None + fva = flux_variability_analysis( + m, reaction_list=reaction_ids, fraction_of_optimum=0.0, processes=1 + ) + return { + rid for rid in reaction_ids + if rid in fva.index + and (fva.at[rid, "maximum"] > threshold or fva.at[rid, "minimum"] < -threshold) + } + + +def _penalties( + template_ids: list[str], scores: dict[str, float] | None, penalty: float +) -> dict[str, float]: + """Penalty per template reaction: ``-score`` for negative scores. + + A score that is not negative takes the penalty of the least negative score, or 0 when + no score is negative (RAVEN's ``getMinNrFluxes``). + """ + raw = {rid: scores[rid] if scores and rid in scores else -penalty for rid in template_ids} + negatives = [v for v in raw.values() if v < 0] + fallback = -max(negatives) if negatives else 0.0 + return {rid: -v if v < 0 else fallback for rid, v in raw.items()} + + def _solve_min_templates( working: cobra.Model, template_ids: list[str], @@ -103,11 +149,9 @@ def _solve_min_templates( for met in working.metabolites: working.constraints[met.id].ub = None - def pen(rid: str) -> float: - return -scores[rid] if scores and rid in scores else penalty - + penalties = _penalties(template_ids, scores, penalty) working.objective = prob.Objective( - sum(pen(rid) * indicators[rid] for rid in template_ids), direction="min" + sum(penalties[rid] * indicators[rid] for rid in template_ids), direction="min" ) working.slim_optimize() if working.solver.status != "optimal": @@ -133,16 +177,16 @@ def connect_blocked_reactions( scores: dict[str, float] | None = None, penalty: float = 1.0, allow_net_production: bool = False, - eps: float = 1.0, + threshold: float = 1.0, + min_flux: float = 0.1, ) -> GapFillResult: """Add template reactions so blocked draft reactions can carry flux. - Finds reactions that - cannot carry flux in ``model`` under its current bounds, then adds the minimum-penalty set of template - reactions that lets the (irreversible) ones carry flux, and returns the filled - model. Like RAVEN, only irreversible blocked reactions are forced — reversible - ones can carry flux trivially in the split formulation, so forcing them is - uninformative. + Finds reactions that cannot carry a flux larger than ``threshold`` in ``model`` under its current bounds, + then adds the minimum-penalty set of template reactions that lets the (irreversible) + ones carry flux, and returns the filled model. Like RAVEN, only irreversible blocked + reactions are forced — reversible ones can carry flux trivially in the split + formulation, so forcing them is uninformative. For the *other* gap-filling flavour — adding reactions to make the model's objective feasible — use ``cobra.flux_analysis.gapfill`` after aligning the @@ -151,39 +195,35 @@ def connect_blocked_reactions( The draft is expected to have exchange reactions for its nutrients (otherwise most reactions are trivially blocked). - ``eps`` is the minimum flux a blocked reaction must be able to carry to count as - rescuable, and the flux forced through it while solving. ``allow_net_production`` + ``threshold`` is the flux a reaction must be able to exceed to count as unblocked, in + the draft and, for a blocked reaction, with the templates added. ``min_flux`` is the + flux forced through each rescuable reaction while solving. ``allow_net_production`` relaxes the steady-state balance from ``S @ v == 0`` to ``S @ v >= 0`` (metabolites - may accumulate), which can rescue reactions whose real fix is a missing - consumption route rather than a genuine connectivity gap. + may accumulate) in both checks and in the MILP, which can rescue reactions whose real + fix is a missing consumption route rather than a genuine connectivity gap. """ if penalty <= 0: raise ValueError(f"penalty must be > 0 (a non-positive cost malforms the MILP); got {penalty}.") templates = _as_models(templates) - blocked = set(find_blocked_reactions(model, processes=1)) - candidates = [r for r in blocked if model.reactions.get_by_id(r).lower_bound >= 0] - working, template_ids = _merge_templates(model, templates) - target: list[str] = [] - if candidates: - fva = flux_variability_analysis( - working, reaction_list=candidates, fraction_of_optimum=0.0, processes=1 - ) - # A reaction can be missing from the FVA frame if the solver dropped it - # (e.g. the reaction was eliminated upstream); treat that as "unreachable" - # rather than letting the KeyError propagate. - target = [ - r for r in candidates - if r in fva.index and fva.at[r, "maximum"] > eps - ] + draft_ids = [r.id for r in model.reactions] + unblocked = _carrying_flux( + model, draft_ids, threshold=threshold, net_production=allow_net_production + ) + blocked = {rid for rid in draft_ids if rid not in unblocked} + candidates = [r for r in draft_ids if r in blocked and model.reactions.get_by_id(r).lower_bound >= 0] + rescuable = _carrying_flux( + working, candidates, threshold=threshold, net_production=allow_net_production + ) + target = [r for r in candidates if r in rescuable] cannot = sorted(blocked - set(target)) if not target: return GapFillResult([], [], cannot, model.copy()) for rid in target: - working.reactions.get_by_id(rid).lower_bound = eps + working.reactions.get_by_id(rid).lower_bound = min_flux chosen = _solve_min_templates( working, template_ids, scores=scores, penalty=penalty, allow_net_production=allow_net_production, diff --git a/tests/test_gapfilling.py b/tests/test_gapfilling.py index 3c37c576..f55bd18e 100644 --- a/tests/test_gapfilling.py +++ b/tests/test_gapfilling.py @@ -123,3 +123,65 @@ def test_unconnectable_reaction_reported_not_added(): assert res.added_reactions == [] assert res.newly_connected == [] assert "r1" in res.cannot_connect + + +# --------------------------------------------------------------------------- # +# RAVEN's blocked-reaction threshold, forced flux, net production and scores +# --------------------------------------------------------------------------- # +def _chain(exchange_ub, template_specs, draft_extra=()): + """Draft EX_A -> A -> B (R1) -> EX_B plus R2: C -> D and R3: E -> F with exits for D and F.""" + def rxn(rid, subs, prods, ub=1000.0): + r = cobra.Reaction(rid, lower_bound=0, upper_bound=ub) + r.add_metabolites({**{_met(s): -1 for s in subs}, **{_met(p): 1 for p in prods}}) + return r + + draft = cobra.Model("draft") + draft.add_reactions([ + rxn("EX_A", [], ["A"], exchange_ub), rxn("R1", ["A"], ["B"]), rxn("EX_B", ["B"], []), + rxn("R2", ["C"], ["D"]), rxn("EX_D", ["D"], []), + rxn("R3", ["E"], ["F"]), rxn("EX_F", ["F"], []), *draft_extra, + ]) + template = cobra.Model("template") + template.add_reactions([rxn(rid, subs, prods) for rid, subs, prods in template_specs]) + return draft, template + + +def test_flux_limit_below_threshold_counts_as_blocked(): + # Nothing can carry more than 0.5, so every draft reaction is blocked and none can be connected. + draft, template = _chain(0.5, [("T1", ["A"], ["C"]), ("T2", ["B"], ["E"])]) + res = connect_blocked_reactions(draft, template) + assert res.added_reactions == [] + assert set(res.cannot_connect) == {r.id for r in draft.reactions} + + +def test_forced_flux_is_small_enough_for_a_shared_supply(): + # R2 and R3 each need a template route from A; forcing 0.1 through each fits in the supply of 1.5. + draft, template = _chain(1.5, [("T1", ["A"], ["C"]), ("T2", ["A"], ["E"])]) + res = connect_blocked_reactions(draft, template) + assert res.added_reactions == ["T1", "T2"] + assert set(res.newly_connected) == {"EX_D", "EX_F", "R2", "R3"} + + +def test_net_production_applies_to_the_connectivity_check(): + # D has no exit: R2 can only carry flux when metabolites may accumulate. + def rxn(rid, subs, prods): + r = cobra.Reaction(rid, lower_bound=0, upper_bound=1000) + r.add_metabolites({**{_met(s): -1 for s in subs}, **{_met(p): 1 for p in prods}}) + return r + + draft = cobra.Model("draft") + draft.add_reactions([rxn("EX_A", [], ["A"]), rxn("R1", ["A"], ["B"]), rxn("EX_B", ["B"], []), rxn("R2", ["C"], ["D"])]) + template = cobra.Model("template") + template.add_reactions([rxn("T1", ["A"], ["C"])]) + assert connect_blocked_reactions(draft, template).cannot_connect == ["R2"] + res = connect_blocked_reactions(draft, template, allow_net_production=True) + assert res.added_reactions == ["T1"] + assert res.newly_connected == ["R2"] + + +def test_positive_score_takes_the_least_negative_penalty(): + # T3 (D -> B) scores +0.5: it is priced like the least negative score, not rewarded. + templates = [("T1", ["A"], ["C"]), ("T2", ["B"], ["C"]), ("T3", ["D"], ["B"])] + draft, template = _chain(1000.0, templates) + res = connect_blocked_reactions(draft, template, scores={"T1": -1.0, "T2": -1.0, "T3": 0.5}) + assert "T3" not in res.added_reactions From 91ff9a2fd793ceef60c73f30799d6d71593f6dbf Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 14:23:50 +0200 Subject: [PATCH 05/71] Score expression data like RAVEN's scoreModel A zero reference level scores 0/0 as -5 and a positive level as +10, a missing level gives the gene no score, and a reaction whose genes have no score takes the no-gene score instead of NaN. --- src/raven_toolbox/init/score.py | 26 ++++++++++++++++++-------- tests/test_init_score.py | 17 +++++++++++++++++ 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/src/raven_toolbox/init/score.py b/src/raven_toolbox/init/score.py index d0b9b32c..5c7557e1 100644 --- a/src/raven_toolbox/init/score.py +++ b/src/raven_toolbox/init/score.py @@ -4,8 +4,8 @@ to per-reaction scores by walking each reaction's GPR: genes joined by **OR** (isozymes) are combined with ``isozyme_scoring`` (default ``max``); genes joined by **AND** (complexes) with ``complex_scoring`` (default ``min``). Genes missing from -``gene_scores`` are *omitted*; a reaction with no genes — or whose genes are all -missing — gets ``no_gene_score`` (default −2). These reaction scores feed +``gene_scores``, or scored NaN, are *omitted*; a reaction with no genes — or whose genes are all +missing or NaN — gets ``no_gene_score`` (default −2). These reaction scores feed :func:`raven_toolbox.init.ftinit`. Upstream — the omics-data → gene-score step (thresholding, expression levels) — lives @@ -36,13 +36,22 @@ def gene_scores_from_expression( HPA are alternative upstream sources). ``reference`` is either a per-gene reference level (e.g. the cross-sample mean) or a single threshold for all genes: a gene expressed above its reference scores positive, below it negative. The - score is clamped to ``[min_score, max_score]``; non-positive level/reference (and - missing reference) → ``min_score`` (RAVEN maps these NaNs to -5). + score is clamped to ``[min_score, max_score]``. A zero level against a positive + reference scores ``min_score``; a positive level against a zero reference scores + ``max_score``; zero against zero, a negative level or reference, and a missing + reference score ``min_score``. A gene whose level is missing (``None`` or NaN) gets + no entry, so :func:`score_reactions_from_genes` leaves it out of the rule. """ scores: dict[str, float] = {} for gene, level in expression.items(): + if level is None or math.isnan(level): + continue ref = reference if isinstance(reference, (int, float)) else reference.get(gene) - if not level or not ref or level <= 0 or ref <= 0: + if ref is None or math.isnan(ref) or level < 0 or ref < 0 or (level == 0 and ref == 0): + scores[gene] = min_score + elif ref == 0: + scores[gene] = max_score + elif level == 0: scores[gene] = min_score else: scores[gene] = max(min(factor * math.log(level / ref), max_score), min_score) @@ -51,7 +60,8 @@ def gene_scores_from_expression( def _score_node(node, gene_scores: Mapping[str, float], iso, cplx) -> float | None: if isinstance(node, ast.Name): - return gene_scores.get(node.id) # None if the gene has no score + score = gene_scores.get(node.id) + return None if score is None or math.isnan(score) else score # None if the gene has no score if isinstance(node, ast.BoolOp): agg = iso if isinstance(node.op, ast.Or) else cplx vals = [s for v in node.values if (s := _score_node(v, gene_scores, iso, cplx)) is not None] @@ -71,8 +81,8 @@ def score_reactions_from_genes( Isozymes (GPR **OR**) are combined with ``isozyme_scoring`` (default ``max``); complex subunits (**AND**) with ``complex_scoring`` (default ``min``). A reaction - with no genes, or whose genes are all missing from ``gene_scores``, gets - ``no_gene_score`` (default -2). + with no genes, or whose genes are all missing from ``gene_scores`` or scored NaN, + gets ``no_gene_score`` (default -2). """ iso, cplx = resolve_aggregators(isozyme_scoring, complex_scoring) diff --git a/tests/test_init_score.py b/tests/test_init_score.py index 8c9a0c12..4a818ddf 100644 --- a/tests/test_init_score.py +++ b/tests/test_init_score.py @@ -67,3 +67,20 @@ def test_expression_per_gene_reference(): expr = {"g": 20.0} s = gene_scores_from_expression(expr, {"g": 5.0}) assert s["g"] == pytest.approx(5 * math.log(4)) + + +def test_expression_zero_reference_scores_like_raven(): + s = gene_scores_from_expression({"zero": 0.0, "one": 1.0, "ten": 10.0}, 0.0) + assert s == {"zero": -5.0, "one": 10.0, "ten": 10.0} + + +def test_expression_missing_level_gets_no_score(): + s = gene_scores_from_expression({"g": math.nan, "h": None, "k": 2.0}, 1.0) # type: ignore[dict-item] + assert set(s) == {"k"} + + +def test_nan_gene_score_is_left_out_of_the_rule(gpr_model): + scores = score_reactions_from_genes(gpr_model, {"g1": math.nan, "g2": math.nan, "g3": math.nan}) + assert set(scores.values()) == {-2.0} + partial = score_reactions_from_genes(gpr_model, {"g1": 1.0, "g2": math.nan, "g3": 3.0}) + assert not any(math.isnan(v) for v in partial.values()) From 8a63b01fa8f982aa3b8ddea3949fafa755124121 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 22 Oct 2026 14:26:44 +0200 Subject: [PATCH 06/71] Fill gaps with a MILP the way RAVEN's gapFillMILP does The growth floor is one constraint on the objective, a template reaction with a positive lower bound is only forced when it is selected, finite bounds are used as the coupling constants, and a template reaction with the id of a draft reaction is renamed. --- src/raven_toolbox/gapfilling/milp.py | 109 ++++++++++++++++----------- tests/test_gapfilling_milp.py | 79 +++++++++++++++++++ 2 files changed, 145 insertions(+), 43 deletions(-) diff --git a/src/raven_toolbox/gapfilling/milp.py b/src/raven_toolbox/gapfilling/milp.py index aa11f4eb..91ab7614 100644 --- a/src/raven_toolbox/gapfilling/milp.py +++ b/src/raven_toolbox/gapfilling/milp.py @@ -2,9 +2,15 @@ Finds the smallest set of template reactions that makes the objective (biomass) >= *min_growth*, following the database repair of Kumar et al. 2007 -(BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212). Reactions already in -the draft keep their bounds: a reaction with the wrong direction is a curation -issue and is not changed here. +(BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212). The floor applies to the +objective as a whole (the sum of coefficient times flux), as in RAVEN. Reactions +already in the draft keep their bounds: a reaction with the wrong direction is a +curation issue and is not changed here. + +A template reaction that has the id of a draft reaction (or of an earlier template +reaction) is renamed to ``_