From 2da2af0ab76fc8befd8e6dbb1a0a685d4ea67558 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 5 Nov 2026 12:44:27 +0100 Subject: [PATCH 1/2] Drop the complexes without an ortholog under the drop policy The policy drop removes each complex that has a subunit without an ortholog and removes the reaction when no complex remains. It removed the reaction as soon as one subunit had no ortholog. --- CHANGELOG.md | 4 +++ .../reconstruction/homology/homology.py | 34 +++++++++---------- tests/test_reconstruction_homology.py | 25 ++++++++++++-- 3 files changed, 43 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a975bb1..448ee4c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,10 @@ ## Unreleased +### Breaking changes + +* `get_model_from_homology` with `complex_policy="drop"` removes each complex that has a subunit without an ortholog, and removes the reaction only when no complex remains (was: the reaction was removed as soon as one subunit had no ortholog). For the rule `(A and B) or (C and D)` without an ortholog of `B`, the draft has `C and D`. + ### Data * The KEGG assets are the `kegg120` release of `toolbox-assets` (was: `kegg118`): the reference model and tables, the taxonomy and the profile HMM libraries of 4246 prokaryotic and 2724 eukaryotic KOs. The baked registry and `data/manifest.json` list them. diff --git a/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index 2239f16..6e631ae 100644 --- a/src/raven_toolbox/reconstruction/homology/homology.py +++ b/src/raven_toolbox/reconstruction/homology/homology.py @@ -91,27 +91,23 @@ def _rewrite_node(node, ortho: dict): return _group("and", children) -def _contains_unmapped(expr) -> bool: - if isinstance(expr, _Unmapped): - return True - return isinstance(expr, tuple) and any(_contains_unmapped(c) for c in expr[1]) - - def _apply_policy(expr, policy: str, model_id: str): """Resolve the unmapped leaves of a rewritten tree under ``complex_policy``. - ``flag`` turns each into an ``OLD__`` placeholder gene, ``keep`` removes - them (a group left empty disappears with them) and ``drop`` removes the whole reaction - as soon as one is left. Returns ``None`` when nothing remains. + ``flag`` turns each into an ``OLD__`` placeholder gene and ``keep`` removes + them (a group left empty disappears with them). ``drop`` removes each complex that has + one: an AND group with an unmapped subunit is removed as a whole, and an OR group keeps + its remaining alternatives. Returns ``None`` when nothing remains. """ - if policy == "drop": - return None if _contains_unmapped(expr) else expr if isinstance(expr, _Unmapped): return f"OLD_{model_id}_{expr.gene}" if policy == "flag" else None if isinstance(expr, str): return expr op, children = expr - return _group(op, [c for c in (_apply_policy(c, policy, model_id) for c in children) if c is not None]) + resolved = [_apply_policy(c, policy, model_id) for c in children] + if policy == "drop" and op == "and" and any(c is None for c in resolved): + return None + return _group(op, [c for c in resolved if c is not None]) def _prune_genes(node, removed: set): @@ -334,11 +330,15 @@ def get_model_from_homology( complex_policy : str What happens to an AND-linked subunit with no ortholog. Default ``"flag"``. ``"flag"`` keeps the reaction with a placeholder gene (``OLD__``) for - later curator review, ``"keep"`` removes that subunit, ``"drop"`` drops the whole - reaction. An isozyme group counts as a subunit without an ortholog only when none of - its members has one, and then one placeholder stands for it. ``"drop"`` removes the - reaction whenever a placeholder would remain, also when another isozyme branch is - complete. Any other value raises :class:`ValueError`. + later curator review, ``"keep"`` removes that subunit, ``"drop"`` removes each complex + that has such a subunit and removes the reaction when no complex remains. An isozyme + group counts as a subunit without an ortholog only when none of its members has one, + and then one placeholder stands for it. Under ``"drop"``, the rule + ``(A and B) or (C and D)`` with no ortholog of ``B`` becomes ``C and D``, and the rule + ``(A or B) and C`` with no ortholog of ``A`` or ``B`` removes the reaction. A complex + that contains another complete complex is kept, and a gene shared by several complexes + removes only the complexes in which it has no ortholog. Any other value raises + :class:`ValueError`. only_genes_in_models : bool Restricts the hits table to hits whose template gene appears in the given ``models`` before mapping orthologs, so a gene the templates lack cannot be a best hit. It raises diff --git a/tests/test_reconstruction_homology.py b/tests/test_reconstruction_homology.py index d2d9d0e..fef710d 100644 --- a/tests/test_reconstruction_homology.py +++ b/tests/test_reconstruction_homology.py @@ -130,9 +130,28 @@ def test_flagged_subunit_stays_beside_a_complete_isozyme_branch(): assert _policy_draft("A or (B and C)", ["A", "C"], "keep") == {"R1": "nA or nC"} -def test_complex_policy_drop_removes_reaction_with_any_flagged_subunit(): - assert _policy_draft("A or (B and C)", ["A", "C"], "drop") == {} - assert _policy_draft("(A and B) or C", ["C"], "drop") == {} +def test_complex_policy_drop_removes_each_complex_with_a_flagged_subunit(): + assert _policy_draft("A or (B and C)", ["A", "C"], "drop") == {"R1": "nA"} + assert _policy_draft("(A and B) or C", ["C"], "drop") == {"R1": "nC"} + assert _policy_draft("A and B", ["A"], "drop") == {} + assert _policy_draft("(A and B) or (C and D)", ["A", "C", "D"], "drop") == {"R1": "nC and nD"} + + +def test_complex_policy_drop_keeps_complete_complexes_of_a_shared_subunit(): + rule = "(A and B) or (A and C) or (D and E) or (A and B and F)" + assert _policy_draft(rule, ["A", "B", "C", "D", "F"], "drop") == { + "R1": "(nA and nB) or (nA and nC) or (nA and nB and nF)"} + assert _policy_draft("(A and B) or B", ["B"], "drop") == {"R1": "nB"} + assert _policy_draft("(A and B) or B", ["A", "B"], "drop") == {"R1": "(nA and nB) or nB"} + + +def test_complex_policy_drop_keeps_isozyme_groups_inside_a_complex(): + hits = make_ortholog_hits([("A", "n1"), ("A", "n2"), ("B", "n3")], "T", "N") + draft = get_model_from_homology([_one_rule_template("T", "R1", "A and B")], hits, "N", complex_policy="drop") + assert _rules(draft) == {"R1": "(n1 or n2) and n3"} + + +def test_complex_policy_flag_and_keep_on_an_incomplete_complex_beside_an_isozyme(): assert _policy_draft("(A and B) or C", ["C"], "flag") == {"R1": "(OLD_T_A and OLD_T_B) or nC"} assert _policy_draft("(A and B) or C", ["C"], "keep") == {"R1": "nC"} From e5af32f02b8e7451b60c8539f2a5bf6f0de8857d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 6 Nov 2026 13:54:48 +0100 Subject: [PATCH 2/2] Treat a gene taken by an earlier template as a subunit without an ortholog In the preferred-order step, a template gene that maps to a gene an earlier template already maps loses its mappings and follows complex_policy. For the rule (a and b) or x with a taken, flag gives (OLD_T2_a and nb) or nx. A reaction with no gene that keeps a mapping is not transferred. --- CHANGELOG.md | 1 + .../reconstruction/homology/homology.py | 73 +++++++------------ tests/test_reconstruction_homology.py | 51 ++++++++++++- 3 files changed, 75 insertions(+), 50 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 448ee4c..a1c7ace 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Breaking changes +* `get_model_from_homology` with `preferred_order` treats a template gene that maps to a gene an earlier template already maps as a subunit without an ortholog, and `complex_policy` decides what happens to it (was: the gene and every part of the rule that needs it were removed). For the rule `(a and b) or x` with `a` taken by an earlier template, `flag` gives `(OLD_T2_a and nb) or nx`, `keep` gives `nb or nx` and `drop` gives `nx`. A reaction that has genes, none of which keeps a mapping, is not transferred. * `get_model_from_homology` with `complex_policy="drop"` removes each complex that has a subunit without an ortholog, and removes the reaction only when no complex remains (was: the reaction was removed as soon as one subunit had no ortholog). For the rule `(A and B) or (C and D)` without an ortholog of `B`, the draft has `C and D`. ### Data diff --git a/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index 6e631ae..f9c8d4d 100644 --- a/src/raven_toolbox/reconstruction/homology/homology.py +++ b/src/raven_toolbox/reconstruction/homology/homology.py @@ -110,21 +110,6 @@ def _apply_policy(expr, policy: str, model_id: str): return _group(op, [c for c in resolved if c is not None]) -def _prune_genes(node, removed: set): - """Remove template genes from a GPR AST, as RAVEN's removeGenes does. - - Returns the pruned tree, or ``None`` when the rule cannot be satisfied without them - (a removed gene in an AND, or every branch of an OR removed). - """ - if isinstance(node, ast.Name): - return None if node.id in removed else node - kids = [_prune_genes(c, removed) for c in node.values] - if isinstance(node.op, ast.And): - return None if any(k is None for k in kids) else node - kept = [k for k in kids if k is not None] - return ast.BoolOp(op=ast.Or(), values=kept) if kept else None - - def _group(op: str, children: list): """Flatten same-operator groups, drop duplicates, and unwrap a single child.""" flat: list = [] @@ -146,23 +131,18 @@ def _render(expr, parent: str | None = None) -> str: return f"({text})" if parent is not None and parent != op else text -def _rewrite_gpr(rxn, ortho: dict, policy: str, model_id: str, removed: frozenset = frozenset()): +def _rewrite_gpr(rxn, ortho: dict, policy: str, model_id: str): """Return the rewritten GPR string, or None if the reaction is unsupported. - ``removed`` holds template genes taken out of the template before the rewrite (genes a - preferred model already covers); a reaction that cannot work without them is unsupported. + A template gene that ``ortho`` does not map, including one that a preferred model already + covers, is a subunit without an ortholog and follows ``policy``. """ if not rxn.gene_reaction_rule: return None # A reaction is only transferred if at least one of its genes has an ortholog. - if not any(g.id in ortho or g.id in removed for g in rxn.genes): + if not any(g.id in ortho for g in rxn.genes): return None - body = rxn.gpr.body - if removed: - body = _prune_genes(body, set(removed)) - if body is None: - return None - result = _apply_policy(_rewrite_node(body, ortho), policy, model_id) + result = _apply_policy(_rewrite_node(rxn.gpr.body, ortho), policy, model_id) return None if result is None else _render(result) @@ -235,26 +215,23 @@ def _ortholog_map( return ortho -def _apply_preferred_order(ortho: dict, order: list[str]) -> tuple[dict, dict]: +def _apply_preferred_order(ortho: dict, order: list[str]) -> dict: """Give each new gene to the first model (in order) that maps it, as RAVEN does. A template gene of a later model that maps to any new gene an earlier model already - maps is removed from that model altogether, together with its other mappings. The - new genes of every model count as taken, including those of its removed genes. - Returns the pruned mapping and ``{model_id: {removed template genes}}``. + maps loses all its mappings and counts as a gene without an ortholog. The new genes of + every model count as taken, including those of the genes that lose their mappings. + Returns the pruned mapping. """ taken: set = set() pruned: dict = {mid: {} for mid in ortho} - removed: dict = {mid: set() for mid in ortho} for position, model_id in enumerate(order): per_model = ortho.get(model_id, {}) for template_gene, new_genes in per_model.items(): - if position > 0 and taken.intersection(new_genes): - removed[model_id].add(template_gene) - else: + if not (position > 0 and taken.intersection(new_genes)): pruned[model_id][template_gene] = new_genes taken.update(ng for new_genes in per_model.values() for ng in new_genes) - return pruned, removed + return pruned def get_model_from_homology( @@ -303,11 +280,14 @@ def get_model_from_homology( preferred_order : Sequence[str] or None Breaks ties when more than one template maps the same new-organism gene: the gene is mapped to the earlier template in the list, and a gene of a later template that maps - to any gene an earlier template already maps is removed from that template - altogether, with its other mappings. A reaction of the later template that cannot work - without the removed gene is not transferred. Templates are merged in the order they - were given in ``models``, whatever the preferred order is, and a template that - ``preferred_order`` leaves out contributes nothing. The tie-break applies only when + to any gene an earlier template already maps loses all its mappings. It counts as a + subunit without an ortholog and follows ``complex_policy``: for the rule + ``(a and b) or x`` with ``a`` taken by an earlier template, ``"flag"`` gives + ``(OLD__a and nb) or nx``, ``"keep"`` gives ``nb or nx`` and ``"drop"`` gives + ``nx``. A reaction of the later template that has genes, none of which keeps a mapping, is + not transferred. Templates are merged in the order they were given in ``models``, + whatever the preferred order is, and a template that ``preferred_order`` leaves out + contributes nothing. The tie-break applies only when more than one model is given. Default ``None``, which uses the order of ``models``. RAVEN's ``preferredOrder`` decides which template receives a gene in the same way and merges the models in the order of ``models``. @@ -473,12 +453,11 @@ def get_model_from_homology( ) order = [str(x) for x in preferred_order] if preferred_order else model_ids - removed: dict = {} if preferred_order and len(models) > 1: - ortho, removed = _apply_preferred_order(ortho, order) + ortho = _apply_preferred_order(ortho, order) draft = _transfer(model_by_id, order, ortho, model_for, model_ids, complex_policy, - keep_gene_free, preserve_notes, removed) + keep_gene_free, preserve_notes) candidates = None if review_identity is not None: @@ -491,12 +470,11 @@ def get_model_from_homology( model_genes=model_genes, max_evalue=max_evalue, min_align_len=min_align_len, min_identity=review_identity, ) - loose_removed: dict = {} if preferred_order and len(models) > 1: - loose_ortho, loose_removed = _apply_preferred_order(loose_ortho, order) + loose_ortho = _apply_preferred_order(loose_ortho, order) loose = _transfer( model_by_id, order, loose_ortho, model_for, model_ids, complex_policy, - keep_gene_free, preserve_notes, loose_removed, + keep_gene_free, preserve_notes, ) extra = {r.id for r in loose.reactions} - {r.id for r in draft.reactions} candidates = _candidate_evidence( @@ -509,14 +487,13 @@ def get_model_from_homology( def _transfer(model_by_id, order, ortho, model_for, model_ids, complex_policy, - keep_gene_free=False, preserve_notes=False, removed=None) -> cobra.Model: + keep_gene_free=False, preserve_notes=False) -> cobra.Model: """Assemble the draft from the reactions of each template whose GPRs survive rewriting. The templates are merged in the order of ``model_ids`` (the order they were given), whatever the preferred order is; reaction ids shared between templates are resolved by that merge. Templates left out of ``order`` contribute nothing. """ - removed = removed or {} transferred = [] for mid in model_ids: model = model_by_id.get(mid) @@ -526,7 +503,7 @@ def _transfer(model_by_id, order, ortho, model_for, model_ids, complex_policy, m = model.copy() keep: dict[str, str] = {} for rxn in m.reactions: - new_gpr = _rewrite_gpr(rxn, per_model, complex_policy, mid, frozenset(removed.get(mid, ()))) + new_gpr = _rewrite_gpr(rxn, per_model, complex_policy, mid) if new_gpr is None and keep_gene_free and not rxn.gene_reaction_rule: new_gpr = "" if new_gpr is not None: diff --git a/tests/test_reconstruction_homology.py b/tests/test_reconstruction_homology.py index fef710d..6222a4e 100644 --- a/tests/test_reconstruction_homology.py +++ b/tests/test_reconstruction_homology.py @@ -272,7 +272,7 @@ def test_preferred_order_without_shared_genes_merges_in_given_order(): assert plain["R1"] == "N1" -def test_preferred_order_keeps_a_reaction_that_survives_without_the_removed_gene(): +def test_preferred_order_retains_the_orphan_subunit_of_a_complex(): t1 = _one_rule_template("T1", "R_T1", "Y") t2 = _one_rule_template("T2", "R_T2", "X or Z") add_reactions_from_equations(t2, [{"id": "R_T3", "equation": "b_c --> a_c", "gene_reaction_rule": "X and Z"}]) @@ -281,8 +281,55 @@ def test_preferred_order_keeps_a_reaction_that_survives_without_the_removed_gene make_ortholog_hits([("X", "N1"), ("Z", "N3")], "T2", "N")], ignore_index=True, ) + flag = get_model_from_homology([t1, t2], hits, "N", preferred_order=["T1", "T2"]) + assert _rules(flag) == {"R_T1": "N1", "R_T2": "N3", "R_T3": "OLD_T2_X and N3"} + keep = get_model_from_homology([t1, t2], hits, "N", preferred_order=["T1", "T2"], complex_policy="keep") + assert _rules(keep) == {"R_T1": "N1", "R_T2": "N3", "R_T3": "N3"} + drop = get_model_from_homology([t1, t2], hits, "N", preferred_order=["T1", "T2"], complex_policy="drop") + assert _rules(drop) == {"R_T1": "N1", "R_T2": "N3"} + + +def test_preferred_order_follows_the_policy_for_a_taken_gene_in_a_complex_beside_an_isozyme(): + t1 = _one_rule_template("T1", "R_T1", "Y") + t2 = _one_rule_template("T2", "R_T2", "(a and b) or x") + hits = pd.concat( + [make_ortholog_hits([("Y", "N1")], "T1", "N"), + make_ortholog_hits([("a", "N1"), ("b", "Nb"), ("x", "Nx")], "T2", "N")], + ignore_index=True, + ) + rules = { + policy: _rules(get_model_from_homology([t1, t2], hits, "N", preferred_order=["T1", "T2"], + complex_policy=policy))["R_T2"] + for policy in ("flag", "keep", "drop") + } + assert rules == {"flag": "(OLD_T2_a and Nb) or Nx", "keep": "Nb or Nx", "drop": "Nx"} + + +def test_preferred_order_chains_the_conflicts_of_three_templates(): + t1 = _one_rule_template("T1", "R_T1", "Y") + t2 = _one_rule_template("T2", "R_T2", "X") + t3 = _one_rule_template("T3", "R_T3", "(p and q) or r") + hits = pd.concat( + [make_ortholog_hits([("Y", "N1")], "T1", "N"), make_ortholog_hits([("X", "N2")], "T2", "N"), + make_ortholog_hits([("p", "N1"), ("q", "N2"), ("r", "N3")], "T3", "N")], + ignore_index=True, + ) + flag = get_model_from_homology([t1, t2, t3], hits, "N", preferred_order=["T1", "T2", "T3"]) + assert _rules(flag)["R_T3"] == "(OLD_T3_p and OLD_T3_q) or N3" + keep = get_model_from_homology([t1, t2, t3], hits, "N", preferred_order=["T1", "T2", "T3"], + complex_policy="keep") + assert _rules(keep)["R_T3"] == "N3" + + +def test_preferred_order_does_not_transfer_a_reaction_whose_genes_were_all_taken(): + t1 = _one_rule_template("T1", "R_T1", "Y") + t2 = _one_rule_template("T2", "R_T2", "a and b") + hits = pd.concat( + [make_ortholog_hits([("Y", "N1")], "T1", "N"), make_ortholog_hits([("a", "N1")], "T2", "N")], + ignore_index=True, + ) res = get_model_from_homology([t1, t2], hits, "N", preferred_order=["T1", "T2"]) - assert _rules(res) == {"R_T1": "N1", "R_T2": "N3"} # R_T3 needs the removed X + assert _rules(res) == {"R_T1": "N1"} # --- gene-free reactions, notes, duplicate genes, metadata ------------------------