Treat a gene taken by an earlier template as a subunit without an ortholog

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

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_<model>_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:

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 ------------------------