From c210f21df28367f31e107031cd799a74dbc20bf0 Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Sun, 27 Sep 2026 22:56:39 -0400 Subject: [PATCH 1/7] pools.csv + pool_transitions.csv (deltasignal specs/039): interconversion pools found at NODE level (the loop A->F->B->R->A must exist, so uuid-separated copies are not re-merged); base = resting form (the modified side is the product of the donor-consuming direction: ATP/GTP/SAM/acetyl-CoA/NAD+/Ub; else more components; else more modified residues). Additive files; network unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 125 +++++++++++++++++++++++++++++++++ src/neo4j_connector.py | 65 +++++++++++++++++ src/pathway_generator.py | 14 ++++ tests/test_pools.py | 58 +++++++++++++++ 4 files changed, 262 insertions(+) create mode 100644 tests/test_pools.py diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 25ea14b..102121f 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3409,6 +3409,131 @@ def parts(node_id: str): logger.info(f"Exported {len(rows)} drug-derived entities to {output_file}") +def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, + uuid_mapping: Dict[str, str], pairs, mod_counts: Dict[str, Tuple[int, int]], + donor_reactions: Optional[Set[str]] = None): + """Interconversion pools at NODE level (deltasignal specs/039). + + A curated pair (F: A -> B, R: B -> A) is a pool transition only where the + same A node feeds a copy of F and receives a copy of R's output, and the + same B node is produced by that F and consumed by that R: the loop + A -> F -> B -> R -> A exists in this network. Occurrences the per-position + uuids separated are not re-merged. Forms joined by such loops form one pool. + + Returns (forms, transitions): + forms: [(pool_id, node_uuid, stable_id, is_base)] + transitions: [(pool_id, from_uuid, to_uuid, reaction_uuid, reaction_stid)] + The base form (the unmodified, resting form) is decided per pair, then per pool: + the modified form is the product of the direction that consumes a group donor + (ATP, GTP, SAM, acetyl-CoA, NAD+, ubiquitin); failing that, the form with more + components (a complex against its parts), then more modified residues. The + pool's base is the form no pair marks as modified; ties fall back to the + smaller stId and are counted by the caller. + """ + donor_reactions = donor_reactions or set() + if pathway_logic_network.empty or not pairs: + return [], [] + ent = _uuid_to_stable_id_map(pathway_logic_network, uuid_mapping) + vr = {str(u): str(r) for u, r in zip(reaction_id_map["uid"], reaction_id_map["reactome_id"])} + ins: Dict[str, Set[str]] = {} + outs: Dict[str, Set[str]] = {} + for _, e in pathway_logic_network.iterrows(): + s, t, et = str(e["source_id"]), str(e["target_id"]), e.get("edge_type") + if et == "input" and t in vr: + ins.setdefault(t, set()).add(s) + elif et == "output" and s in vr: + outs.setdefault(s, set()).add(t) + copies: Dict[str, List[str]] = {} + for u, r in vr.items(): + copies.setdefault(r, []).append(u) + trans = set() + for f, r, a, b in pairs: + for fu in sorted(copies.get(f, [])): + for ru in sorted(copies.get(r, [])): + for au in sorted(ins.get(fu, set()) & outs.get(ru, set())): + if ent.get(au) != a: + continue + for bu in sorted(outs.get(fu, set()) & ins.get(ru, set())): + if ent.get(bu) != b: + continue + trans.add((au, bu, fu, f)) + trans.add((bu, au, ru, r)) + if not trans: + return [], [] + parent: Dict[str, str] = {} + + def find(x): + parent.setdefault(x, x) + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + for a_u, b_u, _, _ in trans: + ra, rb = find(a_u), find(b_u) + if ra != rb: + parent[max(ra, rb)] = min(ra, rb) + groups: Dict[str, List[str]] = {} + for n in sorted(parent): + groups.setdefault(find(n), []).append(n) + + def modified_side(a_u, b_u, f_stid, r_stid): + """The node the pair marks as modified, or None when undecidable.""" + fd, rd = f_stid in donor_reactions, r_stid in donor_reactions + if fd != rd: + return b_u if fd else a_u + ma, ca = mod_counts.get(ent.get(a_u, ""), (0, 0)) + mb, cb = mod_counts.get(ent.get(b_u, ""), (0, 0)) + if ca != cb: + return b_u if cb > ca else a_u + if ma != mb: + return b_u if mb > ma else a_u + return None + marked: Dict[str, int] = {} + fwd = {(a_u, b_u): rs for a_u, b_u, _, rs in trans} + for (a_u, b_u), fs in fwd.items(): + rs = fwd.get((b_u, a_u)) + if rs is None or a_u > b_u: + continue + side = modified_side(a_u, b_u, fs, rs) + if side is not None: + marked[side] = marked.get(side, 0) + 1 + + def base_key(u): + s = ent.get(u, "") + mods, comps = mod_counts.get(s, (0, 0)) + return (marked.get(u, 0), mods, comps, s, u) + forms, transitions = [], [] + order = sorted(groups.values(), key=lambda g: min(ent.get(u, "") + u for u in g)) + for i, members in enumerate(order, start=1): + pid = f"pool{i}" + base = min(members, key=base_key) + for u in sorted(members): + forms.append((pid, u, ent.get(u, ""), u == base)) + for a_u, b_u, rx, rs in sorted(trans): + if a_u in members: + transitions.append((pid, a_u, b_u, rx, rs)) + return forms, transitions + + +def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, + uuid_mapping: Dict[str, str], pools_file: str, transitions_file: str) -> None: + """Write pools.csv and pool_transitions.csv (deltasignal specs/039): each + protein's interconversion pool, for a consumer that solves it at steady state + (pi = pi P) instead of multiplying around the loop. Both files are written + even when empty, so a consumer can tell "no pools" from "no table".""" + from src.neo4j_connector import (get_interconversion_pairs, get_form_modification_counts, + get_donor_reactions) + pairs = get_interconversion_pairs(pathway_id) + mods = get_form_modification_counts({x for p in pairs for x in (p[2], p[3])}) + donors = get_donor_reactions({x for p in pairs for x in (p[0], p[1])}, _UBIQUITIN_STIDS) + forms, transitions = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, pairs, mods, donors) + pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "is_base"]).to_csv(pools_file, index=False) + pd.DataFrame(transitions, columns=["pool_id", "from_uuid", "to_uuid", "reaction_uuid", "reaction_stid"]).to_csv( + transitions_file, index=False) + logger.info(f"Exported {len({f[0] for f in forms})} pools ({len(forms)} forms, " + f"{len(transitions)} transitions) to {pools_file}") + + def export_node_resolution(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, diff --git a/src/neo4j_connector.py b/src/neo4j_connector.py index cdc3b45..b9399f1 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -932,6 +932,71 @@ def derived(s: str, depth: int = 0) -> bool: for s in ids if derived(s)} +# Group donors: a reaction consuming one ADDS a group to its substrate, so its +# product is the modified form of an interconversion pool (specs/039). ChEBI +# ids: ATP, GTP, S-adenosyl-L-methionine, acetyl-CoA, NAD+. +_DONOR_CHEBI = ("30616", "37565", "15414", "59789", "15351", "57288", "15846", "57540") + + +def get_donor_reactions(reaction_ids, ubiquitin_stids=()) -> Set[str]: + """Reactions (stIds) that consume a group donor or ubiquitin.""" + ids = sorted({r for r in reaction_ids if r}) + if not ids: + return set() + query = """ + UNWIND $ids AS rid MATCH (r:ReactionLikeEvent {stId: rid})-[:input]->(i:PhysicalEntity) + OPTIONAL MATCH (i)-[:referenceEntity]->(re:ReferenceMolecule) + WITH rid, i, re WHERE re.identifier IN $chebi OR i.stId IN $ub + RETURN DISTINCT rid + """ + try: + rows = get_graph().run(query, ids=ids, chebi=list(_DONOR_CHEBI), ub=list(ubiquitin_stids)).data() + except Exception: + logger.error("Error in get_donor_reactions", **_traceback_kwargs()) + raise + return {r["rid"] for r in rows} + + +def get_interconversion_pairs(pathway_id: str) -> List[Tuple[str, str, str, str]]: + """Curated interconversion pairs in a pathway (deltasignal specs/039): + (forward reaction, reverse reaction, A, B) with F: A -> B and R: B -> A, both + events of the pathway, A and B distinct non-small-molecule entities. + Returned as stIds, sorted, each unordered pair once per orientation.""" + query = """ + MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(r1:ReactionLikeEvent)-[:input]->(a:PhysicalEntity), + (r1)-[:output]->(b:PhysicalEntity) + MATCH (p)-[:hasEvent*]->(r2:ReactionLikeEvent)-[:input]->(b), (r2)-[:output]->(a) + WHERE r1 <> r2 AND a <> b AND NOT a:SimpleEntity AND NOT b:SimpleEntity + RETURN DISTINCT r1.stId AS f, r2.stId AS r, a.stId AS a, b.stId AS b + """ + try: + rows = get_graph().run(query, pid=pathway_id).data() + except Exception: + logger.error("Error in get_interconversion_pairs", **_traceback_kwargs()) + raise + return sorted((r["f"], r["r"], r["a"], r["b"]) for r in rows) + + +def get_form_modification_counts(stable_ids) -> Dict[str, Tuple[int, int]]: + """stId -> (modified residues summed over its protein leaves, component + count), to pick a pool's least-modified (base) form (deltasignal specs/039).""" + ids = sorted({s for s in stable_ids if s}) + if not ids: + return {} + query = """ + UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) + RETURN s, + size([(e)-[:hasComponent|hasMember|hasCandidate*0..6]->(l:EntityWithAccessionedSequence)-[:hasModifiedResidue]->(m) | m]) AS mods, + size([(e)-[:hasComponent*1..6]->(c) | c]) AS comps + """ + try: + rows = get_graph().run(query, ids=ids).data() + except Exception: + logger.error("Error in get_form_modification_counts", **_traceback_kwargs()) + raise + return {r["s"]: (int(r["mods"]), int(r["comps"])) for r in rows} + + def get_reactome_release() -> Optional[int]: """The Reactome release number of the connected graph, or None. diff --git a/src/pathway_generator.py b/src/pathway_generator.py index 8b9d873..94be301 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -12,6 +12,7 @@ create_pathway_logic_network, export_cofactors, export_drugs, + export_pools, export_containment, export_entity_reaction_proxy_mapping, export_node_reaction_context, @@ -510,6 +511,19 @@ def generate_pathway_file( except Exception as e: logger.error(f"Failed to write drugs.csv: {e}", exc_info=True) + # Interconversion pools (deltasignal specs/039), in their own block. + try: + export_pools( + pathway_id, + result.logic_network, + result.reaction_id_map, + result.uuid_mapping, + str(pathway_output_dir / "pools.csv"), + str(pathway_output_dir / "pool_transitions.csv"), + ) + except Exception as e: + logger.error(f"Failed to write pools.csv: {e}", exc_info=True) + logger.info(f"Output directory: {pathway_output_dir}") except (ConnectionError, ValueError) as e: diff --git a/tests/test_pools.py b/tests/test_pools.py new file mode 100644 index 0000000..5a92358 --- /dev/null +++ b/tests/test_pools.py @@ -0,0 +1,58 @@ +"""Interconversion pools (deltasignal specs/039): each protein's forms joined by +curated forward/reverse reactions, found at NODE level and oriented so the base +is the resting (unmodified) form. Shape: RAS. A = RAS:GDP, B = RAS:GTP; +F = "GEF exchange" consumes GTP (a donor), R = "GAP hydrolysis".""" +import pandas as pd + +import src.logic_network_generator as m + +# _uuid_to_stable_id_map detects the mapping direction by uuid SHAPE +U = {k: f"00000000-0000-0000-0000-{i:012d}" for i, k in enumerate(["a", "a2", "b", "vf", "vr"], 1)} + +A, B, F, R = "R-HSA-GDP", "R-HSA-GTP", "R-HSA-F", "R-HSA-R" +PAIRS = [(F, R, A, B), (R, F, B, A)] + + +def net(loop=True): + # a -> f -> b -> r -> a (the node-level loop); with loop=False, R's output is a + # different copy of A (a2), which the uuids keep apart + a_out = "a" if loop else "a2" + e = lambda s, t, et: {"source_id": U[s], "target_id": U[t], "edge_type": et} + edges = pd.DataFrame([e("a", "vf", "input"), e("vf", "b", "output"), + e("b", "vr", "input"), e("vr", a_out, "output")]) + rmap = pd.DataFrame({"uid": [U["vf"], U["vr"]], "reactome_id": [F, R]}) + umap = {U["a"]: A, U["a2"]: A, U["b"]: B} + return edges, rmap, umap + + +def test_a_node_level_loop_is_one_pool_with_the_resting_form_as_base(): + edges, rmap, umap = net() + forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {F}) + assert sorted((u, base) for _, u, _, base in forms) == sorted([(U["a"], True), (U["b"], False)]) + assert sorted((f, t, rx) for _, f, t, rx, _ in trans) == sorted([(U["a"], U["b"], U["vf"]), (U["b"], U["a"], U["vr"])]) + assert {p for p, *_ in forms} == {"pool1"} + + +def test_copies_the_uuids_separated_are_not_pooled(): + edges, rmap, umap = net(loop=False) + forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}) + assert forms == [] and trans == [] + + +def test_orientation_by_donor_then_components_then_residues(): + edges, rmap, umap = net() + # the donor sits on the REVERSE reaction: then GDP is the modified side + forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {R}) + assert dict((u, b) for _, u, _, b in forms) == {U["a"]: False, U["b"]: True} + # no donor: the form with more components is the modified one (a complex) + forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 3)}, set()) + assert dict((u, b) for _, u, _, b in forms)[U["a"]] is True + # no donor, same components: more modified residues is modified + forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (2, 1), B: (0, 1)}, set()) + assert dict((u, b) for _, u, _, b in forms)[U["b"]] is True + + +def test_no_pairs_or_empty_network_give_no_pools(): + edges, rmap, umap = net() + assert m.find_pools(edges, rmap, umap, [], {}, set()) == ([], []) + assert m.find_pools(pd.DataFrame(columns=["source_id", "target_id", "edge_type"]), rmap, umap, PAIRS, {}, set()) == ([], []) From 1358d039ed526cfe1c80b50eca5d502c64d3616d Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Sun, 27 Sep 2026 23:22:39 -0400 Subject: [PATCH 2/7] specs/039 review fixes: pools are one-protein interconversions only - get_interconversion_pairs: each reaction's only non-small-molecule input is the source form and only output the other form (ubiquitin allowed as a co-substrate), and the two forms share a reference entity. Excludes enzyme binding cycles (E + S -> E:S), release, degradation and two-form reactions, which were 78 of the 102 shipped pools. - find_pools: a reaction node converting two form pairs is dropped with its loops (its node would be written by two fluxes); counted in stats. - orientation: residues first (as pre-registered), then donor, then components (specs/039 amendment 1); undecided pools counted as ties. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 50 ++++++++++++++++++++++++---------- src/neo4j_connector.py | 25 ++++++++++++++--- tests/test_pools.py | 36 +++++++++++++++++++----- 3 files changed, 86 insertions(+), 25 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 102121f..cb5deac 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3411,7 +3411,7 @@ def parts(node_id: str): def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, uuid_mapping: Dict[str, str], pairs, mod_counts: Dict[str, Tuple[int, int]], - donor_reactions: Optional[Set[str]] = None): + donor_reactions: Optional[Set[str]] = None, stats: Optional[Dict[str, int]] = None): """Interconversion pools at NODE level (deltasignal specs/039). A curated pair (F: A -> B, R: B -> A) is a pool transition only where the @@ -3423,13 +3423,22 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram Returns (forms, transitions): forms: [(pool_id, node_uuid, stable_id, is_base)] transitions: [(pool_id, from_uuid, to_uuid, reaction_uuid, reaction_stid)] - The base form (the unmodified, resting form) is decided per pair, then per pool: - the modified form is the product of the direction that consumes a group donor - (ATP, GTP, SAM, acetyl-CoA, NAD+, ubiquitin); failing that, the form with more - components (a complex against its parts), then more modified residues. The - pool's base is the form no pair marks as modified; ties fall back to the - smaller stId and are counted by the caller. + A reaction node that would be a transition more than once (it converts two + different form pairs) is dropped, with every loop through it: its node would + be written by two fluxes (review of specs/039, finding 2). + + The base form (the unmodified, resting form) is decided per pair, then per pool + (specs/039 amendment 1): the modified form is the one with more modified + residues (as pre-registered); failing that, the product of the direction that + consumes a group donor (ATP, GTP, SAM, acetyl-CoA, NAD+, ubiquitin: RAS:GDP + and RAS:GTP carry no residue difference); failing that, the form with more + components. The pool's base is the form no pair marks as modified; a pool no + pair decides falls back to the smaller stId and is counted in + ``stats["ties"]``, with the dropped reactions in ``stats["multi_use"]``. """ + stats = stats if stats is not None else {} + stats["ties"] = 0 + stats["multi_use"] = 0 donor_reactions = donor_reactions or set() if pathway_logic_network.empty or not pairs: return [], [] @@ -3458,6 +3467,15 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram continue trans.add((au, bu, fu, f)) trans.add((bu, au, ru, r)) + uses: Dict[str, Set[Tuple[str, str]]] = {} + for a_u, b_u, rx, _ in trans: + uses.setdefault(rx, set()).add((a_u, b_u)) + multi = {rx for rx, s in uses.items() if len(s) > 1} + if multi: + stats["multi_use"] = len(multi) + bad = {(a_u, b_u) for a_u, b_u, rx, _ in trans if rx in multi} + bad |= {(b, a) for a, b in bad} + trans = {t for t in trans if (t[0], t[1]) not in bad} if not trans: return [], [] parent: Dict[str, str] = {} @@ -3478,15 +3496,15 @@ def find(x): def modified_side(a_u, b_u, f_stid, r_stid): """The node the pair marks as modified, or None when undecidable.""" + ma, ca = mod_counts.get(ent.get(a_u, ""), (0, 0)) + mb, cb = mod_counts.get(ent.get(b_u, ""), (0, 0)) + if ma != mb: + return b_u if mb > ma else a_u fd, rd = f_stid in donor_reactions, r_stid in donor_reactions if fd != rd: return b_u if fd else a_u - ma, ca = mod_counts.get(ent.get(a_u, ""), (0, 0)) - mb, cb = mod_counts.get(ent.get(b_u, ""), (0, 0)) if ca != cb: return b_u if cb > ca else a_u - if ma != mb: - return b_u if mb > ma else a_u return None marked: Dict[str, int] = {} fwd = {(a_u, b_u): rs for a_u, b_u, _, rs in trans} @@ -3507,6 +3525,8 @@ def base_key(u): for i, members in enumerate(order, start=1): pid = f"pool{i}" base = min(members, key=base_key) + if not any(marked.get(u, 0) for u in members): + stats["ties"] += 1 for u in sorted(members): forms.append((pid, u, ent.get(u, ""), u == base)) for a_u, b_u, rx, rs in sorted(trans): @@ -3523,15 +3543,17 @@ def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_ even when empty, so a consumer can tell "no pools" from "no table".""" from src.neo4j_connector import (get_interconversion_pairs, get_form_modification_counts, get_donor_reactions) - pairs = get_interconversion_pairs(pathway_id) + pairs = get_interconversion_pairs(pathway_id, _UBIQUITIN_STIDS) mods = get_form_modification_counts({x for p in pairs for x in (p[2], p[3])}) donors = get_donor_reactions({x for p in pairs for x in (p[0], p[1])}, _UBIQUITIN_STIDS) - forms, transitions = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, pairs, mods, donors) + stats: Dict[str, int] = {} + forms, transitions = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, pairs, mods, donors, stats) pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "is_base"]).to_csv(pools_file, index=False) pd.DataFrame(transitions, columns=["pool_id", "from_uuid", "to_uuid", "reaction_uuid", "reaction_stid"]).to_csv( transitions_file, index=False) logger.info(f"Exported {len({f[0] for f in forms})} pools ({len(forms)} forms, " - f"{len(transitions)} transitions) to {pools_file}") + f"{len(transitions)} transitions; {stats['ties']} oriented by tie-break, " + f"{stats['multi_use']} multi-use reactions dropped) to {pools_file}") def export_node_resolution(pathway_id: str, diff --git a/src/neo4j_connector.py b/src/neo4j_connector.py index b9399f1..7a03cb2 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -957,20 +957,37 @@ def get_donor_reactions(reaction_ids, ubiquitin_stids=()) -> Set[str]: return {r["rid"] for r in rows} -def get_interconversion_pairs(pathway_id: str) -> List[Tuple[str, str, str, str]]: +def get_interconversion_pairs(pathway_id: str, ubiquitin_stids=()) -> List[Tuple[str, str, str, str]]: """Curated interconversion pairs in a pathway (deltasignal specs/039): (forward reaction, reverse reaction, A, B) with F: A -> B and R: B -> A, both - events of the pathway, A and B distinct non-small-molecule entities. - Returned as stIds, sorted, each unordered pair once per orientation.""" + events of the pathway, and A, B two FORMS OF ONE PROTEIN: + + - each reaction converts one form into the other: its only non-small-molecule + input is the source form and its only non-small-molecule output is the + other form (ubiquitin counts as a co-substrate). That excludes binding + (E + S -> E:S), release, degradation, and reactions taking two forms; + - the two forms share a protein (a reference entity of an EWAS leaf). + + Review of specs/039: without these, 78 of 102 shipped pools were enzyme + binding cycles, non-shared-protein pairs or degradation "reverse" steps. + Returned as stIds, sorted, each pair once per orientation.""" query = """ MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(r1:ReactionLikeEvent)-[:input]->(a:PhysicalEntity), (r1)-[:output]->(b:PhysicalEntity) MATCH (p)-[:hasEvent*]->(r2:ReactionLikeEvent)-[:input]->(b), (r2)-[:output]->(a) WHERE r1 <> r2 AND a <> b AND NOT a:SimpleEntity AND NOT b:SimpleEntity + AND size([(r1)-[:input]->(x) WHERE x <> a AND NOT x:SimpleEntity AND NOT x.stId IN $ub | x]) = 0 + AND size([(r1)-[:output]->(y) WHERE y <> b AND NOT y:SimpleEntity AND NOT y.stId IN $ub | y]) = 0 + AND size([(r2)-[:input]->(x) WHERE x <> b AND NOT x:SimpleEntity AND NOT x.stId IN $ub | x]) = 0 + AND size([(r2)-[:output]->(y) WHERE y <> a AND NOT y:SimpleEntity AND NOT y.stId IN $ub | y]) = 0 + WITH DISTINCT r1, r2, a, b, + [(a)-[:hasComponent|hasMember|hasCandidate*0..6]->(:EntityWithAccessionedSequence)-[:referenceEntity]->(re) | re.dbId] AS ra, + [(b)-[:hasComponent|hasMember|hasCandidate*0..6]->(:EntityWithAccessionedSequence)-[:referenceEntity]->(re) | re.dbId] AS rb + WHERE any(x IN ra WHERE x IN rb) RETURN DISTINCT r1.stId AS f, r2.stId AS r, a.stId AS a, b.stId AS b """ try: - rows = get_graph().run(query, pid=pathway_id).data() + rows = get_graph().run(query, pid=pathway_id, ub=list(ubiquitin_stids)).data() except Exception: logger.error("Error in get_interconversion_pairs", **_traceback_kwargs()) raise diff --git a/tests/test_pools.py b/tests/test_pools.py index 5a92358..d383f80 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -13,11 +13,14 @@ PAIRS = [(F, R, A, B), (R, F, B, A)] +def e(s, t, et): + return {"source_id": U[s], "target_id": U[t], "edge_type": et} + + def net(loop=True): # a -> f -> b -> r -> a (the node-level loop); with loop=False, R's output is a # different copy of A (a2), which the uuids keep apart a_out = "a" if loop else "a2" - e = lambda s, t, et: {"source_id": U[s], "target_id": U[t], "edge_type": et} edges = pd.DataFrame([e("a", "vf", "input"), e("vf", "b", "output"), e("b", "vr", "input"), e("vr", a_out, "output")]) rmap = pd.DataFrame({"uid": [U["vf"], U["vr"]], "reactome_id": [F, R]}) @@ -39,17 +42,36 @@ def test_copies_the_uuids_separated_are_not_pooled(): assert forms == [] and trans == [] -def test_orientation_by_donor_then_components_then_residues(): +def test_orientation_by_residues_then_donor_then_components(): edges, rmap, umap = net() - # the donor sits on the REVERSE reaction: then GDP is the modified side + # residues decide first, whatever the donor says + forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (2, 1), B: (0, 1)}, {F}) + assert dict((u, b) for _, u, _, b in forms)[U["b"]] is True + # no residue difference: the donor sits on the REVERSE reaction, so GDP is modified forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {R}) assert dict((u, b) for _, u, _, b in forms) == {U["a"]: False, U["b"]: True} - # no donor: the form with more components is the modified one (a complex) + # no residue or donor difference: the form with more components is modified forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 3)}, set()) assert dict((u, b) for _, u, _, b in forms)[U["a"]] is True - # no donor, same components: more modified residues is modified - forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (2, 1), B: (0, 1)}, set()) - assert dict((u, b) for _, u, _, b in forms)[U["b"]] is True + + +def test_undecided_pool_is_counted_as_a_tie(): + edges, rmap, umap = net() + st = {} + forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 1)}, set(), st) + assert st["ties"] == 1 and len(forms) == 2 + st = {} + m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 1)}, {F}, st) + assert st["ties"] == 0 + + +def test_a_reaction_converting_two_pairs_is_dropped(): + # vf also converts a2 -> b: its node would be written by two fluxes + edges, rmap, umap = net() + edges = pd.concat([edges, pd.DataFrame([e("a2", "vf", "input"), e("vr", "a2", "output")])], ignore_index=True) + st = {} + forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}, st) + assert forms == [] and trans == [] and st["multi_use"] >= 1 def test_no_pairs_or_empty_network_give_no_pools(): From f4282062bd592e80dd8f049ede76c85b30ab08f6 Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Sun, 27 Sep 2026 23:55:06 -0400 Subject: [PATCH 3/7] pools: export whether each transition is catalysed (specs/039 amendment 2) A consumer weighs an intrinsic (uncatalysed) step beside a catalysed one for the same pair; the equal split made SOS1 KO read RAS:GTP 0.53 (Fable review). Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 12 +++++++++--- tests/test_pools.py | 10 +++++++++- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index cb5deac..227bc04 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3422,7 +3422,10 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram Returns (forms, transitions): forms: [(pool_id, node_uuid, stable_id, is_base)] - transitions: [(pool_id, from_uuid, to_uuid, reaction_uuid, reaction_stid)] + transitions: [(pool_id, from_uuid, to_uuid, reaction_uuid, reaction_stid, catalysed)] + ``catalysed`` says whether the reaction node has a catalyst edge; a consumer + weighs an uncatalysed (intrinsic) step beside a catalysed one for the same + pair (specs/039 amendment 2). A reaction node that would be a transition more than once (it converts two different form pairs) is dropped, with every loop through it: its node would be written by two fluxes (review of specs/039, finding 2). @@ -3446,8 +3449,11 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram vr = {str(u): str(r) for u, r in zip(reaction_id_map["uid"], reaction_id_map["reactome_id"])} ins: Dict[str, Set[str]] = {} outs: Dict[str, Set[str]] = {} + catalysed: Set[str] = set() for _, e in pathway_logic_network.iterrows(): s, t, et = str(e["source_id"]), str(e["target_id"]), e.get("edge_type") + if et == "catalyst" and t in vr: + catalysed.add(t) if et == "input" and t in vr: ins.setdefault(t, set()).add(s) elif et == "output" and s in vr: @@ -3531,7 +3537,7 @@ def base_key(u): forms.append((pid, u, ent.get(u, ""), u == base)) for a_u, b_u, rx, rs in sorted(trans): if a_u in members: - transitions.append((pid, a_u, b_u, rx, rs)) + transitions.append((pid, a_u, b_u, rx, rs, rx in catalysed)) return forms, transitions @@ -3549,7 +3555,7 @@ def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_ stats: Dict[str, int] = {} forms, transitions = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, pairs, mods, donors, stats) pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "is_base"]).to_csv(pools_file, index=False) - pd.DataFrame(transitions, columns=["pool_id", "from_uuid", "to_uuid", "reaction_uuid", "reaction_stid"]).to_csv( + pd.DataFrame(transitions, columns=["pool_id", "from_uuid", "to_uuid", "reaction_uuid", "reaction_stid", "catalysed"]).to_csv( transitions_file, index=False) logger.info(f"Exported {len({f[0] for f in forms})} pools ({len(forms)} forms, " f"{len(transitions)} transitions; {stats['ties']} oriented by tie-break, " diff --git a/tests/test_pools.py b/tests/test_pools.py index d383f80..36947d9 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -32,7 +32,7 @@ def test_a_node_level_loop_is_one_pool_with_the_resting_form_as_base(): edges, rmap, umap = net() forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {F}) assert sorted((u, base) for _, u, _, base in forms) == sorted([(U["a"], True), (U["b"], False)]) - assert sorted((f, t, rx) for _, f, t, rx, _ in trans) == sorted([(U["a"], U["b"], U["vf"]), (U["b"], U["a"], U["vr"])]) + assert sorted((f, t, rx) for _, f, t, rx, _, _ in trans) == sorted([(U["a"], U["b"], U["vf"]), (U["b"], U["a"], U["vr"])]) assert {p for p, *_ in forms} == {"pool1"} @@ -78,3 +78,11 @@ def test_no_pairs_or_empty_network_give_no_pools(): edges, rmap, umap = net() assert m.find_pools(edges, rmap, umap, [], {}, set()) == ([], []) assert m.find_pools(pd.DataFrame(columns=["source_id", "target_id", "edge_type"]), rmap, umap, PAIRS, {}, set()) == ([], []) + + +def test_transitions_say_whether_they_are_catalysed(): + edges, rmap, umap = net() + edges = pd.concat([edges, pd.DataFrame([{"source_id": U["a2"], "target_id": U["vf"], "edge_type": "catalyst"}])], + ignore_index=True) + _, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}) + assert {rx: c for _, _, _, rx, _, c in trans} == {U["vf"]: True, U["vr"]: False} From fa5921f7c824a9f89e8087867aed4bb9be1f507f Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Mon, 28 Sep 2026 00:32:38 -0400 Subject: [PATCH 4/7] pools: multi-step modification cycles (deltasignal specs/039 amendment 2) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detection now follows one reference entity R through every R-step (a reaction whose exactly one input and exactly one output contain R, both at stoichiometry 1, and differ) instead of one-reaction A <-> B pairs, so a cycle curated as bind / modify / release (RAS:GTP + GAP -> RAS:GTP:GAP -> RAS:GDP + GAP; S -> S:E -> S*:E -> S* -> S*:P -> S:P -> S) is one pool. Neo4j: get_reaction_participants (inputs, outputs, catalysts per reaction), get_form_profiles (proteins, modification signature per R = residues on R's leaves + small molecules of the complex R's leaf sits in through sets, slots, components, residues; cached across pathways), get_set_descendants. r_steps computes the steps and whether each is enzyme-driven (a catalyst or joining input containing a protein and not R; a catalyst that is a form of R is self-catalysis, non-enzyme). get_interconversion_pairs and get_form_modification_counts are replaced. Node level (find_pools): the R-graph is built from reaction nodes whose R-containing input and output nodes map to the step's forms; per R, a strongly connected set with >= 2 forms and >= 2 signatures is a pool (one signature is a carrier loop); states are the fewest-slot forms per signature, the rest intermediates; transitions are simple paths state -> intermediates* -> another state of <= 6 steps, longer dropped and counted; base state by residues, donor, components (ties counted); pools identical for co-travelling proteins are one pool, a node claimed by two different pools is removed from both and everything recomputed; carriers are the enzyme free-form nodes with the release reaction nodes that produce them. A pool with more than 20,000 uuid paths (the copies per step multiply through a shared intermediate) is dropped and counted. Files: pools.csv (pool_id,node_uuid,stable_id,role,is_base), pool_transitions.csv (pool_id,path_id,step,source_uuid,target_uuid, reaction_uuid,reaction_stid,enzyme_driven — the step's flag; a path is enzyme-driven if any step is), pool_carriers.csv (pool_id,carrier_uuid, release_reaction_uuid); all written even when empty. The log reports pools, states, intermediates, paths, multi-step pools, dropped long paths, shared nodes, ties, carriers, carrier loops, merged proteins and autocatalysis by source form / product form / other form. tests/test_pools.py: 12 tests (RAS shape with GEF, self-catalysed hydrolysis and GAP bind-release; six-form ring; an enzyme's own loop is not a pool; shared-node drop; co-travelling merge; path cap; orientation and ties; set member mapping; stoichiometry). Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 481 ++++++++++++++++++++++++--------- src/neo4j_connector.py | 162 ++++++++--- src/pathway_generator.py | 1 + tests/test_pools.py | 306 ++++++++++++++++----- 4 files changed, 721 insertions(+), 229 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 227bc04..87f3104 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3409,157 +3409,382 @@ def parts(node_id: str): logger.info(f"Exported {len(rows)} drug-derived entities to {output_file}") +POOL_MAX_STEPS = 6 # a state -> intermediates -> state path longer than this is dropped +POOL_PATH_BUDGET = 20000 # a pool with more uuid paths than this is dropped (a safety cap, counted) +POOL_STATS_KEYS = ("pools", "states", "intermediates", "paths", "multi_step_pools", "long_paths_dropped", + "shared_nodes_dropped", "ties", "carriers", "autocat_source", "autocat_product", + "autocat_other", "carrier_loops", "no_path_pools", "off_path_intermediates", + "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins") + + def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, - uuid_mapping: Dict[str, str], pairs, mod_counts: Dict[str, Tuple[int, int]], - donor_reactions: Optional[Set[str]] = None, stats: Optional[Dict[str, int]] = None): - """Interconversion pools at NODE level (deltasignal specs/039). - - A curated pair (F: A -> B, R: B -> A) is a pool transition only where the - same A node feeds a copy of F and receives a copy of R's output, and the - same B node is produced by that F and consumed by that R: the loop - A -> F -> B -> R -> A exists in this network. Occurrences the per-position - uuids separated are not re-merged. Forms joined by such loops form one pool. - - Returns (forms, transitions): - forms: [(pool_id, node_uuid, stable_id, is_base)] - transitions: [(pool_id, from_uuid, to_uuid, reaction_uuid, reaction_stid, catalysed)] - ``catalysed`` says whether the reaction node has a catalyst edge; a consumer - weighs an uncatalysed (intrinsic) step beside a catalysed one for the same - pair (specs/039 amendment 2). - A reaction node that would be a transition more than once (it converts two - different form pairs) is dropped, with every loop through it: its node would - be written by two fluxes (review of specs/039, finding 2). - - The base form (the unmodified, resting form) is decided per pair, then per pool - (specs/039 amendment 1): the modified form is the one with more modified - residues (as pre-registered); failing that, the product of the direction that - consumes a group donor (ATP, GTP, SAM, acetyl-CoA, NAD+, ubiquitin: RAS:GDP - and RAS:GTP carry no residue difference); failing that, the form with more - components. The pool's base is the form no pair marks as modified; a pool no - pair decides falls back to the smaller stId and is counted in - ``stats["ties"]``, with the dropped reactions in ``stats["multi_use"]``. + uuid_mapping: Dict[str, str], steps, profiles: Dict[str, Dict[str, Any]], + donor_reactions: Optional[Set[str]] = None, + set_members: Optional[Dict[str, Set[str]]] = None, + stats: Optional[Dict[str, int]] = None): + """Modification-cycle pools at NODE level (deltasignal specs/039 amendment 2). + + ``steps`` are the pathway's R-steps at stId level, ``(reaction, R, input + form, output form, enzyme_driven)``: a reaction whose exactly one input and + exactly one output contain the reference entity R, both at stoichiometry 1, + and which differ (see :func:`r_steps`). ``profiles`` is + :func:`get_form_profiles` over every node's entity (and a variant's chosen + members); ``set_members`` says which nodes a set-valued participant became. + + A step exists at node level only as a reaction node whose R-containing + input node and output node map to the step's forms (the node's entity, its + variant parent, or a member of the step's set). Per R, the R-graph over + form nodes is condensed; a strongly connected set with >= 2 forms and >= 2 + modification signatures is a pool candidate (one signature is a carrier + loop: an enzyme's E -> E:S -> E). In each signature the forms with the + fewest non-small slots are states, the rest intermediates. Transitions are + simple directed paths state -> intermediates* -> another state of at most + ``POOL_MAX_STEPS`` steps; each uuid path is one parallel transition, and + longer ones are dropped and counted. Because an intermediate is one node + however many reaction copies enter and leave it, the uuid paths through it + are the product of the copies per step (RAF: 48 GAP binding copies x 48 + release copies = 2,304 two-step paths); a pool with more than + ``POOL_PATH_BUDGET`` paths is dropped and counted (``budget_dropped``). The base state follows amendment 1 + (residues, then the donor-consuming direction, then components; undecided + pools fall back to the smaller stId and are counted). Pools identical for + several R (the members of one set travelling together) are one pool; a node + claimed by two DIFFERENT pools is removed from both, everything is + recomputed, and the node is counted. Carriers are the enzyme's free-form + nodes (a joining input of an enzyme-driven step, containing a protein and + not R, that a step of the pool releases) with the reaction nodes releasing + them. + + Returns (forms, transitions, carriers): + forms: [(pool_id, node_uuid, stable_id, role, is_base)], role state | intermediate + transitions: [(pool_id, path_id, step, source_uuid, target_uuid, reaction_uuid, reaction_stid, enzyme_driven)] + ``enzyme_driven`` is the STEP's flag; a path is enzyme-driven if any step is. + carriers: [(pool_id, carrier_uuid, release_reaction_uuid)] + ``stats`` receives the counts in ``POOL_STATS_KEYS``; the autocatalysis + counts are distinct (reaction node, catalyst) pairs where a form of the pool + catalyses a step of its own pool, split by whether the catalyst is the + step's source form, its product form, or another form. """ + import networkx as nx stats = stats if stats is not None else {} - stats["ties"] = 0 - stats["multi_use"] = 0 + for k in POOL_STATS_KEYS: + stats[k] = 0 donor_reactions = donor_reactions or set() - if pathway_logic_network.empty or not pairs: - return [], [] + set_members = set_members or {} + steps = list(steps) + if pathway_logic_network.empty or not steps: + return [], [], [] ent = _uuid_to_stable_id_map(pathway_logic_network, uuid_mapping) vr = {str(u): str(r) for u, r in zip(reaction_id_map["uid"], reaction_id_map["reactome_id"])} - ins: Dict[str, Set[str]] = {} - outs: Dict[str, Set[str]] = {} - catalysed: Set[str] = set() + ins: Dict[str, Dict[str, float]] = {} + outs: Dict[str, Dict[str, float]] = {} + cats: Dict[str, Set[str]] = {} + has_st = "stoichiometry" in pathway_logic_network.columns for _, e in pathway_logic_network.iterrows(): s, t, et = str(e["source_id"]), str(e["target_id"]), e.get("edge_type") - if et == "catalyst" and t in vr: - catalysed.add(t) + st = 1.0 + if has_st: + try: + st = float(e["stoichiometry"]) + except (TypeError, ValueError): + st = 1.0 + if st != st: # NaN + st = 1.0 if et == "input" and t in vr: - ins.setdefault(t, set()).add(s) + ins.setdefault(t, {})[s] = st elif et == "output" and s in vr: - outs.setdefault(s, set()).add(t) + outs.setdefault(s, {})[t] = st + elif et == "catalyst" and t in vr: + cats.setdefault(t, set()).add(s) copies: Dict[str, List[str]] = {} for u, r in vr.items(): copies.setdefault(r, []).append(u) - trans = set() - for f, r, a, b in pairs: - for fu in sorted(copies.get(f, [])): - for ru in sorted(copies.get(r, [])): - for au in sorted(ins.get(fu, set()) & outs.get(ru, set())): - if ent.get(au) != a: - continue - for bu in sorted(outs.get(fu, set()) & ins.get(ru, set())): - if ent.get(bu) != b: + + def parts(u: str) -> Tuple[str, List[str]]: + s = ent.get(u, "") + if "::variant::" not in s: + return s, [] + parent, members = s.split("::variant::", 1) + return parent, [m for m in members.split("_") if m] + + def pe_of(u: str) -> str: + return parts(u)[0] + + prot_cache: Dict[str, frozenset] = {} + + def proteins(u: str) -> frozenset: + if u not in prot_cache: + parent, members = parts(u) + if not members: + out = set(profiles.get(parent, {}).get("proteins", ())) + else: + out = set(profiles.get(parent, {}).get("fixed", ())) + for m in members: + out |= set(profiles.get(m, {}).get("proteins", ())) + prot_cache[u] = frozenset(out) + return prot_cache[u] + + def prof(u: str, key: str, default=0): + return profiles.get(pe_of(u), {}).get(key, default) + + def sig(u: str, r): + return profiles.get(pe_of(u), {}).get("sig", {}).get(r, frozenset()) + + def maps(u: str, pe: str) -> bool: + p = pe_of(u) + return p == pe or p in set_members.get(pe, ()) + + # uuid-level steps: (R, in node, out node, reaction node) -> (reaction stId, enzyme_driven) + ustep: Dict[Tuple[Any, str, str, str], Tuple[str, bool]] = {} + for rx, r, a, b, enz in steps: + for rx_u in sorted(copies.get(rx, [])): + cin = [u for u, st in ins.get(rx_u, {}).items() if st == 1 and maps(u, a) and r in proteins(u)] + cout = [u for u, st in outs.get(rx_u, {}).items() if st == 1 and maps(u, b) and r in proteins(u)] + if len(cin) != 1 or len(cout) != 1: + if len(cin) > 1 or len(cout) > 1: + stats["ambiguous_copies"] += 1 + continue + if cin[0] != cout[0]: + ustep[(r, cin[0], cout[0], rx_u)] = (rx, bool(enz)) + if not ustep: + return [], [], [] + + def build(removed: Set[str]): + local = {k: 0 for k in POOL_STATS_KEYS} + by_r: Dict[Any, List[Tuple[str, str, str, str, bool]]] = {} + for (r, a_u, b_u, rx_u), (rs, enz) in ustep.items(): + if a_u in removed or b_u in removed or rx_u in removed: + continue + by_r.setdefault(r, []).append((a_u, b_u, rx_u, rs, enz)) + pools: Dict[tuple, Dict[str, Any]] = {} + for r in sorted(by_r): + g = nx.DiGraph() + g.add_edges_from((a_u, b_u) for a_u, b_u, _, _, _ in by_r[r]) + for scc in nx.strongly_connected_components(g): + if len(scc) < 2: + continue + by_sig: Dict[frozenset, List[str]] = {} + for u in scc: + by_sig.setdefault(sig(u, r), []).append(u) + if len(by_sig) < 2: + local["carrier_loops"] += 1 + continue + states: Set[str] = set() + for members in by_sig.values(): + fewest = min(prof(u, "slots") for u in members) + states |= {u for u in members if prof(u, "slots") == fewest} + adj: Dict[str, List[Tuple[str, str, str, bool]]] = {} + for a_u, b_u, rx_u, rs, enz in sorted(by_r[r]): + if a_u in scc and b_u in scc: + adj.setdefault(a_u, []).append((b_u, rx_u, rs, enz)) + paths: List[Tuple[Tuple[str, str, str, str, bool], ...]] = [] + over_budget = False + + def walk(u, so_far, seen): + nonlocal over_budget + for v, rx_u, rs, enz in adj.get(u, ()): + if over_budget: + return + if v in seen: + continue + new = so_far + ((u, v, rx_u, rs, enz),) + if v in states: + paths.append(new) + if len(paths) > POOL_PATH_BUDGET: + over_budget = True + return continue - trans.add((au, bu, fu, f)) - trans.add((bu, au, ru, r)) - uses: Dict[str, Set[Tuple[str, str]]] = {} - for a_u, b_u, rx, _ in trans: - uses.setdefault(rx, set()).add((a_u, b_u)) - multi = {rx for rx, s in uses.items() if len(s) > 1} - if multi: - stats["multi_use"] = len(multi) - bad = {(a_u, b_u) for a_u, b_u, rx, _ in trans if rx in multi} - bad |= {(b, a) for a, b in bad} - trans = {t for t in trans if (t[0], t[1]) not in bad} - if not trans: - return [], [] - parent: Dict[str, str] = {} - - def find(x): - parent.setdefault(x, x) - while parent[x] != x: - parent[x] = parent[parent[x]] - x = parent[x] - return x - for a_u, b_u, _, _ in trans: - ra, rb = find(a_u), find(b_u) - if ra != rb: - parent[max(ra, rb)] = min(ra, rb) - groups: Dict[str, List[str]] = {} - for n in sorted(parent): - groups.setdefault(find(n), []).append(n) - - def modified_side(a_u, b_u, f_stid, r_stid): - """The node the pair marks as modified, or None when undecidable.""" - ma, ca = mod_counts.get(ent.get(a_u, ""), (0, 0)) - mb, cb = mod_counts.get(ent.get(b_u, ""), (0, 0)) - if ma != mb: - return b_u if mb > ma else a_u - fd, rd = f_stid in donor_reactions, r_stid in donor_reactions + if len(new) >= POOL_MAX_STEPS: + local["long_paths_dropped"] += 1 + continue + walk(v, new, seen | {v}) + for s0 in sorted(states): + walk(s0, (), {s0}) + if over_budget: + local["budget_dropped"] += 1 + continue + if not paths: + local["no_path_pools"] += 1 + continue + on_path: Set[str] = {n for p in paths for st in p for n in (st[0], st[1])} + kept_states = states & on_path + if len({sig(u, r) for u in kept_states}) < 2: + local["no_path_pools"] += 1 + continue + inter = (scc - states) & on_path + local["off_path_states"] += len(states - on_path) + local["off_path_intermediates"] += len(scc - states - on_path) + key = (frozenset(kept_states), frozenset(inter), + frozenset(tuple((a_u, b_u, rx_u) for a_u, b_u, rx_u, _, _ in p) for p in paths)) + if key in pools: + pools[key]["r"].add(r) + local["merged_proteins"] += 1 + else: + pools[key] = {"r": {r}, "states": kept_states, "inter": inter, "paths": paths, + "nodes": kept_states | inter | {st[2] for p in paths for st in p}} + return list(pools.values()), local + + removed: Set[str] = set() + dropped_shared = 0 + while True: + pools, local = build(removed) + claimed: Dict[str, int] = {} + for p in pools: + for n in p["nodes"]: + claimed[n] = claimed.get(n, 0) + 1 + shared = {n for n, c in claimed.items() if c > 1} + if not shared: + break + removed |= shared + dropped_shared += len(shared) + stats.update(local) + stats["shared_nodes_dropped"] = dropped_shared + + def path_key(p): + return tuple((ent.get(a, ""), a, ent.get(b, ""), b, rx) for a, b, rx, _, _ in p) + + def modified_side(x, y, xy, yx): + """The state a pair marks as modified, or None when undecidable.""" + mx, my = prof(x, "mods"), prof(y, "mods") + if mx != my: + return y if my > mx else x + fd = any(rs in donor_reactions for p in xy for _, _, _, rs, _ in p) + rd = any(rs in donor_reactions for p in yx for _, _, _, rs, _ in p) if fd != rd: - return b_u if fd else a_u - if ca != cb: - return b_u if cb > ca else a_u + return y if fd else x + cx, cy = prof(x, "comps"), prof(y, "comps") + if cx != cy: + return y if cy > cx else x return None - marked: Dict[str, int] = {} - fwd = {(a_u, b_u): rs for a_u, b_u, _, rs in trans} - for (a_u, b_u), fs in fwd.items(): - rs = fwd.get((b_u, a_u)) - if rs is None or a_u > b_u: - continue - side = modified_side(a_u, b_u, fs, rs) - if side is not None: - marked[side] = marked.get(side, 0) + 1 - def base_key(u): - s = ent.get(u, "") - mods, comps = mod_counts.get(s, (0, 0)) - return (marked.get(u, 0), mods, comps, s, u) - forms, transitions = [], [] - order = sorted(groups.values(), key=lambda g: min(ent.get(u, "") + u for u in g)) - for i, members in enumerate(order, start=1): + forms, transitions, carriers = [], [], [] + pools.sort(key=lambda p: min(ent.get(u, "") + u for u in p["states"])) + for i, p in enumerate(pools, start=1): pid = f"pool{i}" - base = min(members, key=base_key) - if not any(marked.get(u, 0) for u in members): + states, inter, paths, rset = p["states"], p["inter"], sorted(p["paths"], key=path_key), p["r"] + between: Dict[Tuple[str, str], List[tuple]] = {} + for path in paths: + between.setdefault((path[0][0], path[-1][1]), []).append(path) + marked: Dict[str, int] = {} + for x in sorted(states): + for y in sorted(states): + if x >= y or not (between.get((x, y)) or between.get((y, x))): + continue + side = modified_side(x, y, between.get((x, y), []), between.get((y, x), [])) + if side is not None: + marked[side] = marked.get(side, 0) + 1 + base = min(states, key=lambda u: (marked.get(u, 0), prof(u, "mods"), prof(u, "comps"), ent.get(u, ""), u)) + if not marked: stats["ties"] += 1 - for u in sorted(members): - forms.append((pid, u, ent.get(u, ""), u == base)) - for a_u, b_u, rx, rs in sorted(trans): - if a_u in members: - transitions.append((pid, a_u, b_u, rx, rs, rx in catalysed)) - return forms, transitions + for u in sorted(states | inter, key=lambda u: (ent.get(u, ""), u)): + forms.append((pid, u, ent.get(u, ""), "state" if u in states else "intermediate", u == base)) + for n, path in enumerate(paths, start=1): + for k, (a_u, b_u, rx_u, rs, enz) in enumerate(path, start=1): + transitions.append((pid, f"{pid}_p{n}", k, a_u, b_u, rx_u, rs, enz)) + rx_on_paths = sorted({st[2] for path in paths for st in path}) + found: Set[Tuple[str, str]] = set() + autocat: Set[Tuple[str, str, str]] = set() + for path in paths: + for a_u, b_u, rx_u, _, enz in path: + for c in cats.get(rx_u, ()): + if c in p["nodes"]: # a form of the pool catalyses a step of its own pool + autocat.add((rx_u, c, "autocat_source" if c == a_u else "autocat_product" if c == b_u + else "autocat_other")) + if not enz: + continue + for j in ins.get(rx_u, {}): + if j == a_u or j in p["nodes"] or pe_of(j) in _UBIQUITIN_STIDS: + continue + pj = proteins(j) + if not pj or pj & rset: + continue + for rel in rx_on_paths: + if j in outs.get(rel, {}): + found.add((j, rel)) + carriers.extend((pid, j, rel) for j, rel in sorted(found, key=lambda t: (ent.get(t[0], ""), t[0], t[1]))) + for _, _, kind in autocat: + stats[kind] += 1 + stats["pools"] += 1 + stats["states"] += len(states) + stats["intermediates"] += len(inter) + stats["paths"] += len(paths) + stats["multi_step_pools"] += int(any(len(path) > 1 for path in paths)) + stats["carriers"] += len({j for j, _ in found}) + return forms, transitions, carriers + + +def r_steps(participants, profiles: Dict[str, Dict[str, Any]], ubiquitin_stids=()) -> List[Tuple[str, Any, str, str, bool]]: + """The R-steps of a pathway at stId level (deltasignal specs/039 amendment + 2), from :func:`get_reaction_participants` rows: for each reaction and each + reference entity R, the reaction is a step if exactly one input and exactly + one output contain R, both at stoichiometry 1, and they differ. A step is + enzyme-driven if a catalyst or a joining input (another input; ubiquitin + is a co-substrate) contains a protein and does not contain R: a catalyst + that is itself a form of R (RAS:GTP's intrinsic hydrolysis, curated with + the set p21 RAS:GTP as its own catalyst) is self-catalysis, which counts as + non-enzyme. Returns sorted ``(reaction, R, input, output, enzyme_driven)``.""" + by_rx: Dict[str, Dict[str, List[Tuple[str, int]]]] = {} + for rx, role, pe, st in participants: + by_rx.setdefault(rx, {}).setdefault(role, []).append((pe, st)) + out = [] + for rx in sorted(by_rx): + roles = by_rx[rx] + rs: Set[Any] = set() + for pe, _ in roles.get("input", []): + rs |= set(profiles.get(pe, {}).get("proteins", ())) + for r in sorted(rs): + a = [(pe, st) for pe, st in roles.get("input", []) if r in profiles.get(pe, {}).get("proteins", ())] + b = [(pe, st) for pe, st in roles.get("output", []) if r in profiles.get(pe, {}).get("proteins", ())] + if len(a) != 1 or len(b) != 1 or a[0][1] != 1 or b[0][1] != 1 or a[0][0] == b[0][0]: + continue + helpers = [pe for pe, _ in roles.get("catalyst", [])] + helpers += [pe for pe, _ in roles.get("input", []) if pe != a[0][0] and pe not in ubiquitin_stids] + enz = any(profiles.get(pe, {}).get("proteins") and r not in profiles[pe]["proteins"] for pe in helpers) + out.append((rx, r, a[0][0], b[0][0], enz)) + return out def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, - uuid_mapping: Dict[str, str], pools_file: str, transitions_file: str) -> None: - """Write pools.csv and pool_transitions.csv (deltasignal specs/039): each - protein's interconversion pool, for a consumer that solves it at steady state - (pi = pi P) instead of multiplying around the loop. Both files are written - even when empty, so a consumer can tell "no pools" from "no table".""" - from src.neo4j_connector import (get_interconversion_pairs, get_form_modification_counts, - get_donor_reactions) - pairs = get_interconversion_pairs(pathway_id, _UBIQUITIN_STIDS) - mods = get_form_modification_counts({x for p in pairs for x in (p[2], p[3])}) - donors = get_donor_reactions({x for p in pairs for x in (p[0], p[1])}, _UBIQUITIN_STIDS) + uuid_mapping: Dict[str, str], pools_file: str, transitions_file: str, + carriers_file: str) -> None: + """Write pools.csv, pool_transitions.csv and pool_carriers.csv (deltasignal + specs/039 amendment 2): each protein's modification-cycle pool, for a + consumer that solves it at steady state (pi = pi P) instead of multiplying + around the loop. All three files are written even when empty, so a consumer + can tell "no pools" from "no table".""" + from src.neo4j_connector import (get_reaction_participants, get_form_profiles, + get_donor_reactions, get_set_descendants) + participants = get_reaction_participants(pathway_id) + ent = _uuid_to_stable_id_map(pathway_logic_network, uuid_mapping) + wanted = {pe for _, _, pe, _ in participants} + for s in ent.values(): + if "::variant::" in s: + parent, members = s.split("::variant::", 1) + wanted.add(parent) + wanted.update(m for m in members.split("_") if m) + else: + wanted.add(s) + profiles = get_form_profiles(wanted) + steps = r_steps(participants, profiles, _UBIQUITIN_STIDS) + donors = get_donor_reactions({s[0] for s in steps}, _UBIQUITIN_STIDS) + set_members = get_set_descendants({x for s in steps for x in (s[2], s[3])}) stats: Dict[str, int] = {} - forms, transitions = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, pairs, mods, donors, stats) - pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "is_base"]).to_csv(pools_file, index=False) - pd.DataFrame(transitions, columns=["pool_id", "from_uuid", "to_uuid", "reaction_uuid", "reaction_stid", "catalysed"]).to_csv( - transitions_file, index=False) - logger.info(f"Exported {len({f[0] for f in forms})} pools ({len(forms)} forms, " - f"{len(transitions)} transitions; {stats['ties']} oriented by tie-break, " - f"{stats['multi_use']} multi-use reactions dropped) to {pools_file}") + forms, transitions, carriers = find_pools(pathway_logic_network, reaction_id_map, uuid_mapping, steps, + profiles, donors, set_members, stats) + pd.DataFrame(forms, columns=["pool_id", "node_uuid", "stable_id", "role", "is_base"]).to_csv(pools_file, index=False) + pd.DataFrame(transitions, columns=["pool_id", "path_id", "step", "source_uuid", "target_uuid", "reaction_uuid", + "reaction_stid", "enzyme_driven"]).to_csv(transitions_file, index=False) + pd.DataFrame(carriers, columns=["pool_id", "carrier_uuid", "release_reaction_uuid"]).to_csv(carriers_file, index=False) + logger.info( + f"Exported {stats['pools']} pools ({stats['states']} states, {stats['intermediates']} intermediates, " + f"{stats['paths']} paths, {stats['multi_step_pools']} multi-step; {len(steps)} R-steps at stId level; " + f"{stats['long_paths_dropped']} paths over {POOL_MAX_STEPS} steps dropped, " + f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['ties']} oriented by tie-break, " + f"{stats['carriers']} carriers, {stats['carrier_loops']} carrier loops, " + f"{stats['merged_proteins']} co-travelling proteins merged, {stats['no_path_pools']} candidates with no path, " + f"{stats['off_path_states']}/{stats['off_path_intermediates']} states/intermediates off every path, " + f"{stats['ambiguous_copies']} ambiguous reaction copies, {stats['budget_dropped']} over the path budget; " + f"autocatalysis by source form {stats['autocat_source']}, by product form {stats['autocat_product']}, " + f"by another form {stats['autocat_other']}) to {pools_file}") def export_node_resolution(pathway_id: str, diff --git a/src/neo4j_connector.py b/src/neo4j_connector.py index 7a03cb2..2873847 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -957,61 +957,143 @@ def get_donor_reactions(reaction_ids, ubiquitin_stids=()) -> Set[str]: return {r["rid"] for r in rows} -def get_interconversion_pairs(pathway_id: str, ubiquitin_stids=()) -> List[Tuple[str, str, str, str]]: - """Curated interconversion pairs in a pathway (deltasignal specs/039): - (forward reaction, reverse reaction, A, B) with F: A -> B and R: B -> A, both - events of the pathway, and A, B two FORMS OF ONE PROTEIN: - - - each reaction converts one form into the other: its only non-small-molecule - input is the source form and its only non-small-molecule output is the - other form (ubiquitin counts as a co-substrate). That excludes binding - (E + S -> E:S), release, degradation, and reactions taking two forms; - - the two forms share a protein (a reference entity of an EWAS leaf). - - Review of specs/039: without these, 78 of 102 shipped pools were enzyme - binding cycles, non-shared-protein pairs or degradation "reverse" steps. - Returned as stIds, sorted, each pair once per orientation.""" - query = """ - MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(r1:ReactionLikeEvent)-[:input]->(a:PhysicalEntity), - (r1)-[:output]->(b:PhysicalEntity) - MATCH (p)-[:hasEvent*]->(r2:ReactionLikeEvent)-[:input]->(b), (r2)-[:output]->(a) - WHERE r1 <> r2 AND a <> b AND NOT a:SimpleEntity AND NOT b:SimpleEntity - AND size([(r1)-[:input]->(x) WHERE x <> a AND NOT x:SimpleEntity AND NOT x.stId IN $ub | x]) = 0 - AND size([(r1)-[:output]->(y) WHERE y <> b AND NOT y:SimpleEntity AND NOT y.stId IN $ub | y]) = 0 - AND size([(r2)-[:input]->(x) WHERE x <> b AND NOT x:SimpleEntity AND NOT x.stId IN $ub | x]) = 0 - AND size([(r2)-[:output]->(y) WHERE y <> a AND NOT y:SimpleEntity AND NOT y.stId IN $ub | y]) = 0 - WITH DISTINCT r1, r2, a, b, - [(a)-[:hasComponent|hasMember|hasCandidate*0..6]->(:EntityWithAccessionedSequence)-[:referenceEntity]->(re) | re.dbId] AS ra, - [(b)-[:hasComponent|hasMember|hasCandidate*0..6]->(:EntityWithAccessionedSequence)-[:referenceEntity]->(re) | re.dbId] AS rb - WHERE any(x IN ra WHERE x IN rb) - RETURN DISTINCT r1.stId AS f, r2.stId AS r, a.stId AS a, b.stId AS b +def get_reaction_participants(pathway_id: str) -> List[Tuple[str, str, str, int]]: + """(reaction stId, role, entity stId, stoichiometry) for every input, output + and catalyst of every reaction-like event of a pathway (deltasignal + specs/039 amendment 2). Roles are ``input``, ``output`` and ``catalyst``; + a catalyst's stoichiometry is 1. One round trip per role.""" + io_query = """ + MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(rx:ReactionLikeEvent) + WITH DISTINCT rx + MATCH (rx)-[e:input|output]->(pe:PhysicalEntity) + RETURN rx.stId AS rx, type(e) AS role, pe.stId AS pe, coalesce(e.stoichiometry, 1) AS st + """ + cat_query = """ + MATCH (p:Pathway {stId: $pid})-[:hasEvent*]->(rx:ReactionLikeEvent) + WITH DISTINCT rx + MATCH (rx)-[:catalystActivity]->(:CatalystActivity)-[:physicalEntity]->(pe:PhysicalEntity) + RETURN DISTINCT rx.stId AS rx, 'catalyst' AS role, pe.stId AS pe, 1 AS st + """ + try: + rows = get_graph().run(io_query, pid=pathway_id).data() + rows += get_graph().run(cat_query, pid=pathway_id).data() + except Exception: + logger.error("Error in get_reaction_participants", **_traceback_kwargs()) + raise + return sorted({(r["rx"], r["role"], r["pe"], int(r["st"] or 1)) for r in rows}) + + +_form_profile_cache: Dict[str, Dict[str, Any]] = {} + + +def get_form_profiles(stable_ids) -> Dict[str, Dict[str, Any]]: + """stId -> what a form is, for pool detection (deltasignal specs/039 + amendment 2). Per entity: + + - ``proteins``: reference entities (dbIds) of its EWAS leaves, through + components, members and candidates; + - ``fixed``: the subset reached through components only (no set on the + way), which every variant of the entity carries; + - ``sig``: per reference entity R, its modification signature: the + modified residues on R's leaves plus the small molecules (ChEBI id, else + name) that are direct components of the complex R's leaf sits in, looking + through sets (GTP attaches to RAS in p21 RAS:GTP = [GTP, {H,K,N}RAS], and + not to the GAP in p21 RAS:GTP:RAS GAPs, which nests it); + - ``mods``: modified residues summed over all leaves (amendment 1's base + rule), ``comps``: components at any depth, ``slots``: non-small-molecule + positions (a set inside a complex is one position; a bare protein or set + is one). + + Cached per process: the 92 catalog pathways share most entities.""" + ids = sorted({s for s in stable_ids if s}) + missing = [s for s in ids if s not in _form_profile_cache] + leaves_query = """ + UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) + OPTIONAL MATCH p = (e)-[:hasComponent|hasMember|hasCandidate*0..8]->(l:EntityWithAccessionedSequence)-[:referenceEntity]->(re) + RETURN s, re.dbId AS r, + CASE WHEN p IS NULL THEN [] ELSE + [n IN nodes(p) WHERE n:PhysicalEntity | + [n.stId, CASE WHEN n:Complex THEN 'C' WHEN n:EntitySet THEN 'S' ELSE 'L' END]] END AS chain, + CASE WHEN l IS NULL THEN [] ELSE [(l)-[:hasModifiedResidue]->(m) | m.displayName] END AS res + """ + counts_query = """ + UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) + RETURN s, + size([(e)-[:hasComponent*1..6]->(c) | c]) AS comps, + size([(e)-[:hasComponent*0..6]->(x) WHERE NOT x:Complex AND NOT x:SimpleEntity | x]) AS slots + """ + smalls_query = """ + UNWIND $ids AS s MATCH (c:Complex {stId: s})-[:hasComponent]->(sm:SimpleEntity) + OPTIONAL MATCH (sm)-[:referenceEntity]->(rm) + RETURN s, coalesce(rm.identifier, head(sm.name), sm.displayName) AS mol """ try: - rows = get_graph().run(query, pid=pathway_id, ub=list(ubiquitin_stids)).data() + for i in range(0, len(missing), 200): + chunk = missing[i:i + 200] + prof: Dict[str, Dict[str, Any]] = { + s: {"proteins": set(), "fixed": set(), "sig": {}, "mods": 0, "comps": 0, "slots": 0, + "_res": {}, "_cx": {}} for s in chunk} + leaf_rows = get_graph().run(leaves_query, ids=chunk).data() + for r in get_graph().run(counts_query, ids=chunk).data(): + prof[r["s"]]["comps"] = int(r["comps"]) + prof[r["s"]]["slots"] = int(r["slots"]) + complexes = set() + for r in leaf_rows: + if r["r"] is None: + continue + pr, rid, chain, res = prof[r["s"]], int(r["r"]), r["chain"], list(r["res"]) + pr["proteins"].add(rid) + if all(k != "S" for _, k in chain[:-1]): + pr["fixed"].add(rid) + pr["mods"] += len(res) + pr["_res"].setdefault(rid, set()).update(res) + # the complex R's leaf sits in: the nearest Complex above it with + # only set steps between + cx = None + for st_id, kind in reversed(chain[:-1]): + if kind == "C": + cx = st_id + break + if kind != "S": + break + if cx is not None: + pr["_cx"].setdefault(rid, set()).add(cx) + complexes.add(cx) + smalls: Dict[str, Set[str]] = {} + for j in range(0, len(complexes), 200): + for r in get_graph().run(smalls_query, ids=sorted(complexes)[j:j + 200]).data(): + smalls.setdefault(r["s"], set()).add(str(r["mol"])) + for s, pr in prof.items(): + for rid in pr["proteins"]: + sig = {("res", x) for x in pr["_res"].get(rid, ())} + for cx in pr["_cx"].get(rid, ()): + sig |= {("mol", x) for x in smalls.get(cx, ())} + pr["sig"][rid] = frozenset(sig) + del pr["_res"], pr["_cx"] + _form_profile_cache[s] = pr except Exception: - logger.error("Error in get_interconversion_pairs", **_traceback_kwargs()) + logger.error("Error in get_form_profiles", **_traceback_kwargs()) raise - return sorted((r["f"], r["r"], r["a"], r["b"]) for r in rows) + return {s: _form_profile_cache[s] for s in ids if s in _form_profile_cache} -def get_form_modification_counts(stable_ids) -> Dict[str, Tuple[int, int]]: - """stId -> (modified residues summed over its protein leaves, component - count), to pick a pool's least-modified (base) form (deltasignal specs/039).""" +def get_set_descendants(stable_ids) -> Dict[str, Set[str]]: + """stId -> the stIds an EntitySet reaches through members and candidates + (deltasignal specs/039 amendment 2): a set-valued participant appears in + the network as its member nodes, and this says which nodes those are.""" ids = sorted({s for s in stable_ids if s}) if not ids: return {} query = """ - UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) - RETURN s, - size([(e)-[:hasComponent|hasMember|hasCandidate*0..6]->(l:EntityWithAccessionedSequence)-[:hasModifiedResidue]->(m) | m]) AS mods, - size([(e)-[:hasComponent*1..6]->(c) | c]) AS comps + UNWIND $ids AS s MATCH (e:EntitySet {stId: s}) + RETURN s, [(e)-[:hasMember|hasCandidate*1..6]->(m) | m.stId] AS members """ try: rows = get_graph().run(query, ids=ids).data() except Exception: - logger.error("Error in get_form_modification_counts", **_traceback_kwargs()) + logger.error("Error in get_set_descendants", **_traceback_kwargs()) raise - return {r["s"]: (int(r["mods"]), int(r["comps"])) for r in rows} + return {r["s"]: {str(m) for m in r["members"] if m} for r in rows} def get_reactome_release() -> Optional[int]: diff --git a/src/pathway_generator.py b/src/pathway_generator.py index 94be301..897b6bf 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -520,6 +520,7 @@ def generate_pathway_file( result.uuid_mapping, str(pathway_output_dir / "pools.csv"), str(pathway_output_dir / "pool_transitions.csv"), + str(pathway_output_dir / "pool_carriers.csv"), ) except Exception as e: logger.error(f"Failed to write pools.csv: {e}", exc_info=True) diff --git a/tests/test_pools.py b/tests/test_pools.py index 36947d9..d78c001 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -1,88 +1,272 @@ -"""Interconversion pools (deltasignal specs/039): each protein's forms joined by -curated forward/reverse reactions, found at NODE level and oriented so the base -is the resting (unmodified) form. Shape: RAS. A = RAS:GDP, B = RAS:GTP; -F = "GEF exchange" consumes GTP (a donor), R = "GAP hydrolysis".""" +"""Modification-cycle pools (deltasignal specs/039 amendment 2): a protein's +forms joined by R-steps, found at NODE level; states are the least-bound form +of each modification signature, the enzyme complexes on the way are +intermediates, and a transition is a path state -> intermediates* -> state. +Shapes: RAS (GEF / intrinsic / GAP bind-hydrolyse-release), a six-form +kinase/phosphatase ring, an enzyme's own binding loop.""" import pandas as pd import src.logic_network_generator as m # _uuid_to_stable_id_map detects the mapping direction by uuid SHAPE -U = {k: f"00000000-0000-0000-0000-{i:012d}" for i, k in enumerate(["a", "a2", "b", "vf", "vr"], 1)} +NAMES = ["g", "t", "tx", "gap", "gef", "vgef", "vhyd", "vbind", "vrel", "a2", + "s", "se", "sxe", "sx", "sxp", "sp", "e", "p", "v1", "v2", "v3", "v4", "v5", "v6", + "x", "y", "xy", "va", "vb", "i1", "i2", "i3", "i4", "i5", "i6", "v7", "v8", "vback"] +U = {k: f"00000000-0000-0000-0000-{i:012d}" for i, k in enumerate(NAMES, 1)} +RAS, GAP_P, GEF_P, SUB, KIN, PHOS = 1, 2, 3, 10, 11, 12 # reference entities -A, B, F, R = "R-HSA-GDP", "R-HSA-GTP", "R-HSA-F", "R-HSA-R" -PAIRS = [(F, R, A, B), (R, F, B, A)] +def e(s, t, et, st=1): + return {"source_id": U[s], "target_id": U[t], "edge_type": et, "stoichiometry": st} -def e(s, t, et): - return {"source_id": U[s], "target_id": U[t], "edge_type": et} +def prof(proteins, sig=None, slots=1, mods=0, comps=0, fixed=None): + return {"proteins": set(proteins), "fixed": set(fixed if fixed is not None else proteins), + "sig": {r: frozenset(sig or ()) for r in proteins}, "mods": mods, "comps": comps, "slots": slots} -def net(loop=True): - # a -> f -> b -> r -> a (the node-level loop); with loop=False, R's output is a - # different copy of A (a2), which the uuids keep apart - a_out = "a" if loop else "a2" - edges = pd.DataFrame([e("a", "vf", "input"), e("vf", "b", "output"), - e("b", "vr", "input"), e("vr", a_out, "output")]) - rmap = pd.DataFrame({"uid": [U["vf"], U["vr"]], "reactome_id": [F, R]}) - umap = {U["a"]: A, U["a2"]: A, U["b"]: B} - return edges, rmap, umap +def network(edges, reactions, umap): + return (pd.DataFrame(edges), + pd.DataFrame({"uid": [U[k] for k in reactions], "reactome_id": [f"R-HSA-{k}" for k in reactions]}), + {U[k]: s for k, s in umap.items()}) -def test_a_node_level_loop_is_one_pool_with_the_resting_form_as_base(): - edges, rmap, umap = net() - forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {F}) - assert sorted((u, base) for _, u, _, base in forms) == sorted([(U["a"], True), (U["b"], False)]) - assert sorted((f, t, rx) for _, f, t, rx, _, _ in trans) == sorted([(U["a"], U["b"], U["vf"]), (U["b"], U["a"], U["vr"])]) - assert {p for p, *_ in forms} == {"pool1"} +# --- RAS: GDP <-> GTP by GEF (enzyme), intrinsic exchange / hydrolysis (self), GAP bind + release ---------- -def test_copies_the_uuids_separated_are_not_pooled(): - edges, rmap, umap = net(loop=False) - forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}) - assert forms == [] and trans == [] +def ras_network(hydrolysis_catalyst="t"): + edges = [e("g", "vgef", "input"), e("vgef", "t", "output"), e("gef", "vgef", "catalyst"), + e("t", "vhyd", "input"), e("vhyd", "g", "output"), e(hydrolysis_catalyst, "vhyd", "catalyst"), + e("t", "vbind", "input"), e("gap", "vbind", "input"), e("vbind", "tx", "output"), + e("tx", "vrel", "input"), e("vrel", "g", "output"), e("vrel", "gap", "output")] + net, rmap, umap = network(edges, ["vgef", "vhyd", "vbind", "vrel"], + {"g": "R-HSA-G", "t": "R-HSA-T", "tx": "R-HSA-TX", "gap": "R-HSA-GAP", "gef": "R-HSA-GEF"}) + profiles = {"R-HSA-G": prof([RAS], [("mol", "GDP")], slots=1, comps=2), + "R-HSA-T": prof([RAS], [("mol", "GTP")], slots=1, comps=2), + "R-HSA-TX": {"proteins": {RAS, GAP_P}, "fixed": {RAS, GAP_P}, + "sig": {RAS: frozenset([("mol", "GTP")]), GAP_P: frozenset()}, "mods": 0, "comps": 3, "slots": 2}, + "R-HSA-GAP": prof([GAP_P]), "R-HSA-GEF": prof([GEF_P])} + participants = [("R-HSA-vgef", "input", "R-HSA-G", 1), ("R-HSA-vgef", "output", "R-HSA-T", 1), + ("R-HSA-vgef", "catalyst", "R-HSA-GEF", 1), + ("R-HSA-vhyd", "input", "R-HSA-T", 1), ("R-HSA-vhyd", "output", "R-HSA-G", 1), + ("R-HSA-vhyd", "catalyst", "R-HSA-T", 1), + ("R-HSA-vbind", "input", "R-HSA-T", 1), ("R-HSA-vbind", "input", "R-HSA-GAP", 1), + ("R-HSA-vbind", "output", "R-HSA-TX", 1), + ("R-HSA-vrel", "input", "R-HSA-TX", 1), ("R-HSA-vrel", "output", "R-HSA-G", 1), + ("R-HSA-vrel", "output", "R-HSA-GAP", 1)] + return net, rmap, umap, profiles, participants -def test_orientation_by_residues_then_donor_then_components(): - edges, rmap, umap = net() - # residues decide first, whatever the donor says - forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (2, 1), B: (0, 1)}, {F}) - assert dict((u, b) for _, u, _, b in forms)[U["b"]] is True - # no residue difference: the donor sits on the REVERSE reaction, so GDP is modified - forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 2), B: (0, 2)}, {R}) - assert dict((u, b) for _, u, _, b in forms) == {U["a"]: False, U["b"]: True} - # no residue or donor difference: the form with more components is modified - forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 3)}, set()) - assert dict((u, b) for _, u, _, b in forms)[U["a"]] is True +def test_r_steps_flag_self_catalysis_as_non_enzyme(): + _, _, _, profiles, participants = ras_network() + steps = m.r_steps(participants, profiles) + by = {(rx, r): (a, b, enz) for rx, r, a, b, enz in steps} + assert by[("R-HSA-vgef", RAS)] == ("R-HSA-G", "R-HSA-T", True) # GEF catalyst: enzyme + assert by[("R-HSA-vhyd", RAS)] == ("R-HSA-T", "R-HSA-G", False) # catalysed by RAS:GTP itself + assert by[("R-HSA-vbind", RAS)] == ("R-HSA-T", "R-HSA-TX", True) # joining input: the GAP + assert by[("R-HSA-vrel", RAS)] == ("R-HSA-TX", "R-HSA-G", False) # release: nothing joins + assert by[("R-HSA-vbind", GAP_P)] == ("R-HSA-GAP", "R-HSA-TX", True) # for the GAP, RAS:GTP joins + assert ("R-HSA-vrel", GAP_P) in by + + +def test_ras_pool_has_gdp_base_gtp_state_and_the_gap_complex_as_intermediate(): + net, rmap, umap, profiles, participants = ras_network() + st = {} + forms, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, + {"R-HSA-vgef"}, {}, st) + assert [(u, role, base) for _, u, _, role, base in forms] == [ + (U["g"], "state", True), (U["t"], "state", False), (U["tx"], "intermediate", False)] + paths = {} + for _, pid, step, a, b, rx, rs, enz in trans: + paths.setdefault(pid, []).append((step, a, b, rx, rs, enz)) + shapes = sorted(tuple((a, b, rs, enz) for _, a, b, _, rs, enz in sorted(p)) for p in paths.values()) + assert shapes == sorted([ + ((U["g"], U["t"], "R-HSA-vgef", True),), + ((U["t"], U["g"], "R-HSA-vhyd", False),), + ((U["t"], U["tx"], "R-HSA-vbind", True), (U["tx"], U["g"], "R-HSA-vrel", False)), + ]) + assert carriers == [("pool1", U["gap"], U["vrel"])] + assert st["pools"] == 1 and st["states"] == 2 and st["intermediates"] == 1 and st["paths"] == 3 + assert st["multi_step_pools"] == 1 and st["carriers"] == 1 and st["ties"] == 0 + assert st["carrier_loops"] == 1 # the GAP's own GAP -> RAS:GTP:GAP -> GAP loop + assert st["autocat_source"] == 1 and st["autocat_product"] == 0 and st["autocat_other"] == 0 + + +def test_a_step_whose_node_loop_is_broken_is_not_a_transition(): + net, rmap, umap, profiles, participants = ras_network() + # the release reaction returns a DIFFERENT copy of RAS:GDP + net.loc[(net["source_id"] == U["vrel"]) & (net["target_id"] == U["g"]), "target_id"] = U["a2"] + umap[U["a2"]] = "R-HSA-G" + st = {} + forms, trans, _ = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert {u for _, u, _, _, _ in forms} == {U["g"], U["t"]} and st["paths"] == 2 and st["intermediates"] == 0 -def test_undecided_pool_is_counted_as_a_tie(): - edges, rmap, umap = net() +def test_a_form_node_that_is_a_member_of_the_steps_set_maps_to_it(): + net, rmap, umap, profiles, participants = ras_network() + participants = [(rx, role, "R-HSA-GSET" if pe == "R-HSA-G" else pe, s) for rx, role, pe, s in participants] + profiles["R-HSA-GSET"] = profiles["R-HSA-G"] + steps = m.r_steps(participants, profiles) + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, {})[0] == [] + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {"R-HSA-GSET": {"R-HSA-G"}}, {}) + assert len(forms) == 3 + + +# --- six-form ring S -> S:E -> S*:E -> S* -> S*:P -> S:P -> S ----------------------------------------------- + +def ring_network(): + edges = [e("s", "v1", "input"), e("e", "v1", "input"), e("v1", "se", "output"), + e("se", "v2", "input"), e("v2", "sxe", "output"), + e("sxe", "v3", "input"), e("v3", "sx", "output"), e("v3", "e", "output"), + e("sx", "v4", "input"), e("p", "v4", "input"), e("v4", "sxp", "output"), + e("sxp", "v5", "input"), e("v5", "sp", "output"), + e("sp", "v6", "input"), e("v6", "s", "output"), e("v6", "p", "output")] + net, rmap, umap = network(edges, ["v1", "v2", "v3", "v4", "v5", "v6"], + {"s": "R-HSA-S", "se": "R-HSA-SE", "sxe": "R-HSA-SXE", "sx": "R-HSA-SX", + "sxp": "R-HSA-SXP", "sp": "R-HSA-SP", "e": "R-HSA-E", "p": "R-HSA-P"}) + ph = [("res", "phospho-S at 1")] + profiles = {"R-HSA-S": prof([SUB]), "R-HSA-SX": prof([SUB], ph, mods=1), + "R-HSA-E": prof([KIN]), "R-HSA-P": prof([PHOS]), + "R-HSA-SE": {"proteins": {SUB, KIN}, "fixed": {SUB, KIN}, "sig": {SUB: frozenset(), KIN: frozenset()}, + "mods": 0, "comps": 2, "slots": 2}, + "R-HSA-SXE": {"proteins": {SUB, KIN}, "fixed": {SUB, KIN}, "sig": {SUB: frozenset(ph), KIN: frozenset()}, + "mods": 1, "comps": 2, "slots": 2}, + "R-HSA-SXP": {"proteins": {SUB, PHOS}, "fixed": {SUB, PHOS}, "sig": {SUB: frozenset(ph), PHOS: frozenset()}, + "mods": 1, "comps": 2, "slots": 2}, + "R-HSA-SP": {"proteins": {SUB, PHOS}, "fixed": {SUB, PHOS}, "sig": {SUB: frozenset(), PHOS: frozenset()}, + "mods": 0, "comps": 2, "slots": 2}} + participants = [("R-HSA-v1", "input", "R-HSA-S", 1), ("R-HSA-v1", "input", "R-HSA-E", 1), ("R-HSA-v1", "output", "R-HSA-SE", 1), + ("R-HSA-v2", "input", "R-HSA-SE", 1), ("R-HSA-v2", "output", "R-HSA-SXE", 1), + ("R-HSA-v3", "input", "R-HSA-SXE", 1), ("R-HSA-v3", "output", "R-HSA-SX", 1), ("R-HSA-v3", "output", "R-HSA-E", 1), + ("R-HSA-v4", "input", "R-HSA-SX", 1), ("R-HSA-v4", "input", "R-HSA-P", 1), ("R-HSA-v4", "output", "R-HSA-SXP", 1), + ("R-HSA-v5", "input", "R-HSA-SXP", 1), ("R-HSA-v5", "output", "R-HSA-SP", 1), + ("R-HSA-v6", "input", "R-HSA-SP", 1), ("R-HSA-v6", "output", "R-HSA-S", 1), ("R-HSA-v6", "output", "R-HSA-P", 1)] + return net, rmap, umap, profiles, participants + + +def test_six_form_ring_gives_two_states_four_intermediates_two_paths_two_carriers(): + net, rmap, umap, profiles, participants = ring_network() st = {} - forms, _ = m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 1)}, set(), st) - assert st["ties"] == 1 and len(forms) == 2 + forms, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + roles = {u: (role, base) for _, u, _, role, base in forms} + assert roles == {U["s"]: ("state", True), U["sx"]: ("state", False), + U["se"]: ("intermediate", False), U["sxe"]: ("intermediate", False), + U["sxp"]: ("intermediate", False), U["sp"]: ("intermediate", False)} + paths = {} + for _, pid, step, a, b, rx, rs, enz in trans: + paths.setdefault(pid, []).append((step, a, b, rx, enz)) + assert len(paths) == 2 and st["paths"] == 2 and st["multi_step_pools"] == 1 + for p in paths.values(): + assert [s for s, *_ in sorted(p)] == [1, 2, 3] + assert any(enz for *_, enz in p) # the binding step joins the enzyme + assert [enz for *_, enz in sorted(p)] == [True, False, False] + assert {(p[0][1], p[-1][2]) for p in (sorted(v) for v in paths.values())} == {(U["s"], U["sx"]), (U["sx"], U["s"])} + assert sorted(carriers) == sorted([("pool1", U["e"], U["v3"]), ("pool1", U["p"], U["v6"])]) + assert st["carriers"] == 2 and st["carrier_loops"] == 2 and st["ties"] == 0 + + +def test_an_enzymes_own_binding_loop_is_not_a_pool(): + # E -> E:S -> E: one signature for E, so a carrier loop, not a pool + edges = [e("e", "v1", "input"), e("v1", "se", "output"), e("se", "v3", "input"), e("v3", "e", "output")] + net, rmap, umap = network(edges, ["v1", "v3"], {"e": "R-HSA-E", "se": "R-HSA-SE"}) + profiles = {"R-HSA-E": prof([KIN]), "R-HSA-SE": prof([KIN, SUB], slots=2, comps=2)} + steps = [("R-HSA-v1", KIN, "R-HSA-E", "R-HSA-SE", True), ("R-HSA-v3", KIN, "R-HSA-SE", "R-HSA-E", False)] st = {} - m.find_pools(edges, rmap, umap, PAIRS, {A: (0, 1), B: (0, 1)}, {F}, st) - assert st["ties"] == 0 + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["carrier_loops"] == 1 and st["pools"] == 0 + +# --- shared nodes, path cap, orientation -------------------------------------------------------------------- -def test_a_reaction_converting_two_pairs_is_dropped(): - # vf also converts a2 -> b: its node would be written by two fluxes - edges, rmap, umap = net() - edges = pd.concat([edges, pd.DataFrame([e("a2", "vf", "input"), e("vr", "a2", "output")])], ignore_index=True) +def test_a_node_claimed_by_two_pools_is_dropped_and_counted(): + # X + Y -> XY (both get modified) -> X + Y: the complex is a form of X's pool and of Y's + edges = [e("x", "va", "input"), e("y", "va", "input"), e("va", "xy", "output"), + e("xy", "vb", "input"), e("vb", "x", "output"), e("vb", "y", "output")] + net, rmap, umap = network(edges, ["va", "vb"], {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY"}) + profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), + "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 2, "comps": 2, "slots": 2, + "sig": {1: frozenset([("res", "p1")]), 2: frozenset([("res", "p2")])}}} + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-X", False), + ("R-HSA-va", 2, "R-HSA-Y", "R-HSA-XY", True), ("R-HSA-vb", 2, "R-HSA-XY", "R-HSA-Y", False)] + st = {} + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["shared_nodes_dropped"] == 3 and st["pools"] == 0 # XY and both reaction nodes + + +def test_co_travelling_proteins_of_one_set_are_one_pool(): + # p21 RAS:GDP <-> RAS:GTP: the H/K/NRAS graphs are identical, so one pool, counted as merged + net, rmap, umap, profiles, participants = ras_network() + for s in ("R-HSA-G", "R-HSA-T"): + profiles[s] = prof([RAS, 4], profiles[s]["sig"][RAS], comps=2) + profiles["R-HSA-TX"]["proteins"].add(4) + profiles["R-HSA-TX"]["sig"][4] = profiles["R-HSA-TX"]["sig"][RAS] + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert st["pools"] == 1 and st["merged_proteins"] == 1 and st["shared_nodes_dropped"] == 0 and len(forms) == 3 + + +def chain_network(n_intermediates): + inter = ["i1", "i2", "i3", "i4", "i5", "i6"][:n_intermediates] + seq = ["s"] + inter + ["sx"] + rxs = ["v1", "v2", "v3", "v4", "v5", "v6", "v7"][:len(seq) - 1] + edges, umap, steps = [], {"s": "R-HSA-S", "sx": "R-HSA-SX"}, [] + profiles = {"R-HSA-S": prof([SUB]), "R-HSA-SX": prof([SUB], [("res", "p")], mods=1)} + for k, rx in enumerate(rxs): + a, b = seq[k], seq[k + 1] + edges += [e(a, rx, "input"), e(rx, b, "output")] + for x in (a, b): + if x in inter: + umap[x] = f"R-HSA-{x.upper()}" + profiles[f"R-HSA-{x.upper()}"] = prof([SUB], slots=2, comps=2) + steps.append((f"R-HSA-{rx}", SUB, umap[a], umap[b], True)) + edges += [e("sx", "vback", "input"), e("vback", "s", "output")] + steps.append(("R-HSA-vback", SUB, "R-HSA-SX", "R-HSA-S", True)) + net, rmap, umap = network(edges, rxs + ["vback"], umap) + return net, rmap, umap, steps, profiles + + +def test_paths_over_six_steps_are_dropped_and_counted(): + net, rmap, umap, steps, profiles = chain_network(5) # 6 steps: kept + st = {} + m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 2 and st["long_paths_dropped"] == 0 and st["intermediates"] == 5 + net, rmap, umap, steps, profiles = chain_network(6) # 7 steps: dropped + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 1 and st["long_paths_dropped"] == 1 + assert st["off_path_intermediates"] == 6 and [r for _, _, _, r, _ in forms] == ["state", "state"] + + +def test_orientation_by_residues_then_donor_then_components_and_ties_are_counted(): + net, rmap, umap, profiles, participants = ras_network() + steps = m.r_steps(participants, profiles) + + def base_of(profiles, donors, st=None): + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, donors, {}, st if st is not None else {}) + return next(u for _, u, _, role, base in forms if base) + # no residue difference: the donor is consumed on the way to GTP, so GDP is base + assert base_of(profiles, {"R-HSA-vgef"}) == U["g"] + assert base_of(profiles, {"R-HSA-vhyd"}) == U["t"] + # residues decide first, whatever the donor says + profiles["R-HSA-G"]["mods"] = 2 + assert base_of(profiles, {"R-HSA-vgef"}) == U["t"] + profiles["R-HSA-G"]["mods"] = 0 + # no residue or donor difference: the form with more components is modified + profiles["R-HSA-T"]["comps"] = 3 + assert base_of(profiles, set()) == U["g"] + profiles["R-HSA-T"]["comps"] = 2 + # nothing decides: smaller stId, counted as a tie st = {} - forms, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}, st) - assert forms == [] and trans == [] and st["multi_use"] >= 1 + assert base_of(profiles, set(), st) == U["g"] and st["ties"] == 1 -def test_no_pairs_or_empty_network_give_no_pools(): - edges, rmap, umap = net() - assert m.find_pools(edges, rmap, umap, [], {}, set()) == ([], []) - assert m.find_pools(pd.DataFrame(columns=["source_id", "target_id", "edge_type"]), rmap, umap, PAIRS, {}, set()) == ([], []) +def test_no_steps_or_empty_network_give_no_pools(): + net, rmap, umap, profiles, participants = ras_network() + assert m.find_pools(net, rmap, umap, [], profiles, set()) == ([], [], []) + empty = pd.DataFrame(columns=["source_id", "target_id", "edge_type"]) + assert m.find_pools(empty, rmap, umap, m.r_steps(participants, profiles), profiles, set()) == ([], [], []) -def test_transitions_say_whether_they_are_catalysed(): - edges, rmap, umap = net() - edges = pd.concat([edges, pd.DataFrame([{"source_id": U["a2"], "target_id": U["vf"], "edge_type": "catalyst"}])], - ignore_index=True) - _, trans = m.find_pools(edges, rmap, umap, PAIRS, {}, {F}) - assert {rx: c for _, _, _, rx, _, c in trans} == {U["vf"]: True, U["vr"]: False} +def test_a_stoichiometry_two_input_is_not_a_step(): + _, _, _, profiles, participants = ras_network() + participants = [(rx, role, pe, 2 if (rx, role) == ("R-HSA-vgef", "input") else s) for rx, role, pe, s in participants] + assert ("R-HSA-vgef", RAS) not in {(rx, r) for rx, r, *_ in m.r_steps(participants, profiles)} From ad5e06e02335c152f76b28e174dc26c7d20281b2 Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Mon, 28 Sep 2026 00:36:48 -0400 Subject: [PATCH 5/7] pools: a step's parallel copies are one step (specs/039 amendment 2 contract) A path is now its NODE sequence plus the reaction stId at each step. The reaction nodes of that stId between the step's two nodes (variant and set copies) are the step's parallel copies, one row each in pool_transitions.csv: rows sharing (pool_id, path_id, step) have the same source_uuid, target_uuid, reaction_stid and enzyme_driven and differ in reaction_uuid. A consumer reads a step as the mean of its copies and a path as the product of its steps, which equals the sum over the uuid paths the previous commit enumerated (RAF: 48 bind copies x 48 release copies were 2,304 two-step paths; now one path of 96 rows). Different reactions between the same two nodes stay different paths (GEF exchange beside intrinsic exchange), because each is weighed on its own as enzyme-driven or not. POOL_PATH_BUDGET is kept as a guard on the path count; nothing hits it now (before: HRR and EGFR pool1 did). The log adds the step-copy row count. Interpretations of the pre-registered rule, accepted by the deltasignal side and recorded here: - "a catalyst or joining input containing a protein other than R" is read as containing a protein and NOT containing R, so the set p21 RAS:GTP catalysing intrinsic hydrolysis is self-catalysis (non-enzyme) for each RAS member; - "the complex whose direct children include R's leaf" looks through set membership (R97 puts RAS in a set inside p21 RAS:GTP); GAPs in RAS:GTP:GAP get an empty signature, as the derivation checked; - enzyme_driven is written per step; a path is enzyme-driven if any step is; - co-travelling reference entities with identical pools (paralogs of one set) are one pool, counted as merged; only DIFFERENT pools sharing a node trigger the shared-node drop; - states and intermediates on no kept path are omitted and counted; a candidate whose kept states have fewer than two signatures is dropped; - a variant node's proteins are its parent's non-set proteins plus its chosen members'; its signature is the parent's; - the base rule's residues are summed over all leaves (amendment 1 as written), not R's leaves only. tests/test_pools.py: 13 tests; the new one asserts that two reaction nodes of one stId collapse into one path with one row per copy while an intrinsic exchange beside the GEF stays a separate path. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 90 ++++++++++++++++++++-------------- tests/test_pools.py | 36 ++++++++++++++ 2 files changed, 89 insertions(+), 37 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 87f3104..746a185 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3410,11 +3410,11 @@ def parts(node_id: str): POOL_MAX_STEPS = 6 # a state -> intermediates -> state path longer than this is dropped -POOL_PATH_BUDGET = 20000 # a pool with more uuid paths than this is dropped (a safety cap, counted) +POOL_PATH_BUDGET = 20000 # a pool with more paths (node + reaction sequences) than this is dropped and counted POOL_STATS_KEYS = ("pools", "states", "intermediates", "paths", "multi_step_pools", "long_paths_dropped", "shared_nodes_dropped", "ties", "carriers", "autocat_source", "autocat_product", "autocat_other", "carrier_loops", "no_path_pools", "off_path_intermediates", - "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins") + "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins", "copy_rows") def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, @@ -3439,11 +3439,14 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram loop: an enzyme's E -> E:S -> E). In each signature the forms with the fewest non-small slots are states, the rest intermediates. Transitions are simple directed paths state -> intermediates* -> another state of at most - ``POOL_MAX_STEPS`` steps; each uuid path is one parallel transition, and - longer ones are dropped and counted. Because an intermediate is one node - however many reaction copies enter and leave it, the uuid paths through it - are the product of the copies per step (RAF: 48 GAP binding copies x 48 - release copies = 2,304 two-step paths); a pool with more than + ``POOL_MAX_STEPS`` steps, longer ones dropped and counted. A path's + identity is its node sequence plus the reaction stId at each step; the + reaction NODES of that stId between the step's two nodes (variant and set + copies) are the step's parallel copies and are listed one row each, so a + consumer reads a step as the mean of its copies and a path as the product + of its steps. Different reactions between the same two nodes (GEF + exchange beside intrinsic exchange) stay different paths, since each is + weighed on its own as enzyme-driven or not. A pool with more than ``POOL_PATH_BUDGET`` paths is dropped and counted (``budget_dropped``). The base state follows amendment 1 (residues, then the donor-consuming direction, then components; undecided pools fall back to the smaller stId and are counted). Pools identical for @@ -3457,7 +3460,9 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram Returns (forms, transitions, carriers): forms: [(pool_id, node_uuid, stable_id, role, is_base)], role state | intermediate transitions: [(pool_id, path_id, step, source_uuid, target_uuid, reaction_uuid, reaction_stid, enzyme_driven)] - ``enzyme_driven`` is the STEP's flag; a path is enzyme-driven if any step is. + one row per (path, step, copy): rows sharing (pool_id, path_id, step) have the same + source, target and reaction_stid and differ in reaction_uuid. ``enzyme_driven`` is the + STEP's flag; a path is enzyme-driven if any step is. carriers: [(pool_id, carrier_uuid, release_reaction_uuid)] ``stats`` receives the counts in ``POOL_STATS_KEYS``; the autocatalysis counts are distinct (reaction node, catalyst) pairs where a form of the pool @@ -3572,21 +3577,28 @@ def build(removed: Set[str]): for members in by_sig.values(): fewest = min(prof(u, "slots") for u in members) states |= {u for u in members if prof(u, "slots") == fewest} - adj: Dict[str, List[Tuple[str, str, str, bool]]] = {} + # a step is (source node, target node, reaction stId); the reaction + # nodes of that stId between those nodes are its parallel copies + copies_of: Dict[Tuple[str, str, str], List[str]] = {} + enz_of: Dict[Tuple[str, str, str], bool] = {} for a_u, b_u, rx_u, rs, enz in sorted(by_r[r]): if a_u in scc and b_u in scc: - adj.setdefault(a_u, []).append((b_u, rx_u, rs, enz)) - paths: List[Tuple[Tuple[str, str, str, str, bool], ...]] = [] + copies_of.setdefault((a_u, b_u, rs), []).append(rx_u) + enz_of[(a_u, b_u, rs)] = enz_of.get((a_u, b_u, rs), False) or enz + adj: Dict[str, List[Tuple[str, str]]] = {} + for a_u, b_u, rs in sorted(copies_of): + adj.setdefault(a_u, []).append((b_u, rs)) + paths: List[Tuple[Tuple[str, str, str, bool, Tuple[str, ...]], ...]] = [] over_budget = False def walk(u, so_far, seen): nonlocal over_budget - for v, rx_u, rs, enz in adj.get(u, ()): + for v, rs in adj.get(u, ()): if over_budget: return if v in seen: continue - new = so_far + ((u, v, rx_u, rs, enz),) + new = so_far + ((u, v, rs, enz_of[(u, v, rs)], tuple(copies_of[(u, v, rs)])),) if v in states: paths.append(new) if len(paths) > POOL_PATH_BUDGET: @@ -3614,13 +3626,13 @@ def walk(u, so_far, seen): local["off_path_states"] += len(states - on_path) local["off_path_intermediates"] += len(scc - states - on_path) key = (frozenset(kept_states), frozenset(inter), - frozenset(tuple((a_u, b_u, rx_u) for a_u, b_u, rx_u, _, _ in p) for p in paths)) + frozenset(tuple((a_u, b_u, rs, cp) for a_u, b_u, rs, _, cp in p) for p in paths)) if key in pools: pools[key]["r"].add(r) local["merged_proteins"] += 1 else: pools[key] = {"r": {r}, "states": kept_states, "inter": inter, "paths": paths, - "nodes": kept_states | inter | {st[2] for p in paths for st in p}} + "nodes": kept_states | inter | {c for p in paths for st in p for c in st[4]}} return list(pools.values()), local removed: Set[str] = set() @@ -3640,15 +3652,15 @@ def walk(u, so_far, seen): stats["shared_nodes_dropped"] = dropped_shared def path_key(p): - return tuple((ent.get(a, ""), a, ent.get(b, ""), b, rx) for a, b, rx, _, _ in p) + return tuple((ent.get(a, ""), a, ent.get(b, ""), b, rs) for a, b, rs, _, _ in p) def modified_side(x, y, xy, yx): """The state a pair marks as modified, or None when undecidable.""" mx, my = prof(x, "mods"), prof(y, "mods") if mx != my: return y if my > mx else x - fd = any(rs in donor_reactions for p in xy for _, _, _, rs, _ in p) - rd = any(rs in donor_reactions for p in yx for _, _, _, rs, _ in p) + fd = any(rs in donor_reactions for p in xy for _, _, rs, _, _ in p) + rd = any(rs in donor_reactions for p in yx for _, _, rs, _, _ in p) if fd != rd: return y if fd else x cx, cy = prof(x, "comps"), prof(y, "comps") @@ -3678,28 +3690,30 @@ def modified_side(x, y, xy, yx): for u in sorted(states | inter, key=lambda u: (ent.get(u, ""), u)): forms.append((pid, u, ent.get(u, ""), "state" if u in states else "intermediate", u == base)) for n, path in enumerate(paths, start=1): - for k, (a_u, b_u, rx_u, rs, enz) in enumerate(path, start=1): - transitions.append((pid, f"{pid}_p{n}", k, a_u, b_u, rx_u, rs, enz)) - rx_on_paths = sorted({st[2] for path in paths for st in path}) + for k, (a_u, b_u, rs, enz, cp) in enumerate(path, start=1): + for rx_u in sorted(cp): + transitions.append((pid, f"{pid}_p{n}", k, a_u, b_u, rx_u, rs, enz)) + rx_on_paths = sorted({c for path in paths for st in path for c in st[4]}) found: Set[Tuple[str, str]] = set() autocat: Set[Tuple[str, str, str]] = set() for path in paths: - for a_u, b_u, rx_u, _, enz in path: - for c in cats.get(rx_u, ()): - if c in p["nodes"]: # a form of the pool catalyses a step of its own pool - autocat.add((rx_u, c, "autocat_source" if c == a_u else "autocat_product" if c == b_u - else "autocat_other")) - if not enz: - continue - for j in ins.get(rx_u, {}): - if j == a_u or j in p["nodes"] or pe_of(j) in _UBIQUITIN_STIDS: - continue - pj = proteins(j) - if not pj or pj & rset: + for a_u, b_u, _, enz, cp in path: + for rx_u in cp: + for c in cats.get(rx_u, ()): + if c in p["nodes"]: # a form of the pool catalyses a step of its own pool + autocat.add((rx_u, c, "autocat_source" if c == a_u else "autocat_product" if c == b_u + else "autocat_other")) + if not enz: continue - for rel in rx_on_paths: - if j in outs.get(rel, {}): - found.add((j, rel)) + for j in ins.get(rx_u, {}): + if j == a_u or j in p["nodes"] or pe_of(j) in _UBIQUITIN_STIDS: + continue + pj = proteins(j) + if not pj or pj & rset: + continue + for rel in rx_on_paths: + if j in outs.get(rel, {}): + found.add((j, rel)) carriers.extend((pid, j, rel) for j, rel in sorted(found, key=lambda t: (ent.get(t[0], ""), t[0], t[1]))) for _, _, kind in autocat: stats[kind] += 1 @@ -3707,6 +3721,7 @@ def modified_side(x, y, xy, yx): stats["states"] += len(states) stats["intermediates"] += len(inter) stats["paths"] += len(paths) + stats["copy_rows"] += sum(len(cp) for path in paths for *_, cp in path) stats["multi_step_pools"] += int(any(len(path) > 1 for path in paths)) stats["carriers"] += len({j for j, _ in found}) return forms, transitions, carriers @@ -3776,7 +3791,8 @@ def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_ pd.DataFrame(carriers, columns=["pool_id", "carrier_uuid", "release_reaction_uuid"]).to_csv(carriers_file, index=False) logger.info( f"Exported {stats['pools']} pools ({stats['states']} states, {stats['intermediates']} intermediates, " - f"{stats['paths']} paths, {stats['multi_step_pools']} multi-step; {len(steps)} R-steps at stId level; " + f"{stats['paths']} paths in {stats['copy_rows']} step-copy rows, {stats['multi_step_pools']} multi-step; " + f"{len(steps)} R-steps at stId level; " f"{stats['long_paths_dropped']} paths over {POOL_MAX_STEPS} steps dropped, " f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['ties']} oriented by tie-break, " f"{stats['carriers']} carriers, {stats['carrier_loops']} carrier loops, " diff --git a/tests/test_pools.py b/tests/test_pools.py index d78c001..8f753b0 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -270,3 +270,39 @@ def test_a_stoichiometry_two_input_is_not_a_step(): _, _, _, profiles, participants = ras_network() participants = [(rx, role, pe, 2 if (rx, role) == ("R-HSA-vgef", "input") else s) for rx, role, pe, s in participants] assert ("R-HSA-vgef", RAS) not in {(rx, r) for rx, r, *_ in m.r_steps(participants, profiles)} + + +def test_copies_of_one_reaction_collapse_into_one_path_with_one_row_per_copy(): + # RAS with the GAP binding curated for two GAP variants (two reaction NODES of + # one stId) and both release copies, plus an intrinsic exchange beside the GEF: + # copies of one stId are one step (k rows); different reactions stay apart. + net, rmap, umap, profiles, participants = ras_network() + extra = [e("t", "v7", "input"), e("a2", "v7", "input"), e("v7", "tx", "output"), + e("tx", "v8", "input"), e("v8", "g", "output"), e("v8", "a2", "output"), + e("g", "vback", "input"), e("vback", "t", "output")] + net = pd.concat([net, pd.DataFrame(extra)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v7"], U["v8"], U["vback"]], + "reactome_id": ["R-HSA-vbind", "R-HSA-vrel", "R-HSA-vint"]})], ignore_index=True) + umap[U["a2"]] = "R-HSA-GAP2" + profiles["R-HSA-GAP2"] = prof([GAP_P]) + participants += [("R-HSA-vint", "input", "R-HSA-G", 1), ("R-HSA-vint", "output", "R-HSA-T", 1)] + st = {} + _, trans, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + rows = {} + for _, pid, step, a, b, rx, rs, enz in trans: + rows.setdefault((pid, step), []).append((a, b, rx, rs, enz)) + by_shape = {} + for (pid, step), rs_ in rows.items(): + assert len({(a, b, rs, enz) for a, b, _, rs, enz in rs_}) == 1 # one source, target, stId, flag per step + a, b, _, rs, enz = rs_[0] + by_shape.setdefault(pid, []).append((step, a, b, rs, enz, sorted(rx for _, _, rx, _, _ in rs_))) + shapes = sorted(tuple(x[1:] for x in sorted(v)) for v in by_shape.values()) + assert shapes == sorted([ + ((U["g"], U["t"], "R-HSA-vgef", True, [U["vgef"]]),), + ((U["g"], U["t"], "R-HSA-vint", False, [U["vback"]]),), # not merged into the GEF step + ((U["t"], U["g"], "R-HSA-vhyd", False, [U["vhyd"]]),), + ((U["t"], U["tx"], "R-HSA-vbind", True, sorted([U["vbind"], U["v7"]])), + (U["tx"], U["g"], "R-HSA-vrel", False, sorted([U["vrel"], U["v8"]]))), + ]) + assert st["paths"] == 4 and st["copy_rows"] == 7 and st["multi_step_pools"] == 1 + assert sorted(carriers) == sorted([("pool1", U["gap"], U["vrel"]), ("pool1", U["a2"], U["v8"])]) From bdaed8edf23d878a54fceec6f8950fe7dfd5fb3b Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Mon, 28 Sep 2026 01:15:22 -0400 Subject: [PATCH 6/7] pools: core-only states (amendment 3), no carrier a pool manages, long paths counted whole Pre-arm review of specs/039 amendment 2 found two generator defects. 1. States from annotation-dependent signatures. Activated FGFR4 <-> FGFR4:PLCG1 <-> FGFR4:p-4Y-PLCG1 (R-HSA-5654743) was a three-state pool of FGFR4, and EGF:p-EGFR:p-ERBB2 / :PLCG1 / p-ERBB2 heterodimers:PTK6 (R-HSA-1227986) three states of ERBB2, because a bound partner's phosphorylation is a signature change for the partner and the receptor's forms differed only by what was bound. Rule (amendment 3): the CORE of a strongly connected set is the proteins every form of it carries; a form is eligible to be a state only if its proteins are exactly the core (small molecules allowed); states are per signature among those; the set is a pool only if its core-only forms carry >= 2 signatures; every other form is an intermediate. An enzyme's own loop has one core-only form, so it is not a pool, as before. On the pools039b build: 5654743 dropped; 1227986 now has two states (p-6Y-EGFR:p-6Y,Y1112-ERBB2 base, p-7Y,Y1112-ERBB2) with the PLCG1 and PTK6 complexes as intermediates; the other 39 pools unchanged. 2. A carrier that is a node of another pool (R-HSA-74752: SHC1's carrier Insulin:p-6Y-IR is the receptor pool's intermediate) is no longer a carrier: nodes any pool manages (states, intermediates, step reaction nodes) are excluded from every carrier list and counted (carrier_conflicts; 1 on the build). 3. long_paths_dropped counts dropped paths (state -> ... -> state over 6 steps), not truncated prefixes. tests/test_pools.py: 15 tests. The shared-node fixture is now two genuine pools (X and Y both phosphorylated inside X:Y); new tests cover the FGFR4 shape (not a pool) against the receptor-phosphorylated shape (a pool), the carrier conflict, and whole-path counting. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 45 +++++++++++----- tests/test_pools.py | 99 ++++++++++++++++++++++++++++++---- 2 files changed, 121 insertions(+), 23 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 746a185..0ea31e2 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3414,7 +3414,8 @@ def parts(node_id: str): POOL_STATS_KEYS = ("pools", "states", "intermediates", "paths", "multi_step_pools", "long_paths_dropped", "shared_nodes_dropped", "ties", "carriers", "autocat_source", "autocat_product", "autocat_other", "carrier_loops", "no_path_pools", "off_path_intermediates", - "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins", "copy_rows") + "off_path_states", "ambiguous_copies", "budget_dropped", "merged_proteins", "copy_rows", + "carrier_conflicts") def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, @@ -3434,10 +3435,14 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram A step exists at node level only as a reaction node whose R-containing input node and output node map to the step's forms (the node's entity, its variant parent, or a member of the step's set). Per R, the R-graph over - form nodes is condensed; a strongly connected set with >= 2 forms and >= 2 - modification signatures is a pool candidate (one signature is a carrier - loop: an enzyme's E -> E:S -> E). In each signature the forms with the - fewest non-small slots are states, the rest intermediates. Transitions are + form nodes is condensed. In a strongly connected set with >= 2 forms, the + CORE is the proteins every form carries; only a form made of exactly the + core (small molecules allowed) can be a state, and the set is a pool only + if its core-only forms carry >= 2 modification signatures (amendment 3; + one signature is a carrier loop: an enzyme's E -> E:S -> E, and a set + whose forms differ only by what else is bound is not a pool). In each + signature the core-only forms with the fewest non-small slots are states, + the rest intermediates. Transitions are simple directed paths state -> intermediates* -> another state of at most ``POOL_MAX_STEPS`` steps, longer ones dropped and counted. A path's identity is its node sequence plus the reaction stId at each step; the @@ -3455,7 +3460,8 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram recomputed, and the node is counted. Carriers are the enzyme's free-form nodes (a joining input of an enzyme-driven step, containing a protein and not R, that a step of the pool releases) with the reaction nodes releasing - them. + them; a carrier that any pool manages (state, intermediate or step + node) is excluded and counted (``carrier_conflicts``). Returns (forms, transitions, carriers): forms: [(pool_id, node_uuid, stable_id, role, is_base)], role state | intermediate @@ -3567,9 +3573,15 @@ def build(removed: Set[str]): for scc in nx.strongly_connected_components(g): if len(scc) < 2: continue + # amendment 3: the core is the proteins every form of the set carries; + # only a form made of exactly the core (small molecules allowed) can + # be a state, and the set is a pool only if such forms carry >= 2 + # signatures. An enzyme's loop has one core-only form: not a pool. + core = frozenset.intersection(*(proteins(u) for u in scc)) by_sig: Dict[frozenset, List[str]] = {} for u in scc: - by_sig.setdefault(sig(u, r), []).append(u) + if proteins(u) == core: + by_sig.setdefault(sig(u, r), []).append(u) if len(by_sig) < 2: local["carrier_loops"] += 1 continue @@ -3600,14 +3612,14 @@ def walk(u, so_far, seen): continue new = so_far + ((u, v, rs, enz_of[(u, v, rs)], tuple(copies_of[(u, v, rs)])),) if v in states: - paths.append(new) - if len(paths) > POOL_PATH_BUDGET: + if len(new) > POOL_MAX_STEPS: + local["long_paths_dropped"] += 1 # a whole path, not a prefix + else: + paths.append(new) + if len(paths) + local["long_paths_dropped"] > POOL_PATH_BUDGET: over_budget = True return continue - if len(new) >= POOL_MAX_STEPS: - local["long_paths_dropped"] += 1 - continue walk(v, new, seen | {v}) for s0 in sorted(states): walk(s0, (), {s0}) @@ -3649,6 +3661,7 @@ def walk(u, so_far, seen): removed |= shared dropped_shared += len(shared) stats.update(local) + managed: Set[str] = set().union(*(p["nodes"] for p in pools)) if pools else set() stats["shared_nodes_dropped"] = dropped_shared def path_key(p): @@ -3714,6 +3727,11 @@ def modified_side(x, y, xy, yx): for rel in rx_on_paths: if j in outs.get(rel, {}): found.add((j, rel)) + # a carrier that a pool manages (a state, an intermediate or a step node of + # any pool) is not a carrier: its value is that pool's to write + conflicts = {j for j, _ in found if j in managed} + stats["carrier_conflicts"] += len(conflicts) + found = {(j, rel) for j, rel in found if j not in conflicts} carriers.extend((pid, j, rel) for j, rel in sorted(found, key=lambda t: (ent.get(t[0], ""), t[0], t[1]))) for _, _, kind in autocat: stats[kind] += 1 @@ -3795,7 +3813,8 @@ def export_pools(pathway_id: str, pathway_logic_network: pd.DataFrame, reaction_ f"{len(steps)} R-steps at stId level; " f"{stats['long_paths_dropped']} paths over {POOL_MAX_STEPS} steps dropped, " f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['ties']} oriented by tie-break, " - f"{stats['carriers']} carriers, {stats['carrier_loops']} carrier loops, " + f"{stats['carriers']} carriers ({stats['carrier_conflicts']} excluded as pool nodes), " + f"{stats['carrier_loops']} carrier loops, " f"{stats['merged_proteins']} co-travelling proteins merged, {stats['no_path_pools']} candidates with no path, " f"{stats['off_path_states']}/{stats['off_path_intermediates']} states/intermediates off every path, " f"{stats['ambiguous_copies']} ambiguous reaction copies, {stats['budget_dropped']} over the path budget; " diff --git a/tests/test_pools.py b/tests/test_pools.py index 8f753b0..bd1914f 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -11,7 +11,8 @@ # _uuid_to_stable_id_map detects the mapping direction by uuid SHAPE NAMES = ["g", "t", "tx", "gap", "gef", "vgef", "vhyd", "vbind", "vrel", "a2", "s", "se", "sxe", "sx", "sxp", "sp", "e", "p", "v1", "v2", "v3", "v4", "v5", "v6", - "x", "y", "xy", "va", "vb", "i1", "i2", "i3", "i4", "i5", "i6", "v7", "v8", "vback"] + "x", "y", "xy", "va", "vb", "i1", "i2", "i3", "i4", "i5", "i6", "v7", "v8", "vback", + "xsy", "xs", "ys", "vc", "vd", "ve"] U = {k: f"00000000-0000-0000-0000-{i:012d}" for i, k in enumerate(NAMES, 1)} RAS, GAP_P, GEF_P, SUB, KIN, PHOS = 1, 2, 3, 10, 11, 12 # reference entities @@ -178,18 +179,89 @@ def test_an_enzymes_own_binding_loop_is_not_a_pool(): # --- shared nodes, path cap, orientation -------------------------------------------------------------------- def test_a_node_claimed_by_two_pools_is_dropped_and_counted(): - # X + Y -> XY (both get modified) -> X + Y: the complex is a form of X's pool and of Y's + # X + Y -> X:Y -> X*:Y* -> X* + Y*; X* -> X and Y* -> Y directly. X's pool and + # Y's pool both claim the two complexes and the three reaction nodes. edges = [e("x", "va", "input"), e("y", "va", "input"), e("va", "xy", "output"), - e("xy", "vb", "input"), e("vb", "x", "output"), e("vb", "y", "output")] - net, rmap, umap = network(edges, ["va", "vb"], {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY"}) - profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), - "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 2, "comps": 2, "slots": 2, - "sig": {1: frozenset([("res", "p1")]), 2: frozenset([("res", "p2")])}}} - steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-X", False), - ("R-HSA-va", 2, "R-HSA-Y", "R-HSA-XY", True), ("R-HSA-vb", 2, "R-HSA-XY", "R-HSA-Y", False)] + e("xy", "vb", "input"), e("vb", "xsy", "output"), + e("xsy", "vc", "input"), e("vc", "xs", "output"), e("vc", "ys", "output"), + e("xs", "vd", "input"), e("vd", "x", "output"), e("ys", "ve", "input"), e("ve", "y", "output")] + net, rmap, umap = network(edges, ["va", "vb", "vc", "vd", "ve"], + {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY", "xsy": "R-HSA-XSY", + "xs": "R-HSA-XS", "ys": "R-HSA-YS"}) + p1, p2 = [("res", "p1")], [("res", "p2")] + profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), "R-HSA-XS": prof([1], p1, mods=1), "R-HSA-YS": prof([2], p2, mods=1), + "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 0, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset()}}, + "R-HSA-XSY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 2, "comps": 2, "slots": 2, + "sig": {1: frozenset(p1), 2: frozenset(p2)}}} + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-XS", False), ("R-HSA-vd", 1, "R-HSA-XS", "R-HSA-X", True), + ("R-HSA-va", 2, "R-HSA-Y", "R-HSA-XY", True), ("R-HSA-vb", 2, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 2, "R-HSA-XSY", "R-HSA-YS", False), ("R-HSA-ve", 2, "R-HSA-YS", "R-HSA-Y", True)] + # each alone is a pool + for keep in (1, 2): + st = {} + m.find_pools(net, rmap, umap, [t for t in steps if t[1] == keep], profiles, set(), {}, st) + assert st["pools"] == 1 and st["intermediates"] == 2 st = {} assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) - assert st["shared_nodes_dropped"] == 3 and st["pools"] == 0 # XY and both reaction nodes + assert st["shared_nodes_dropped"] == 5 and st["pools"] == 0 # X:Y, X*:Y* and the three reaction nodes + + +def test_forms_differing_only_by_a_bound_partner_are_not_a_pool(): + # Activated FGFR4 <-> FGFR4:PLCG1 <-> FGFR4:p-PLCG1 (R-HSA-5654743): the PLCG1 + # phosphorylation gives FGFR4 no new signature; the only core-only form is + # FGFR4 itself, so this is a carrier loop, not a pool (amendment 3) + edges = [e("x", "va", "input"), e("y", "va", "input"), e("va", "xy", "output"), + e("xy", "vb", "input"), e("vb", "xsy", "output"), + e("xsy", "vc", "input"), e("vc", "x", "output"), e("vc", "ys", "output")] + net, rmap, umap = network(edges, ["va", "vb", "vc"], + {"x": "R-HSA-X", "y": "R-HSA-Y", "xy": "R-HSA-XY", "xsy": "R-HSA-XSY", "ys": "R-HSA-YS"}) + p2 = [("res", "p2")] + profiles = {"R-HSA-X": prof([1]), "R-HSA-Y": prof([2]), "R-HSA-YS": prof([2], p2, mods=1), + "R-HSA-XY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 0, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset()}}, + "R-HSA-XSY": {"proteins": {1, 2}, "fixed": {1, 2}, "mods": 1, "comps": 2, "slots": 2, + "sig": {1: frozenset(), 2: frozenset(p2)}}} + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-X", False)] + st = {} + assert m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) == ([], [], []) + assert st["carrier_loops"] == 1 + # and the same shape where the RECEPTOR is what gets phosphorylated IS a pool + profiles["R-HSA-XSY"]["sig"] = {1: frozenset(p2), 2: frozenset()} + profiles["R-HSA-XS"] = prof([1], p2, mods=1) + umap[U["ys"]] = "R-HSA-XS" + steps = [("R-HSA-va", 1, "R-HSA-X", "R-HSA-XY", True), ("R-HSA-vb", 1, "R-HSA-XY", "R-HSA-XSY", False), + ("R-HSA-vc", 1, "R-HSA-XSY", "R-HSA-XS", False)] + edges += [e("ys", "vd", "input"), e("vd", "x", "output")] + net, rmap, _ = network(edges, ["va", "vb", "vc", "vd"], {}) + steps.append(("R-HSA-vd", 1, "R-HSA-XS", "R-HSA-X", True)) + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["pools"] == 1 and {u: r for _, u, _, r, _ in forms} == { + U["x"]: "state", U["ys"]: "state", U["xy"]: "intermediate", U["xsy"]: "intermediate"} + + +def test_a_carrier_that_another_pool_manages_is_excluded_and_counted(): + # the ring's kinase E is also a form of its own pool (E <-> E* by an outside + # kinase and phosphatase): E is not a carrier of the substrate's pool + net, rmap, umap, profiles, participants = ring_network() + edges = [e("e", "v7", "input"), e("v7", "i1", "output"), e("i1", "v8", "input"), e("v8", "e", "output")] + net = pd.concat([net, pd.DataFrame(edges)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v7"], U["v8"]], "reactome_id": ["R-HSA-v7", "R-HSA-v8"]})], + ignore_index=True) + umap[U["i1"]] = "R-HSA-ES" + profiles["R-HSA-ES"] = prof([KIN], [("res", "pE")], mods=1) + participants += [("R-HSA-v7", "input", "R-HSA-E", 1), ("R-HSA-v7", "output", "R-HSA-ES", 1), + ("R-HSA-v7", "catalyst", "R-HSA-P", 1), + ("R-HSA-v8", "input", "R-HSA-ES", 1), ("R-HSA-v8", "output", "R-HSA-E", 1), + ("R-HSA-v8", "catalyst", "R-HSA-P", 1)] + st = {} + forms, _, carriers = m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) + assert st["pools"] == 2 and st["shared_nodes_dropped"] == 0 + assert carriers == [(next(p for p, u, *_ in forms if u == U["s"]), U["p"], U["v6"])] + assert st["carrier_conflicts"] == 1 and st["carriers"] == 1 def test_co_travelling_proteins_of_one_set_are_one_pool(): @@ -234,6 +306,13 @@ def test_paths_over_six_steps_are_dropped_and_counted(): forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) assert st["paths"] == 1 and st["long_paths_dropped"] == 1 assert st["off_path_intermediates"] == 6 and [r for _, _, _, r, _ in forms] == ["state", "state"] + # a second reaction for one of its steps is a second dropped PATH, not a prefix + net = pd.concat([net, pd.DataFrame([e("i3", "v8", "input"), e("v8", "i4", "output")])], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v8"]], "reactome_id": ["R-HSA-v8"]})], ignore_index=True) + steps.append(("R-HSA-v8", SUB, "R-HSA-I3", "R-HSA-I4", True)) + st = {} + m.find_pools(net, rmap, umap, steps, profiles, set(), {}, st) + assert st["paths"] == 1 and st["long_paths_dropped"] == 2 def test_orientation_by_residues_then_donor_then_components_and_ties_are_counted(): From 7d6549035f365859d30ae7f1181f23be06c18764 Mon Sep 17 00:00:00 2001 From: Adam Wright <668955+adamjohnwright@users.noreply.github.com> Date: Mon, 28 Sep 2026 02:56:30 -0400 Subject: [PATCH 7/7] pools: mypy clean (networkx stubs ignore, annotated output lists, no reused loop variables) Co-Authored-By: Claude Opus 5.5 (1M context) --- src/logic_network_generator.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 0ea31e2..ba51f0e 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -3475,7 +3475,7 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram catalyses a step of its own pool, split by whether the catalyst is the step's source form, its product form, or another form. """ - import networkx as nx + import networkx as nx # type: ignore[import-untyped] stats = stats if stats is not None else {} for k in POOL_STATS_KEYS: stats[k] = 0 @@ -3681,7 +3681,9 @@ def modified_side(x, y, xy, yx): return y if cy > cx else x return None - forms, transitions, carriers = [], [], [] + forms: List[Tuple] = [] + transitions: List[Tuple] = [] + carriers: List[Tuple] = [] pools.sort(key=lambda p: min(ent.get(u, "") + u for u in p["states"])) for i, p in enumerate(pools, start=1): pid = f"pool{i}" @@ -3702,10 +3704,10 @@ def modified_side(x, y, xy, yx): stats["ties"] += 1 for u in sorted(states | inter, key=lambda u: (ent.get(u, ""), u)): forms.append((pid, u, ent.get(u, ""), "state" if u in states else "intermediate", u == base)) - for n, path in enumerate(paths, start=1): - for k, (a_u, b_u, rs, enz, cp) in enumerate(path, start=1): + for path_no, path in enumerate(paths, start=1): + for step_no, (a_u, b_u, rs, enz, cp) in enumerate(path, start=1): for rx_u in sorted(cp): - transitions.append((pid, f"{pid}_p{n}", k, a_u, b_u, rx_u, rs, enz)) + transitions.append((pid, f"{pid}_p{path_no}", step_no, a_u, b_u, rx_u, rs, enz)) rx_on_paths = sorted({c for path in paths for st in path for c in st[4]}) found: Set[Tuple[str, str]] = set() autocat: Set[Tuple[str, str, str]] = set()