diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index f1335d3..fa175e0 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -499,10 +499,31 @@ def _register_entity_uuid( def _reject_removed_env() -> None: - """Raise if any removed LNG_* flag is set. Safe to call repeatedly.""" + """Raise if any removed LNG_* flag is set, or a live one holds a value it + does not take. Safe to call repeatedly.""" for name, why in sorted(_REMOVED_ENV.items()): if name in os.environ: raise ValueError(f"{name} was removed: {why}") + pool_active_via() + + +_POOL_ACTIVE_VIA = ("direct", "made_from") + + +def pool_active_via() -> str: + """LNG_POOL_ACTIVE_VIA (deltasignal specs/039 amendments 6 and 7): how a + pool's state is found to act downstream, for the base-state orientation. + ``made_from`` (default since deltasignal adopted amendment 7, 2026-09-28): + a catalyst or positive regulator edge into a reaction outside the pool's + steps, directly or through a set-pool node; or an input edge from the state + into a non-step reaction whose output is, or feeds the set-pool node of, + such a catalyst or regulator (CCNA:CDK2 is made into the catalytic + CCNA:p-T160-CDK2 by CAK). ``direct``: the first test only (amendment 6). + Any other value is an error at startup.""" + value = os.environ.get("LNG_POOL_ACTIVE_VIA", "made_from") + if value not in _POOL_ACTIVE_VIA: + raise ValueError(f"LNG_POOL_ACTIVE_VIA must be one of {_POOL_ACTIVE_VIA}, got {value!r}") + return value def _register_phase1( @@ -3415,7 +3436,7 @@ def parts(node_id: str): "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", - "carrier_conflicts", "nonregenerating_paths") + "carrier_conflicts", "nonregenerating_paths", "oriented_by_activity", "activity_fallbacks") def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFrame, @@ -3458,9 +3479,19 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram 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 + ``POOL_PATH_BUDGET`` paths is dropped and counted (``budget_dropped``). The base state is the + resting, INACTIVE form (amendment 6): a state is active if it has a + ``catalyst`` or positive ``regulator`` edge into a reaction node that is + not one of the pool's own steps, directly or through a set-pool node it + feeds by a ``set_member`` edge (a member of a set-valued catalyst); with exactly one active state the base is + the other state (two states) or the non-active state furthest from it in + the state graph (more; a tie falls back). Under ``LNG_POOL_ACTIVE_VIA= + made_from`` (amendment 7) a state is also active if an input edge from it + enters a non-step reaction whose output acts that way on another non-step + reaction. Otherwise amendment 1 decides (residues, then the donor-consuming direction, then components; undecided - pools fall back to the smaller stId and are counted). Pools identical for + pools fall back to the smaller stId and are counted). Counted as + ``oriented_by_activity`` and ``activity_fallbacks``. 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 @@ -3495,6 +3526,10 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram ins: Dict[str, Dict[str, float]] = {} outs: Dict[str, Dict[str, float]] = {} cats: Dict[str, Set[str]] = {} + acts: Dict[str, Set[str]] = {} # node -> reaction nodes it catalyses or positively regulates + set_pools_of: Dict[str, Set[str]] = {} # member node -> the set-pool nodes it feeds (specs/033) + feeds: Dict[str, Set[str]] = {} # node -> reaction nodes it is an input of + active_via = pool_active_via() 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") @@ -3512,6 +3547,13 @@ def find_pools(pathway_logic_network: pd.DataFrame, reaction_id_map: pd.DataFram outs.setdefault(s, {})[t] = st elif et == "catalyst" and t in vr: cats.setdefault(t, set()).add(s) + acts.setdefault(s, set()).add(t) + elif et == "regulator" and t in vr and str(e.get("pos_neg", "pos")) == "pos": + acts.setdefault(s, set()).add(t) + elif et == "set_member": + set_pools_of.setdefault(s, set()).add(t) + if et == "input" and t in vr: + feeds.setdefault(s, set()).add(t) copies: Dict[str, List[str]] = {} for u, r in vr.items(): copies.setdefault(r, []).append(u) @@ -3734,9 +3776,52 @@ def modified_side(x, y, xy, yx): 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 + # amendment 6: phi0 is the share of the ACTIVE form. A state is active if it + # acts downstream of the pool (a catalyst or positive regulator edge into a + # reaction node that is not one of the pool's own steps). With exactly one + # active state the base is the other one (two states) or the non-active + # state furthest from it in the state graph (more); otherwise, or on a + # tie, the residue / donor / component rule decides. + step_nodes = {c for path in paths for st in path for c in st[4]} + + def acts_on(u: str) -> Set[str]: + # a member of a set-valued catalyst reaches the reaction through the + # set's pool node (specs/033), so the pool node's edges are the member's + out = set(acts.get(u, set())) + for sp in set_pools_of.get(u, ()): + out |= acts.get(sp, set()) + return out + active = {u for u in states if acts_on(u) - step_nodes} + if active_via == "made_from": + # amendment 7: the state the catalytic form is made from. An input + # edge into a non-step reaction whose output acts (as above) on + # another non-step reaction, within two reactions. + for u in states - active: + for rx1 in feeds.get(u, set()) - step_nodes: + if any(acts_on(x) - step_nodes - {rx1} for x in outs.get(rx1, {})): + active.add(u) + break + base = None + if len(active) == 1: + a0 = next(iter(active)) + if len(states) == 2: + base = next(u for u in states if u != a0) + else: + sg = nx.Graph() + sg.add_nodes_from(states) + sg.add_edges_from((x, y) for x, y in between) + dist = nx.single_source_shortest_path_length(sg, a0) + far = max(dist.get(u, -1) for u in states if u != a0) + furthest = [u for u in states if u != a0 and dist.get(u, -1) == far] + if len(furthest) == 1: + base = furthest[0] + if base is not None: + stats["oriented_by_activity"] += 1 + else: + stats["activity_fallbacks"] += 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(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 path_no, path in enumerate(paths, start=1): @@ -3850,7 +3935,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['nonregenerating_paths']} non-regenerating multi-step paths dropped, " - f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['ties']} oriented by tie-break, " + f"{stats['shared_nodes_dropped']} shared nodes dropped, {stats['oriented_by_activity']} oriented by the active " + f"state, {stats['activity_fallbacks']} by residues/donor/components of which {stats['ties']} by tie-break, " 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, " diff --git a/src/pathway_generator.py b/src/pathway_generator.py index 897b6bf..fd186e9 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -39,6 +39,7 @@ "LNG_SET_MEMBERS_OR", "LNG_SET_POOL", "LNG_CAP_POOLS", + "LNG_POOL_ACTIVE_VIA", "LNG_DIAGRAM_SET_MEMBER", # Of these three, only LNG_EMIT_ONE_SIDED is a genuine cache gap: it is # read in reaction_generator.decompose_by_reactions, so it changes diff --git a/tests/test_pools.py b/tests/test_pools.py index b4a8b53..96957e5 100644 --- a/tests/test_pools.py +++ b/tests/test_pools.py @@ -5,6 +5,7 @@ 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 pytest import src.logic_network_generator as m @@ -500,3 +501,164 @@ def test_small_molecules_are_exempt_from_regeneration(): st = {} m.find_pools(net, rmap, umap, m.r_steps(participants, profiles), profiles, set(), {}, st) assert st["nonregenerating_paths"] == 1 and st["pools"] == 0 + + +# --- amendment 6: the base is the INACTIVE form ------------------------------------------------------------ + +def two_state(active=(), regulator=None): + """S <-> S* (single steps, donor on the forward step); ``active`` states get a + catalyst edge into a downstream reaction node v7 (or a ``regulator`` edge of + the given sign).""" + edges = [e("s", "v1", "input"), e("v1", "sx", "output"), e("sx", "v2", "input"), e("v2", "s", "output")] + for a in active: + edges.append(e(a, "v7", "regulator" if regulator else "catalyst")) + edges.append(e("i1", "v7", "input")) + net = pd.DataFrame(edges) + if regulator: + net["pos_neg"] = [regulator if r["edge_type"] == "regulator" else "pos" for _, r in net.iterrows()] + _, rmap, umap = network([], ["v1", "v2", "v7"], {"s": "R-HSA-S", "sx": "R-HSA-SX", "i1": "R-HSA-I"}) + profiles = {"R-HSA-S": prof([SUB]), "R-HSA-SX": prof([SUB], [("res", "p")], mods=1), "R-HSA-I": prof([20])} + steps = [("R-HSA-v1", SUB, "R-HSA-S", "R-HSA-SX", True), ("R-HSA-v2", SUB, "R-HSA-SX", "R-HSA-S", True)] + return net, rmap, umap, profiles, steps + + +def base_and_stats(net, rmap, umap, profiles, steps, donors=frozenset({"R-HSA-v1"})): + st = {} + forms, _, _ = m.find_pools(net, rmap, umap, steps, profiles, set(donors), {}, st) + return next(u for _, u, _, _, b in forms if b), st + + +def test_an_inhibitory_phosphorylation_makes_the_modified_state_the_base(): + # the UNMODIFIED state catalyses downstream (CCNA:CDK2 vs CCNA:p-Y15-CDK2): + # it is the active form, so phi0 goes to it and the base is the modified state + base, st = base_and_stats(*two_state(active=("s",))) + assert base == U["sx"] and st["oriented_by_activity"] == 1 and st["activity_fallbacks"] == 0 + # a positive regulator edge counts the same; a negative one does not + base, _ = base_and_stats(*two_state(active=("s",), regulator="pos")) + assert base == U["sx"] + base, st = base_and_stats(*two_state(active=("s",), regulator="neg")) + assert base == U["s"] and st["activity_fallbacks"] == 1 + + +def test_an_activating_phosphorylation_keeps_the_unmodified_base(): + base, st = base_and_stats(*two_state(active=("sx",))) + assert base == U["s"] and st["oriented_by_activity"] == 1 + + +def test_both_or_neither_state_active_falls_back_to_the_residue_rule(): + for active in ((), ("s", "sx")): + base, st = base_and_stats(*two_state(active=active)) + assert base == U["s"] and st["oriented_by_activity"] == 0 and st["activity_fallbacks"] == 1 and st["ties"] == 0 + + +def test_a_states_own_step_does_not_make_it_active(): + # RAS:GTP catalyses its own hydrolysis (a pool step): not "acting downstream" + net, rmap, umap, profiles, participants = ras_network() + base, st = base_and_stats(net, rmap, umap, profiles, m.r_steps(participants, profiles), {"R-HSA-vgef"}) + assert base == U["g"] and st["activity_fallbacks"] == 1 + # RAS:GTP also activating RAF downstream: still base GDP, now by activity + net = pd.concat([net, pd.DataFrame([e("t", "v7", "catalyst"), e("a2", "v7", "input")])], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v7"]], "reactome_id": ["R-HSA-v7"]})], ignore_index=True) + base, st = base_and_stats(net, rmap, umap, profiles, m.r_steps(participants, profiles), {"R-HSA-vgef"}) + assert base == U["g"] and st["oriented_by_activity"] == 1 + + +def test_three_states_take_the_non_active_state_furthest_from_the_active_one(): + # A <-> B <-> C in a line; C active -> base A. With A <-> C as well, B and A tie -> fallback. + def line(extra=()): + edges = [e("s", "v1", "input"), e("v1", "i1", "output"), e("i1", "v2", "input"), e("v2", "s", "output"), + e("i1", "v3", "input"), e("v3", "sx", "output"), e("sx", "v4", "input"), e("v4", "i1", "output"), + e("sx", "v7", "catalyst"), e("i2", "v7", "input")] + list(extra) + rxs = ["v1", "v2", "v3", "v4", "v7"] + (["v5", "v6"] if extra else []) + net, rmap, umap = network(edges, rxs, {"s": "R-HSA-A", "i1": "R-HSA-B", "sx": "R-HSA-C", "i2": "R-HSA-I"}) + profiles = {"R-HSA-A": prof([SUB]), "R-HSA-B": prof([SUB], [("res", "p")], mods=1), + "R-HSA-C": prof([SUB], [("res", "p"), ("res", "q")], mods=2), "R-HSA-I": prof([20])} + steps = [("R-HSA-v1", SUB, "R-HSA-A", "R-HSA-B", True), ("R-HSA-v2", SUB, "R-HSA-B", "R-HSA-A", True), + ("R-HSA-v3", SUB, "R-HSA-B", "R-HSA-C", True), ("R-HSA-v4", SUB, "R-HSA-C", "R-HSA-B", True)] + if extra: + steps += [("R-HSA-v5", SUB, "R-HSA-A", "R-HSA-C", True), ("R-HSA-v6", SUB, "R-HSA-C", "R-HSA-A", True)] + return net, rmap, umap, profiles, steps + base, st = base_and_stats(*line(), donors=set()) + assert base == U["s"] and st["oriented_by_activity"] == 1 and st["states"] == 3 + base, st = base_and_stats(*line([e("s", "v5", "input"), e("v5", "sx", "output"), + e("sx", "v6", "input"), e("v6", "s", "output")]), donors=set()) + assert base == U["s"] and st["activity_fallbacks"] == 1 # a tie: residues decide + + +def test_a_member_of_a_set_valued_catalyst_acts_through_the_set_pool_node(): + # LPIN1 <-> p-S106-LPIN1: the unphosphorylated form is a member of the + # "lipins" set-pool node, and that node catalyses PA dephosphorylation + net, rmap, umap, profiles, steps = two_state() + net = pd.concat([net, pd.DataFrame([e("s", "i2", "set_member"), e("i2", "v7", "catalyst")])], ignore_index=True) + umap[U["i2"]] = "R-HSA-LIPINS" + profiles["R-HSA-LIPINS"] = prof([SUB, 30]) + base, st = base_and_stats(net, rmap, umap, profiles, steps) + assert base == U["sx"] and st["oriented_by_activity"] == 1 + # a set-pool node that acts on nothing outside the pool gives no activity + net.loc[net["edge_type"] == "catalyst", "target_id"] = U["v1"] + base, st = base_and_stats(net, rmap, umap, profiles, steps) + assert base == U["s"] and st["activity_fallbacks"] == 1 + + +# --- amendment 7 (LNG_POOL_ACTIVE_VIA=made_from): the state the catalytic form is made from ------------------ + +def cdk2_shape(): + """CCNA:CDK2 <-> CCNA:p-Y15-CDK2 (the pool); CCNA:CDK2 -> CAK -> CCNA:p-T160-CDK2, + which catalyses a downstream phosphorylation (v8).""" + net, rmap, umap, profiles, steps = two_state() + extra = [e("s", "v5", "input"), e("v5", "i2", "output"), e("i2", "v8", "catalyst"), e("i3", "v8", "input")] + net = pd.concat([net, pd.DataFrame(extra)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v5"], U["v8"]], "reactome_id": ["R-HSA-v5", "R-HSA-v8"]})], + ignore_index=True) + umap[U["i2"]], umap[U["i3"]] = "R-HSA-T160", "R-HSA-ORC" + profiles["R-HSA-T160"] = prof([SUB], [("res", "pT160")], mods=1) + profiles["R-HSA-ORC"] = prof([40]) + return net, rmap, umap, profiles, steps + + +def test_made_from_flips_the_cdk2_shape_and_direct_does_not(monkeypatch): + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "made_from") + base, st = base_and_stats(*cdk2_shape()) + assert base == U["sx"] and st["oriented_by_activity"] == 1 + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "direct") + base, st = base_and_stats(*cdk2_shape()) + assert base == U["s"] and st["activity_fallbacks"] == 1 + monkeypatch.delenv("LNG_POOL_ACTIVE_VIA") + assert base_and_stats(*cdk2_shape())[0] == U["sx"] # the default is made_from + # the catalytic form reached through a set-pool node counts too + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "made_from") + net, rmap, umap, profiles, steps = cdk2_shape() + net.loc[(net["source_id"] == U["i2"]) & (net["edge_type"] == "catalyst"), "edge_type"] = "set_member" + net.loc[(net["source_id"] == U["i2"]) & (net["edge_type"] == "set_member"), "target_id"] = U["i4"] + net = pd.concat([net, pd.DataFrame([e("i4", "v8", "catalyst")])], ignore_index=True) + umap[U["i4"]] = "R-HSA-CDK2S" + profiles["R-HSA-CDK2S"] = prof([SUB, 41]) + assert base_and_stats(net, rmap, umap, profiles, steps)[0] == U["sx"] + + +def test_made_from_does_not_flip_the_ras_shape(monkeypatch): + # RAS:GTP binds RAF as an INPUT; RAS:GTP:RAF catalyses nothing within two reactions + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "made_from") + net, rmap, umap, profiles, participants = ras_network() + extra = [e("t", "v5", "input"), e("a2", "v5", "input"), e("v5", "i2", "output"), + e("i2", "v8", "input"), e("v8", "i3", "output")] + net = pd.concat([net, pd.DataFrame(extra)], ignore_index=True) + rmap = pd.concat([rmap, pd.DataFrame({"uid": [U["v5"], U["v8"]], "reactome_id": ["R-HSA-v5", "R-HSA-v8"]})], + ignore_index=True) + umap[U["a2"]], umap[U["i2"]], umap[U["i3"]] = "R-HSA-RAF", "R-HSA-RASRAF", "R-HSA-PRAF" + profiles.update({"R-HSA-RAF": prof([50]), "R-HSA-RASRAF": prof([RAS, 50], slots=2, comps=2), + "R-HSA-PRAF": prof([50], [("res", "p")], mods=1)}) + base, st = base_and_stats(net, rmap, umap, profiles, m.r_steps(participants, profiles), {"R-HSA-vgef"}) + assert base == U["g"] and st["activity_fallbacks"] == 1 + # a pool step's own reaction (the GAP binding) never counts as the route + assert st["oriented_by_activity"] == 0 + + +def test_unknown_pool_active_via_is_a_startup_error(monkeypatch): + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "indirect") + with pytest.raises(ValueError, match="LNG_POOL_ACTIVE_VIA"): + m._reject_removed_env() + with pytest.raises(ValueError): + m.find_pools(*two_state()[:3], two_state()[4], two_state()[3]) + monkeypatch.setenv("LNG_POOL_ACTIVE_VIA", "made_from") + m._reject_removed_env()