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..f72d019 100644 --- a/docs/studies/multiorganism_validation.md +++ b/docs/studies/multiorganism_validation.md @@ -14,10 +14,20 @@ only a compartment-label map and a metabolite base-key per organism. | organism | kingdom | model | reactions | compartments | chloroplast | wall | certified | reaction agr | gene agr | transports | blocked | |---|---|---|--:|--:|:--:|--:|:--:|--:|--:|--:|--:| -| *S. cerevisiae* | fungal | yeast-GEM | 2569 | 14 | no | ~140 s | **yes** | 72.0 % | **88.7 %** | 1001 | 29.7 % | -| *H. sapiens* | mammalian | Human-GEM | 12854 | 9 | no | 190–464 s | **yes** | 52.8 % | 79.4 % | 3701 | 10.7 % | -| *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 % | +| *S. cerevisiae* | fungal | yeast-GEM | 2569 | 14 | no | ~140 s | **yes** | 72.5 % | **88.7 %** | 967 | 29.7 % | +| *H. sapiens* | mammalian | Human-GEM | 12877 | 9 | no | 190–464 s | **yes** | 56.9 % | 79.4 % | 3608 | 10.7 % | +| *A. thaliana* | plant | AraCore v2.1 | 585 | 6 | **yes** | 35 s | **yes** | **82.1 %** | **93.7 %** | 200 | 0.0 % | +| *C. reinhardtii* | algal | iCre1355 | 2394 | 11 | **yes** | 93 s | **yes** | 54.1 % | 70.5 % | 1083 | 32.9 % | + +All rows are measured on the **deterministic score-aligned placement** (each reaction placed by its own +enzymes' DeepLoc scores; see [yeast_validation.md](yeast_validation.md)), not the earlier arbitrary +co-optimum. **Gene** agreement is unchanged from the arbitrary-placement runs — the tie-break fixes only +the reaction placement, leaving the gene layout untouched — which is why every gene-agreement figure +above matches its pre-tie-break value. Reaction agreement rose or held (Human-GEM 52.8 % → 56.9 % the +largest move; AraCore and iCre1355 were already near their ceilings), and coherent placement needs fewer +transports across the board. The one non-obvious shift is iCre1355's blocked fraction (19.0 % → 32.9 %): +tighter, less transport-heavy placement leaves more reactions dead-ended, though the model still +certifies. Placement is reproducible run-to-run. 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; @@ -26,22 +36,22 @@ Chloroplast Protein Atlas, UniProt experimental across all) replaces it for the ## What P1 establishes -1. **It scales past 4 compartments at genome size.** Human-GEM — 12854 reactions, **9 compartments**, - 5206 reactions re-placed — certifies (materialised biomass above the floor by a real FBA) in +1. **It scales past 4 compartments at genome size.** Human-GEM — 12877 reactions, **9 compartments**, + 5221 reactions re-placed — certifies (materialised biomass above the floor by a real FBA) in ~3–8 min. iCre1355 reaches **11 compartments** (the richest organelle set here) in ~1.5 min. The flux-free master + FBA-certification design does not degrade as compartments and reactions grow; this was the open scaling risk, now closed. 2. **It handles the chloroplast** — the generality hook CarveFungi structurally cannot represent. DeepLoc's "Plastid" label maps to the chloroplast code (`h` in both plant models), and both AraCore - (plant) and iCre1355 (alga) place and certify with a chloroplast, AraCore at **81.8 % reaction / + (plant) and iCre1355 (alga) place and certify with a chloroplast, AraCore at **82.1 % reaction / 93.7 % gene** self-agreement. 3. **It is organism- and predictor-agnostic in practice, not just in principle.** The same function, the same DeepLoc-2 evidence, no algorithm changes — only a `{DeepLoc-label -> compartment-code}` map and a metabolite base-key (name-based for yeast/Human-GEM; bracket-stripped `Glc[c]->Glc` for AraCore; suffix-stripped `10fthf_c->10fthf` for iCre1355) — carry it from a 585-reaction plant core - model to a 12854-reaction human GEM. + model to a 12877-reaction human GEM. 4. **Every model stays functional.** All four materialise a growing model (blocked-reaction fractions - 0–30 %, mostly the models' intrinsic dead-ends), the property that distinguishes this method. + 0–33 %, mostly the models' intrinsic dead-ends), the property that distinguishes this method. ## Honest caveats 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: