diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index f5ddae4..1851928 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -1675,6 +1675,12 @@ def _emit_boundary_decomposition_edges( ) -> None: """Expose the members of every root-input and terminal-output complex. + Under LNG_BOUNDARY_HIERARCHY=1 (the default since deltasignal specs/030) the + ROOT side is decomposed one hasComponent level at a time (see below); the + TERMINAL side is still flat. A produced species that the hierarchy joins + to a root keeps the dissociation sinks it was given as a terminal: they are + redundant readouts, harmless. + Boundary membership is decided **positionally, per network occurrence** — not by a global stId set difference. A node is a *root input* if it is a source but never a target (no reaction produces it); a *terminal output* if @@ -1828,12 +1834,99 @@ def _is_complex(entity_id: str) -> bool: # # Default ON pending the measurement that decides it. boundary_expansion = os.environ.get("LNG_BOUNDARY_EXPANSION", "1") == "1" + # LNG_BOUNDARY_HIERARCHY=1 (deltasignal specs/030): decompose a root complex + # ONE hasComponent level at a time instead of straight to its base leaves. + # A component that already exists as a node -- e.g. ISGF3 [cytosol], which + # the pathway produces -- is joined to the complex and not descended into; + # a nested complex that does not exist (ISGF3:KPNA1) is BUILT as a node and + # decomposed in turn; everything else falls back to its terminal leaves as + # before. Flat decomposition skipped every intermediate complex, so a + # species the pathway produces never reached the root complex that + # contains it: Reactome has no ISGF3 + KPNA1 + KPNB1 binding reaction, the + # root ISGF3:KPNA1:KPNB1 is translocated to the nucleus, and every IFN + # alpha/beta perturbation upstream of ISGF3 was severed there (200 held-out + # cases). The downstream-reuse rule (specs/018) applies unchanged. + hierarchy = os.environ.get("LNG_BOUNDARY_HIERARCHY", "1") == "1" # default since deltasignal specs/030 (+228 held-out) + from src.neo4j_connector import get_complex_components + nested_registry: Dict[tuple, str] = {} seen_edges: Set[tuple] = set() assembly_count = 0 - for complex_uuid in (root_uuids if boundary_expansion else ()): + nested_built = 0 + + # A ROOT copy (the one the benchmark's root-pinning protocol perturbs) is + # preferred over a produced copy of the same species. Joining the first + # eligible copy linked PDGF A/B heterodimer to a produced PDGFB copy while + # the pinned root copy fed only the processing reaction, so a PDGFB + # knockdown never reached the complex. "Root" is `targets` above: produced + # by no reaction of the pathway, depletion not counting as production. A + # live in-degree (review of PR #97) also counted the joins this function + # emits, so a root complex stopped being a root once it had been + # decomposed and the choice followed stid sort order; and it counted + # depletion edges, which root detection deliberately ignores. + def _existing_upstream(comp_stid: str, root_uuid: str): + downstream = _downstream_of(root_uuid) + ok = [c for c in stid_to_existing_uuids.get(comp_stid, []) + if c not in downstream and c != root_uuid] + if not ok: + return None + roots = [c for c in ok if c not in targets] + return (roots or ok)[0] + + def _emit(src: str, dst: str) -> None: + nonlocal assembly_count + if (src, dst) in seen_edges: + return + seen_edges.add((src, dst)) + pathway_logic_network_data.append({ + "source_id": src, "target_id": dst, "pos_neg": "pos", "and_or": "and", + "edge_type": "assembly", "stoichiometry": 1, + }) + assembly_count += 1 + if hierarchy: + # Keep the downstream test current: without this, two joins that + # each pass against the pre-loop snapshot can close a cycle + # together (review of PR #97: a 162-node SCC in DSB Repair). + _succ[src].append(dst) + _reach_cache.clear() + + def _decompose_hier(container_uuid: str, container_stid: str, root_uuid: str, depth: int) -> None: + nonlocal nested_built + for comp in sorted(get_complex_components(container_stid) or {}): + existing = _existing_upstream(comp, root_uuid) + if existing is not None: + _emit(existing, container_uuid) # a species the network already has + elif _is_complex(comp) and depth < 6: + # Built once PER ROOT: its inputs are chosen against that root's + # downstream set, so a copy built for another root may join a + # node that is not valid here, or miss one that is. + key = (root_uuid, comp) + if key not in nested_registry: + nested_registry[key] = str(uuid.uuid4()) + reactome_id_to_uuid[nested_registry[key]] = comp + nested_built += 1 + _decompose_hier(nested_registry[key], comp, root_uuid, depth + 1) + _emit(nested_registry[key], container_uuid) + else: + for leaf in sorted(get_terminal_components(comp)): + _emit(_existing_upstream(leaf, root_uuid) or _leaf_uuid(leaf, root_uuid), container_uuid) + + # Deterministic order under the hierarchy: with the downstream test now + # updated as edges are emitted, which of two jointly cycle-closing joins + # survives depends on order, and set order of uuid4 strings is not stable. + # Copies of one stid are ordered by insertion, NOT by uuid string: a uuid4 + # is new every run, so a string tiebreak picked a different copy run to + # run (review of PR #97; _emit_diagram_set_member_edges has the same rule). + _ordinal = {u: i for i, u in enumerate(reactome_id_to_uuid)} + root_order = (sorted(root_uuids, key=lambda u: (reactome_id_to_uuid.get(u) or "", + _ordinal.get(u, len(_ordinal)))) + if hierarchy else root_uuids) + for complex_uuid in (root_order if boundary_expansion else ()): stid = reactome_id_to_uuid.get(complex_uuid) or "" if not stid or not _is_complex(stid): continue + if hierarchy: + _decompose_hier(complex_uuid, stid, complex_uuid, 0) + continue leaves = get_terminal_components(stid) if leaves == {str(stid)}: # nothing below the complex to expose continue @@ -1875,6 +1968,8 @@ def _is_complex(entity_id: str) -> bool: }) dissociation_count += 1 + if hierarchy: + logger.info(f"Boundary hierarchy: {nested_built} nested complexes built") if assembly_count or dissociation_count: logger.info( f"Boundary expansion (positional): {assembly_count} assembly edges " diff --git a/src/pathway_generator.py b/src/pathway_generator.py index f7edf92..cfe7d4f 100755 --- a/src/pathway_generator.py +++ b/src/pathway_generator.py @@ -46,6 +46,7 @@ # which settings produced a catalog has been worth more than the re-fetch, # and src_sha256 already invalidates on any source change anyway. "LNG_BOUNDARY_EXPANSION", + "LNG_BOUNDARY_HIERARCHY", "LNG_COMPOSITION_EDGES", "LNG_EMIT_ONE_SIDED", # Determinism controls: node ids are uuid4 and several selections iterate diff --git a/tests/test_boundary_hierarchy.py b/tests/test_boundary_hierarchy.py new file mode 100644 index 0000000..b2bb9f3 --- /dev/null +++ b/tests/test_boundary_hierarchy.py @@ -0,0 +1,185 @@ +"""LNG_BOUNDARY_HIERARCHY=1 (deltasignal specs/030): a root complex is decomposed one +hasComponent level at a time. A component the network already has is JOINED; a nested +complex it lacks is BUILT; the rest fall back to terminal leaves. + +Shape: IFN alpha/beta. ISGF3 [cytosol] is produced upstream and consumed by nothing; the +root complex ISGF3:KPNA1:KPNB1 (produced by nothing) is translocated. Reactome has no +binding reaction between them, so flat decomposition (straight to STAT1/STAT2/IRF9/KPNA1/ +KPNB1 leaves) left every perturbation upstream of ISGF3 severed at that point.""" +import pytest +import src.neo4j_connector as nc +import src.logic_network_generator as lng +from src.logic_network_generator import _emit_boundary_decomposition_edges + +K, N, I = "R-HSA-K", "R-HSA-N", "R-HSA-I" # K = ISGF3:KPNA1:KPNB1, N = ISGF3:KPNA1, I = ISGF3 +A1, B1, S1, S2, IR = "R-HSA-KPNA1", "R-HSA-KPNB1", "R-HSA-STAT1", "R-HSA-STAT2", "R-HSA-IRF9" +LABELS = {K: ["Complex"], N: ["Complex"], I: ["Complex"], + **{x: ["EntityWithAccessionedSequence"] for x in (A1, B1, S1, S2, IR)}} +COMPONENTS = {K: {N: 1, B1: 1}, N: {I: 1, A1: 1}, I: {S1: 1, S2: 1, IR: 1}} +LEAVES = {K: {S1, S2, IR, A1, B1}, N: {S1, S2, IR, A1}, I: {S1, S2, IR}} + + +@pytest.fixture +def stub(monkeypatch): + monkeypatch.setattr(nc, "get_labels", lambda s: LABELS.get(s, [])) + monkeypatch.setattr(lng, "get_labels", lambda s: LABELS.get(s, []), raising=False) + monkeypatch.setattr(nc, "get_complex_components", lambda s: COMPONENTS.get(s, {})) + monkeypatch.setattr(lng, "get_terminal_components", lambda s: LEAVES.get(s, {s}), raising=False) + monkeypatch.delenv("LNG_COMPOSITION_EDGES", raising=False) + + +def network(isgf3_downstream_of_root=False): + # u_X -> r0 -> u_I (ISGF3 produced, dead end); u_K (root) -> r1 -> u_Z + data = [ + {"source_id": "u_X", "target_id": "r0", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r0", "target_id": "u_I", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + {"source_id": "u_K", "target_id": "r1", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r1", "target_id": "u_Z", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + ] + if isgf3_downstream_of_root: # ISGF3 is produced FROM the root complex: reusing it would weld a cycle + data[0] = {"source_id": "u_Z", "target_id": "r0", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1} + return data, {"u_X": "R-HSA-X", "u_I": I, "u_K": K, "u_Z": "R-HSA-Z"} + + +def assembly(data): + return [(e["source_id"], e["target_id"]) for e in data if e["edge_type"] == "assembly"] + + +def test_flat_mode_is_leaves_as_before(stub, monkeypatch): + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "0") + data, r2u = network() + _emit_boundary_decomposition_edges(data, r2u) + edges = assembly(data) + assert {t for _, t in edges} == {"u_K"} + assert {r2u[s] for s, _ in edges} == {S1, S2, IR, A1, B1} + assert ("u_I", "u_K") not in edges and all(s != "u_I" for s, _ in edges) + + +def test_hierarchy_builds_the_nested_complex_and_joins_the_produced_species(stub, monkeypatch): + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + data, r2u = network() + _emit_boundary_decomposition_edges(data, r2u) + edges = assembly(data) + nested = [u for u, s in r2u.items() if s == N] + assert len(nested) == 1 # ISGF3:KPNA1 built once + n = nested[0] + assert ("u_I", n) in edges # produced ISGF3 joins it ... + assert {r2u[s] for s, t in edges if t == n} == {I, A1} # ... with KPNA1 + assert {r2u[s] for s, t in edges if t == "u_K"} == {N, B1} # and it joins the root with KPNB1 + assert not any(t == "u_I" for _, t in edges) # ISGF3 is not descended into + + +def test_a_produced_copy_downstream_of_the_root_is_not_reused(stub, monkeypatch): + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + data, r2u = network(isgf3_downstream_of_root=True) + _emit_boundary_decomposition_edges(data, r2u) + edges = assembly(data) + assert all(s != "u_I" for s, _ in edges) # no weld through the root's own output + built_I = [u for u, s in r2u.items() if s == I and u != "u_I"] + assert len(built_I) == 1 # a separate ISGF3 node is built instead + assert {r2u[s] for s, t in edges if t == built_I[0]} == {S1, S2, IR} + + +def test_hierarchy_is_the_default(stub, monkeypatch): + monkeypatch.delenv("LNG_BOUNDARY_HIERARCHY", raising=False) + data, r2u = network() + _emit_boundary_decomposition_edges(data, r2u) + assert ("u_I", next(u for u, s in r2u.items() if s == N)) in assembly(data) + + +# --- review of PR #97 ------------------------------------------------------- + +def test_a_root_copy_is_preferred_over_a_produced_copy(stub, monkeypatch): + # PDGF shape: two copies of KPNB1. u_Bprod is produced by an unrelated + # reaction; u_Broot is a root (what the root-pinning benchmark perturbs). + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + data, r2u = network() + data += [{"source_id": "u_Y", "target_id": "r9", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r9", "target_id": "u_Bprod", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + {"source_id": "u_Broot", "target_id": "r8", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}] + r2u.update({"u_Y": "R-HSA-Y", "u_Bprod": B1, "u_Broot": B1}) + _emit_boundary_decomposition_edges(data, r2u) + into_root = {s for s, t in assembly(data) if t == "u_K"} + assert "u_Broot" in into_root and "u_Bprod" not in into_root + + +def test_two_joins_that_close_a_cycle_together_are_not_both_made(monkeypatch): + # Root R1 contains produced P2; root R2 contains produced P1. R1 -> rx -> P1 and + # R2 -> ry -> P2. Joining P2 -> R1 alone is acyclic, and so is P1 -> R2 alone, + # but both together close R1 -> P1 -> R2 -> P2 -> R1. A pre-loop snapshot + # allows both; the updated test must refuse the second. + R1, R2, P1, P2 = "R-HSA-R1", "R-HSA-R2", "R-HSA-P1", "R-HSA-P2" + labels = {R1: ["Complex"], R2: ["Complex"], P1: ["Complex"], P2: ["Complex"]} + comps = {R1: {P2: 1}, R2: {P1: 1}, P1: {}, P2: {}} + monkeypatch.setattr(nc, "get_labels", lambda s: labels.get(s, [])) + monkeypatch.setattr(lng, "get_labels", lambda s: labels.get(s, []), raising=False) + monkeypatch.setattr(nc, "get_complex_components", lambda s: comps.get(s, {})) + monkeypatch.setattr(lng, "get_terminal_components", lambda s: {s}, raising=False) + monkeypatch.delenv("LNG_COMPOSITION_EDGES", raising=False) + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + e = lambda a, b, t: {"source_id": a, "target_id": b, "pos_neg": "pos", "and_or": "and", "edge_type": t, "stoichiometry": 1} + data = [e("u_R1", "rx", "input"), e("rx", "u_P1", "output"), e("u_R2", "ry", "input"), e("ry", "u_P2", "output")] + r2u = {"u_R1": R1, "u_R2": R2, "u_P1": P1, "u_P2": P2} + _emit_boundary_decomposition_edges(data, r2u) + joins = {(s, t) for s, t in assembly(data)} + assert not {("u_P2", "u_R1"), ("u_P1", "u_R2")} <= joins # never both + + +def test_a_nested_complex_is_built_per_root(stub, monkeypatch): + # Two root occurrences of K each get their own ISGF3:KPNA1 node. + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + data, r2u = network() + data.append({"source_id": "u_K2", "target_id": "r2", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}) + r2u["u_K2"] = K + _emit_boundary_decomposition_edges(data, r2u) + assert len([u for u, s in r2u.items() if s == N]) == 2 + + +# --- second review of PR #97 ------------------------------------------------ + +def _fixture(monkeypatch, labels, comps): + monkeypatch.setattr(nc, "get_labels", lambda s: labels.get(s, [])) + monkeypatch.setattr(lng, "get_labels", lambda s: labels.get(s, []), raising=False) + monkeypatch.setattr(nc, "get_complex_components", lambda s: comps.get(s, {})) + monkeypatch.setattr(lng, "get_terminal_components", + lambda s: set(comps[s]) if comps.get(s) else {s}, raising=False) + monkeypatch.delenv("LNG_COMPOSITION_EDGES", raising=False) + monkeypatch.setenv("LNG_BOUNDARY_HIERARCHY", "1") + + +def _e(a, b, t): + return {"source_id": a, "target_id": b, "pos_neg": "pos", "and_or": "and", "edge_type": t, "stoichiometry": 1} + + +@pytest.mark.parametrize("k1, k2", [("u_z", "u_a"), ("u_a", "u_z")]) +def test_copies_of_one_root_are_ordered_by_insertion_not_uuid_string(monkeypatch, k1, k2): + # Two root copies of K = {A, B}: K1 -> r -> a_prod, K2 -> s -> b_prod. Whichever + # copy is processed first takes the real join (the other is then downstream of + # it), so the order must not follow the uuid STRING, which is new every run. + K_, A_, B_ = "R-HSA-K", "R-HSA-A", "R-HSA-B" + _fixture(monkeypatch, {K_: ["Complex"]}, {K_: {A_: 1, B_: 1}}) + data = [_e(k1, "r", "input"), _e("r", "a_prod", "output"), + _e(k2, "s", "input"), _e("s", "b_prod", "output")] + r2u = {k1: K_, k2: K_, "a_prod": A_, "b_prod": B_} # k1 inserted first + _emit_boundary_decomposition_edges(data, r2u) + joins = set(assembly(data)) + assert ("b_prod", k1) in joins and ("a_prod", k2) not in joins + + +@pytest.mark.parametrize("depleted", [False, True]) +def test_a_root_copy_stays_a_root_after_it_is_decomposed(monkeypatch, depleted): + # S has a root copy (itself a root complex, decomposed FIRST because its stid + # sorts first, which gives it incoming joins) and a produced copy. Root R + # contains S and must still prefer the root copy. A depletion edge into the + # root copy does not make it produced either. + S_, R_, L_ = "R-HSA-AAA", "R-HSA-ZZZ", "R-HSA-L" + _fixture(monkeypatch, {S_: ["Complex"], R_: ["Complex"]}, {S_: {L_: 1}, R_: {S_: 1}}) + data = [_e("s_root", "r1", "input"), _e("u_Y", "r2", "input"), _e("r2", "s_prod", "output"), + _e("u_R", "r3", "input")] + if depleted: + data.append(_e("u_W", "s_root", "depletion")) + # s_prod is registered first, so falling back to list order would pick it. + r2u = {"s_prod": S_, "s_root": S_, "u_R": R_, "u_Y": "R-HSA-Y", "u_W": "R-HSA-W"} + _emit_boundary_decomposition_edges(data, r2u) + into_r = {s for s, t in assembly(data) if t == "u_R"} + assert into_r == {"s_root"} diff --git a/tests/test_boundary_leaf_reuse.py b/tests/test_boundary_leaf_reuse.py index c081bbc..8e3cd37 100644 --- a/tests/test_boundary_leaf_reuse.py +++ b/tests/test_boundary_leaf_reuse.py @@ -20,6 +20,12 @@ def stub(monkeypatch): monkeypatch.setattr(lng, "get_labels", lambda stid: LABELS.get(stid, []), raising=False) monkeypatch.setattr(lng, "get_terminal_components", lambda stid: {P, Q} if stid == C else {stid}, raising=False) monkeypatch.delenv("LNG_COMPOSITION_EDGES", raising=False) + # One hasComponent level, derived from the same stub, so these reuse rules are + # exercised under LNG_BOUNDARY_HIERARCHY=1, the default since specs/030. + def _components(stid): + leaves = lng.get_terminal_components(stid) + return {} if leaves == {stid} else {x: 1 for x in leaves} + monkeypatch.setattr(nc, "get_complex_components", _components) def network(): diff --git a/tests/test_logic_network_generator.py b/tests/test_logic_network_generator.py index fd11cf5..aeb91e7 100644 --- a/tests/test_logic_network_generator.py +++ b/tests/test_logic_network_generator.py @@ -1,6 +1,7 @@ """Tests for logic_network_generator module.""" from typing import Dict, List, Any +import os import sys from pathlib import Path from unittest.mock import patch @@ -400,6 +401,12 @@ def test_non_boundary_entity_gets_separate_uuids(self): assert registry[("A", "vr1", "input")] != registry[("A", "vr2", "input")] +# Pinned to the flat leaf mode: this class stubs get_terminal_components but not +# get_complex_components, which the hierarchical default (specs/030) calls, so +# under the default the lookup would go to the database and return nothing for +# these fake ids. The leaf-reuse rules it pins are covered for the hierarchy in +# tests/test_boundary_hierarchy.py and tests/test_boundary_leaf_reuse.py. +@patch.dict(os.environ, {"LNG_BOUNDARY_HIERARCHY": "0"}) class TestBoundaryLeavesReuseExistingUUIDs: """Boundary expansion must NOT mint a fresh UUID for a leaf if that leaf's stId already has a UUID elsewhere in the network. Otherwise