From 3f9465e845cc226cfe8fc94df8b0f084b9322bce Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 19 Jul 2026 23:59:13 +0200 Subject: [PATCH] Deterministic score-aligned reaction placement in the assignment master The placement master maximises a gene-localisation objective that never mentions the per-reaction placement variable, so each reaction's compartment was a free co-optimum: the pinned solver returned it reproducibly but arbitrarily, giving 52.8% reaction agreement with curated yeast-GEM (an earlier un-pinned build happened to land ~72%). Add a lexicographic second pass: fix the gene layout to the primary optimum (fixing the y solution, not the objective value, so there is no tolerance to tune), then place each reaction in the compartment its own enzymes are predicted to occupy -- the summed DeepLoc score of the reaction's genes, with a small default_compartment prior so genes-free and score-tied reactions fall there deterministically. Yeast reaction agreement rises to 72.5% (1408/1943) and now rests on the localisation evidence rather than a solver tie-break; gene agreement is unchanged (88.7%, 716/807 -- the gene layout is untouched); coherent placement adds fewer transports (1001 -> 967); growth and blocked fraction unchanged; the warm-started second solve adds ~4s. Reproducible across independent runs. Docs updated (yeast_validation, localization_redesign, multiorganism_validation footnote); the multi-organism reaction-agreement rows await re-measurement. --- CHANGELOG.md | 13 ++++++ docs/studies/localization_redesign.md | 7 +++ docs/studies/multiorganism_validation.md | 7 +++ docs/studies/yeast_validation.md | 16 +++++-- src/raven_toolbox/localization/certify.py | 54 ++++++++++++++++++++--- 5 files changed, 88 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f51ea24..50fa6fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,19 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* **Deterministic, score-aligned reaction placement in `assign_compartments`.** The placement master + maximises a gene-localisation objective that never mentions the per-reaction placement variable, so + each reaction's compartment was a free co-optimum: the pinned solver returned it reproducibly but + *arbitrarily*, and reaction agreement with curated yeast-GEM was only **52.8 %** (an earlier un-pinned + build happened to land ~72 %). A lexicographic second pass now fixes the gene layout to the primary + optimum, then places each reaction in the compartment its own enzymes are predicted to occupy — the + summed DeepLoc score of the reaction's genes, with a small `default_compartment` prior so genes-free + and score-tied reactions fall there deterministically. Yeast reaction agreement rises to **72.5 %** + (1408/1943) and now rests on the localisation evidence rather than a solver tie-break; **gene agreement + is unchanged** (88.7 %, 716/807 — the gene layout is untouched); coherent placement fragments fewer + metabolites, so it adds **fewer transports** (1001 → 967); growth and blocked-fraction are unchanged; + the second solve is warm-started and adds ~4 s. Reproducible run-to-run (identical placement across + independent runs). * **`diff_models` compares grRules as logic, not text.** The GPR check now DNF-expands each rule (via the existing `manipulation.gpr_to_dnf`), sorts the genes within each isozyme clause and sorts the clauses, so operand order no longer registers as a difference: `a and b` == `b and a` and `a or b` == `b or a`. The diff --git a/docs/studies/localization_redesign.md b/docs/studies/localization_redesign.md index 966caa6..a2a5e75 100644 --- a/docs/studies/localization_redesign.md +++ b/docs/studies/localization_redesign.md @@ -182,6 +182,13 @@ Reading: - **The honest agreement-vs-transport frontier is now measurable.** 72.5 % agreement costs ~1261 reconnecting transports (983 FVA-usable), against curated yeast-GEM's own 1467 inter-compartment reactions — high agreement honestly requires many transports (there is no free lunch). +- **Reaction placement is deterministic and score-aligned, not an arbitrary co-optimum.** The master's + objective is over gene-localisation variables only, so each reaction's compartment is a free + co-optimum — a pinned solver returns it reproducibly, but *arbitrarily* (an un-tie-broken solve scored + only 52.8 %). A lexicographic second pass fixes the gene layout to the primary optimum, then places + each reaction in the compartment its own enzymes are predicted to occupy (summed DeepLoc score of the + reaction's genes, small default-compartment prior for genes-free and tied reactions). The gene layout + is untouched, so the 72.5 % rests on the localisation evidence itself, not on a solver's tie-break. Every number is authoritative because the growth column *is* the certificate: there is no separate MILP floor to disagree with it. Raw JSON in `.research_tmp/certify_yeast_gem.json`. diff --git a/docs/studies/multiorganism_validation.md b/docs/studies/multiorganism_validation.md index 9e72501..b6e514e 100644 --- a/docs/studies/multiorganism_validation.md +++ b/docs/studies/multiorganism_validation.md @@ -19,6 +19,13 @@ only a compartment-label map and a metabolite base-key per organism. | *A. thaliana* | plant | AraCore v2.1 | 585 | 6 | **yes** | 35 s | **yes** | **81.8 %** | **93.7 %** | 215 | 0.0 % | | *C. reinhardtii* | algal | iCre1355 | 2394 | 11 | **yes** | 93 s | **yes** | 53.4 % | 70.5 % | 1108 | 19.0 % | +> **Reaction-agreement column predates the deterministic placement tie-break** (Unreleased; see +> [yeast_validation.md](yeast_validation.md)). Before it, the per-reaction compartment was a free +> co-optimum the solver resolved arbitrarily; it is now placed by each reaction's own enzyme scores. +> These reaction-agreement figures will be re-measured on the new placement (yeast, re-run, went 72.0 % +> → 72.5 %). The **gene**-agreement column, growth, and blocked fraction are unaffected — the tie-break +> leaves the gene layout untouched. + Agreement is scored against each model's **own** compartment annotation (a *circular* truth — the model being re-placed is also the reference). This measures self-consistency and functionality, not accuracy; the independent-ground-truth work (P2: HPA for human, SUBA5 for Arabidopsis, the Chlamydomonas diff --git a/docs/studies/yeast_validation.md b/docs/studies/yeast_validation.md index a5cc672..044fa6e 100644 --- a/docs/studies/yeast_validation.md +++ b/docs/studies/yeast_validation.md @@ -21,19 +21,29 @@ parent organelle, matching DeepLoc's resolution). | metric | `assign_compartments` | |---|--:| -| reaction-level agreement | **72.0 %** (1400/1943) | +| reaction-level agreement | **72.5 %** (1408/1943) | | gene-level agreement | **88.7 %** (716/807) | -| transports added | 1001 (curated 1467) | +| transports added | 967 (curated 1467) | | blocked re-placed reactions | 29.7 % | | materialised growth | 0.1426 /h (floor 0.040) | | certified (real FBA) | **yes** | Gene agreement (88.7 %) sits well above DeepLoc's ~64 % raw ceiling because functionality corrects -placements; the model stays functional because the method provisions the ~1000 transports the placement +placements; the model stays functional because the method provisions the ~970 transports the placement actually needs, and the blocked 29.7 % is yeast-GEM's *intrinsic* dead-end floor, not a placement artefact. Every number is authoritative because the growth column is a real FBA on the materialised model — the certificate, not a separate solver floor. +Reaction placement is **deterministic and score-aligned.** The placement master maximises a +gene-localisation objective, which never mentions the per-reaction placement variable — so each +reaction's compartment was a free co-optimum: reproducible once the solver is pinned, but an *arbitrary* +vertex (which on its own scored only 52.8 %). A lexicographic second pass fixes the gene layout to the +primary optimum, then places each reaction in the compartment its own enzymes are predicted to occupy +(the summed DeepLoc score of the reaction's genes, with a small default-compartment prior for genes-free +and score-tied reactions). The gene layout — hence gene agreement — is untouched; the 72.5 % reaction +agreement now rests on the localisation evidence itself, not on a solver's arbitrary tie-break, and the +coherent placement needs fewer transports (967 vs 1001). + ## 2. Head-to-head with CarveFungi (Comparison 2) CarveFungi's own `minmax_reduction` carve (CPLEX, unmodified) on its real universal-DB candidate set + diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index a76bcce..9ae76fe 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -158,8 +158,14 @@ def _score(score_df, g: str, c: str) -> float: return 0.0 if v is None or (isinstance(v, float) and math.isnan(v)) else float(v) +# Tie-break weight for the placement master's second pass: small enough never to override a real +# per-reaction score difference (DeepLoc scores are ~[0, 1] summed over a reaction's genes), just enough +# to send a genes-free or score-tied reaction to the default compartment deterministically. +_DEFAULT_COMPARTMENT_PRIOR = 1e-3 + + def _solve_placement_master( - model, movable, genes_in_scope, gene_rxns, score_df, compartments, *, + model, movable, genes_in_scope, gene_rxns, score_df, compartments, default_compartment, *, multi_compartment_penalty, forced, colocation_groups, time_limit, ): """Flux-free score-maximising placement MILP (mono-localisation). @@ -169,6 +175,11 @@ def _solve_placement_master( compartment (``x[a, c] = x[b, c]``), letting the score objective pick which one. No flux and no growth floor at all in the master — the placement can never harvest a compartment's score through leaked flux. + + Solved in two lexicographic passes: the primary objective places genes by score (and penalises + spread), which leaves the per-reaction placement ``x`` a free co-optimum; the second pass then fixes + the gene layout and places each reaction in the compartment its own enzymes score highest for, so the + reaction placement is deterministic and evidence-aligned rather than an arbitrary co-optimal vertex. """ model.solver # noqa: B018 — initialise the solver so model.problem is usable prob = model.problem @@ -212,21 +223,52 @@ def _solve_placement_master( opt.add(list(x.values()) + list(y.values()) + cons) - obj = [] + # Primary objective: place each gene in the compartment(s) its DeepLoc score favours, penalising + # spread across compartments. This fixes the gene layout but leaves the *reaction* placement x a free + # co-optimum (the objective never mentions x), which a deterministic solver resolves to an arbitrary + # vertex -- yeast-GEM reaction agreement 52.8%. A lexicographically-lower second pass then chooses x. + primary = [] for g in genes_sorted: for c in compartments: s = _score(score_df, g, c) if s: - obj.append(mul([Real(s), y[g, c]])) + primary.append(mul([Real(s), y[g, c]])) if multi_compartment_penalty: for v in y.values(): - obj.append(mul([Real(-multi_compartment_penalty), v])) - opt.objective = prob.Objective(add(obj) if obj else Real(0.0), direction="max") + primary.append(mul([Real(-multi_compartment_penalty), v])) + opt.objective = prob.Objective(add(primary) if primary else Real(0.0), direction="max") if time_limit is not None: opt.configuration.timeout = int(time_limit) _pin_deterministic(prob, opt) opt.optimize() + if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit"): + return opt.status, {}, {} + # Lexicographic second pass. Fix the gene layout to the primary optimum -- fixing the *solution* (each + # y binary), not the objective value, so there is no near-optimal tolerance to tune -- then place each + # reaction in the compartment its own enzymes are predicted to occupy: reward x[r, c] by the summed + # gene score of r's genes for c, with a small prior for `default_compartment` so genes-free and + # score-tied reactions fall there deterministically rather than to an arbitrary co-optimal vertex. + # The gene layout (agreement, multi-compartment consolidation) is untouched; only reaction placement, + # which was free, is now meaningful. + try: # read every primal before touching any bound -- the first bound change discards the solution + y_star = {k: 1.0 if (v.primal or 0.0) >= 0.5 else 0.0 for k, v in y.items()} + except (AttributeError, ValueError): + return "no_incumbent", {}, {} + for k, v in y.items(): + v.lb = v.ub = y_star[k] + secondary = [] + for r in movable: + r_genes = sorted({gg.id for gg in r.genes} & genes_in_scope) + for c in compartments: + w = sum(_score(score_df, g, c) for g in r_genes) + if c == default_compartment: + w += _DEFAULT_COMPARTMENT_PRIOR + if w: + secondary.append(mul([Real(w), x[r.id, c]])) + opt.objective = prob.Objective(add(secondary) if secondary else Real(0.0), direction="max") + _pin_deterministic(prob, opt) + opt.optimize() if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit"): return opt.status, {}, {} try: @@ -621,7 +663,7 @@ def assign_compartments( seen.add(signature) status, placements, gene_comps = _solve_placement_master( - model, movable, genes_in_scope, gene_rxns, score_df, compartments, + model, movable, genes_in_scope, gene_rxns, score_df, compartments, default_compartment, multi_compartment_penalty=multi_compartment_penalty, forced=forced, colocation_groups=groups, time_limit=time_limit) if not placements: