diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b8810aa..64ee5c2f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,70 @@ 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, 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. +* `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`). 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`). +* `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). +* `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 + `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 + 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; + `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. 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`). +* `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. +* 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. * `analysis.fseof` (and so `geckopy.ec_fseof`) solves each level with pFBA that holds the objective at 0.999999 of its optimum at that level, as RAVEN's `solveLP` does, where it held the full value. On the scipy solver (HiGHS) the optimum is only accurate to the solver's feasibility tolerance, so the second step was infeasible at the top level of ecYeastGEM and the scan returned 9 of 10 levels; it now returns all 10, as on Gurobi and GLPK. A level at which the biomass maximisation is infeasible still ends the scan. * `analysis.fseof` takes the top of the scan as `max_fraction` of the target flux in the pFBA solution that maximises the target, as RAVEN's `FSEOF` and GECKO's `ecFSEOF` do, where it took the plain maximum (about 1e-6 relative difference). * `analysis.fseof` lists targets of equal score in model order (a stable sort), as MATLAB's `sort` does, where their order was arbitrary. @@ -23,6 +87,25 @@ 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. + 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. +* `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. ## 3.0.0b1 — 2026-09-07 diff --git a/scripts/benchmark_certified_multiorg.py b/scripts/benchmark_certified_multiorg.py index 2e8cf14f..f472fb81 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/scripts/benchmark_replicate_yeast_gem.py b/scripts/benchmark_replicate_yeast_gem.py index a72549e5..091e0c44 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.""" diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index d12d8bf0..6163cca8 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -12,7 +12,9 @@ import cobra import numpy as np import pandas as pd +from cobra.core.gene import GPR +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 @@ -177,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 @@ -434,75 +440,14 @@ 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]: @@ -518,28 +463,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() - 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: - continue - 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: @@ -909,7 +838,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. @@ -965,7 +894,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/src/raven_toolbox/gapfilling/topological.py b/src/raven_toolbox/gapfilling/topological.py index af3d8a83..6d146012 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/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index bd671260..2a725cfd 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -142,6 +142,8 @@ 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, + _pin_fluxes: bool = True, metabolomics: Mapping[str, Iterable[str]] | None = None, prod_weight: float = 0.5, seed: int = _EXTRACT_SEED, @@ -187,7 +189,21 @@ 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`. + 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`), 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 + reaction id of each merge group, so the order does not depend on which member names + the merged reaction. + + ``_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 @@ -502,11 +518,27 @@ 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. + # 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) + 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): + fluxes = _read_solution()[1] kept = free_or_essential | on deleted = [r.id for r in model.reactions if r.id not in kept] @@ -536,7 +568,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 @@ -546,8 +578,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. @@ -556,9 +589,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``, @@ -613,7 +646,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: @@ -629,6 +662,145 @@ def _phase(objective, label: str) -> bool: return ok +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 + 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; 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 ``(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, 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()] + 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, 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]])) + 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) + free_sign.append(rid) + opt.add(extra) + 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 + if time_limit is not None: + 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": + 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, + ) + 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 + + +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. 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), 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: + 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))} + 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 + 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() + _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, + ) + 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]: """Push tiny reaction scores off zero (RAVEN ``ftINIT.m:160-161``). @@ -707,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, - 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, @@ -728,7 +900,8 @@ 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, + _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, @@ -860,7 +1033,14 @@ 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; among those, it prefers the + 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 @@ -925,6 +1105,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: @@ -944,7 +1132,10 @@ 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, + # 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, diff --git a/src/raven_toolbox/init/genes.py b/src/raven_toolbox/init/genes.py index fa78a089..a964d2c0 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/merge.py b/src/raven_toolbox/init/merge.py index ce87fcf5..670e0d47 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/src/raven_toolbox/init/taskfill.py b/src/raven_toolbox/init/taskfill.py index 0d0bd7f5..b502df34 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 ddfd046f..8f8678a8 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 9ccad62c..18b9993d 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 @@ -60,6 +61,7 @@ ec_data_from_yaml_sections, ec_data_to_yaml_sections, ) +from raven_toolbox.utils.gpr import set_gene_rule # A dedicated instance for *writing* only (reading uses _read_yaml below). # Mutating cobra's own YAML instance would change its width for @@ -248,6 +250,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: @@ -255,10 +260,20 @@ 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 +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 read_yaml_model(path: str | Path) -> cobra.Model: """Read a RAVEN/cobrapy YAML model into a ``cobra.Model``. @@ -378,8 +393,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/src/raven_toolbox/localization/assign.py b/src/raven_toolbox/localization/assign.py index 38cf74ec..29fb56f2 100644 --- a/src/raven_toolbox/localization/assign.py +++ b/src/raven_toolbox/localization/assign.py @@ -20,21 +20,51 @@ """ from __future__ import annotations -from collections.abc import Callable, Mapping +from collections.abc import Callable, Iterable, Mapping from dataclasses import dataclass, field 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"] # --------------------------------------------------------------------------- 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 _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: comps = {m.compartment for m in rxn.metabolites if m.compartment} return next(iter(comps)) if len(comps) == 1 else None @@ -122,10 +152,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. + 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() 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. @@ -138,10 +172,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 @@ -174,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): @@ -193,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/localization/certify.py b/src/raven_toolbox/localization/certify.py index 6e7f6421..fd43cd16 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -42,6 +42,7 @@ AssignmentProposal, GrowthCondition, _base_met, + _check_unique_keys, _rxn_compartment, apply_assignment, ) @@ -649,6 +650,7 @@ def assign_compartments( if default_compartment not in compartments: raise ValueError(f"default_compartment={default_compartment!r} not in {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 b03dbaa1..0ec0e46b 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, + _base_met, + _check_unique_keys, + _rxn_compartment, +) from raven_toolbox.localization.scores import LocalizationScores __all__ = ["curation_priority"] @@ -171,6 +176,7 @@ def curation_priority( import pandas as pd 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 1cad0053..3e254e15 100644 --- a/src/raven_toolbox/localization/relocate.py +++ b/src/raven_toolbox/localization/relocate.py @@ -28,6 +28,7 @@ from raven_toolbox.localization.assign import ( AssignmentProposal, _base_met, + _check_unique_keys, apply_assignment, ) from raven_toolbox.localization.curation import _is_impermeant @@ -117,6 +118,7 @@ def relocate_reactions( 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 + _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/add.py b/src/raven_toolbox/manipulation/add.py index a76bc7bf..b546b620 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 01ab3f38..6718fbc5 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 1911fe3c..a051c4c8 100644 --- a/src/raven_toolbox/manipulation/compartments.py +++ b/src/raven_toolbox/manipulation/compartments.py @@ -9,14 +9,21 @@ """ from __future__ import annotations +import copy import warnings from collections.abc import Callable, Iterable import cobra +from raven_toolbox.utils.gpr import set_gene_rule + # 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, @@ -24,7 +31,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. @@ -47,22 +54,48 @@ 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, 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. + 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`. 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 @@ -74,22 +107,26 @@ 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. - canonical: dict[str, cobra.Metabolite] = {} - used_ids: set[str] = set() + # 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. + protos = {} 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] + 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 = dict(proto.notes or {}) + 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 # Rewrite all reactions: replace each metabolite with its canonical, summing @@ -136,7 +173,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 {}) @@ -212,7 +249,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 {}) @@ -231,22 +268,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/src/raven_toolbox/manipulation/expand.py b/src/raven_toolbox/manipulation/expand.py index f509ea41..fe77d56c 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, set_gene_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) @@ -107,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 8e6f8b39..3283e977 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 83665624..504e0146 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/remove.py b/src/raven_toolbox/manipulation/remove.py index 43f966fb..074cc180 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,22 @@ 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. + 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: drop metabolites *and* genes orphaned by the removal. + passes through to cobra: also drop metabolites orphaned by the removal. Returns ------- @@ -97,6 +108,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 +129,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/src/raven_toolbox/manipulation/transfer.py b/src/raven_toolbox/manipulation/transfer.py index 4f744c05..c2570e42 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/blast.py b/src/raven_toolbox/reconstruction/homology/blast.py index 8ffe1e08..afbbb7e4 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/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index e8aea823..7ce6abb1 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/homology/supporth.py b/src/raven_toolbox/reconstruction/homology/supporth.py index ecd1ac09..56d7bd32 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/assemble.py b/src/raven_toolbox/reconstruction/kegg/assemble.py index 5c86775e..44d0c849 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/reconstruction/kegg/hmm.py b/src/raven_toolbox/reconstruction/kegg/hmm.py index 0fe0ed1b..df16624e 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/src/raven_toolbox/reconstruction/kegg/query.py b/src/raven_toolbox/reconstruction/kegg/query.py index f949325c..b49e0172 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/src/raven_toolbox/utils/__init__.py b/src/raven_toolbox/utils/__init__.py index 3919e0a0..58a97bbb 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 86e352ba..a2f3cb30 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,138 @@ 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. +# 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", "|", "||"}) - ``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, 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. """ - if not rule.strip(): + 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 rule is None or _RULE_BLANK.fullmatch(rule): return None - gene_logger = logging.getLogger("cobra.core.gene") - was_disabled = gene_logger.disabled - gene_logger.disabled = True + 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)) + + +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: - 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 + 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) diff --git a/tests/test_curation.py b/tests/test_curation.py index d4d5c5b3..654d0a0b 100644 --- a/tests/test_curation.py +++ b/tests/test_curation.py @@ -1285,6 +1285,44 @@ 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} + + +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_gapfilling_topological.py b/tests/test_gapfilling_topological.py index f8273854..12666560 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"} diff --git a/tests/test_gene_rule_copy.py b/tests/test_gene_rule_copy.py new file mode 100644 index 00000000..de64d064 --- /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 diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 581f5dba..89f80daa 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -220,6 +220,134 @@ 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) + + +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) + + +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.""" + 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. # --------------------------------------------------------------------------- # @@ -314,7 +442,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 diff --git a/tests/test_init_merge.py b/tests/test_init_merge.py index 028a8208..0604c652 100644 --- a/tests/test_init_merge.py +++ b/tests/test_init_merge.py @@ -99,6 +99,44 @@ 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_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()) diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 8d6b89d9..ab736e9b 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,79 @@ 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). 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}, -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, ""), + } + 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"}) + + +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 + + 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 run_spy(*args, **kwargs): + seen.append((step[0], kwargs["_pin_fluxes"])) + return real_run(*args, **kwargs) + + 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} + + # --------------------------------------------------------------------------- # # metabolomics (production-bonus for detected metabolites) — full pipeline. # Exercises what run_ftinit's own metabolomics tests cannot: metabolite *name* diff --git a/tests/test_io_yaml.py b/tests/test_io_yaml.py index d8265480..390f8c7b 100644 --- a/tests/test_io_yaml.py +++ b/tests/test_io_yaml.py @@ -532,3 +532,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 diff --git a/tests/test_localization_certify.py b/tests/test_localization_certify.py index c31883ef..680f253e 100644 --- a/tests/test_localization_certify.py +++ b/tests/test_localization_certify.py @@ -431,3 +431,142 @@ 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 + + +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 269f12e1..72bc8a62 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. @@ -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(): @@ -159,8 +163,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 +174,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 +203,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 +216,95 @@ 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"} + + +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", "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(): diff --git a/tests/test_manipulation_expand.py b/tests/test_manipulation_expand.py index a692ef3f..21c12900 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_manipulation_remove.py b/tests/test_manipulation_remove.py index 3bfe0bef..3162c5c0 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"] diff --git a/tests/test_reconstruction_blast.py b/tests/test_reconstruction_blast.py index 184150d8..30b8c23c 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" diff --git a/tests/test_reconstruction_kegg_hmm.py b/tests/test_reconstruction_kegg_hmm.py index 48bba63a..b1107c4f 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_kegg_query.py b/tests/test_reconstruction_kegg_query.py index bce7003b..6992c3c1 100644 --- a/tests/test_reconstruction_kegg_query.py +++ b/tests/test_reconstruction_kegg_query.py @@ -174,3 +174,46 @@ 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 _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 + + 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): + calls.append(cmd) + open(cmd[cmd.index("--tblout") + 1], "w").close() + 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" diff --git a/tests/test_reconstruction_supporth.py b/tests/test_reconstruction_supporth.py index 7e2e3795..6ff1a067 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") diff --git a/tests/test_utils_gpr.py b/tests/test_utils_gpr.py index 68f00f7d..0c3d2ab9 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,78 @@ 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) + + +# 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