diff --git a/src/logic_network_generator.py b/src/logic_network_generator.py index 693b85e..9656e0b 100755 --- a/src/logic_network_generator.py +++ b/src/logic_network_generator.py @@ -481,6 +481,78 @@ def _register_entity_uuid( return entity_uuid_registry[key] +def _register_phase1( + vr_entities: Dict[str, tuple], + entity_uuid_registry: Dict[tuple, str], + root_input_eids: Set[str], + root_input_uuid_cache: Dict[str, str], + terminal_output_eids: Set[str], + terminal_output_uuid_cache: Dict[str, str], + vr_to_reaction: Optional[Dict[str, str]] = None, + share_variants: bool = False, +) -> Dict[str, int]: + """Phase 1 of UUID assignment: one key per (entity, virtual reaction, role). + + A Reactome reaction with EntitySet participants becomes several VIRTUAL + reactions, one per member combination, and every one of them registers its + own UUID for every participant. An entity common to all variants -- BCDX2 in + HDR's strand-invasion reactions -- therefore exists as one node PER VARIANT: + 33 copies from 5 reactions, each a full copy with the same downstream edges. + Nothing downstream can tell them apart, and any solver rule over the copies + is wrong in one direction or the other: min over copies caps the container + by whichever copy a perturbation missed (missed change); max lets one + elevated copy lift it (false change). Both were measured (deltasignal + specs/016): -83 and -117 held-out. + + `share_variants` collapses exactly that: the SAME entity in the SAME role + across the variants of ONE Reactome reaction shares one UUID. Variants that + differ in WHICH set member they use have different entity ids and are not + conflated. Copies across DIFFERENT reactions stay positional, as designed + (bridging those was measured harmful four times). Boundary entities already + share per stId through their caches and are left to them. + + Returns a small stats dict for the log. + """ + vr_to_reaction = vr_to_reaction or {} + variant_cache: Dict[tuple, str] = {} + shared = 0 + unmapped = 0 + for vr_uid, (input_ids, output_ids, *_) in vr_entities.items(): + rxn = vr_to_reaction.get(str(vr_uid)) + if share_variants and rxn in (None, "", "nan", "None"): + # `reactome_id` read back from a cached CSV with no dtype turns a + # missing id into the STRING "nan"; treated as a real key it would + # collapse every such variant across reactions onto one uuid. + unmapped += 1 + for role, ids, boundary_eids, boundary_cache in ( + ("input", input_ids, root_input_eids, root_input_uuid_cache), + ("output", output_ids, terminal_output_eids, terminal_output_uuid_cache), + ): + for eid in ids: + key = (eid, vr_uid, role) + if (share_variants and rxn not in (None, "", "nan", "None") + and eid not in boundary_eids + and key not in entity_uuid_registry): + vkey = (eid, rxn, role) + if vkey in variant_cache: + entity_uuid_registry[key] = variant_cache[vkey] + shared += 1 + continue + u = _register_entity_uuid(eid, vr_uid, role, entity_uuid_registry, + boundary_eids, boundary_cache) + variant_cache[vkey] = u + continue + _register_entity_uuid(eid, vr_uid, role, entity_uuid_registry, + boundary_eids, boundary_cache) + if share_variants: + logger.info(f"Variant-node sharing: {shared} (entity, reaction, role) registrations " + f"reused a sibling variant's UUID ({len(variant_cache)} distinct)") + if unmapped: + logger.warning(f"Variant-node sharing: {unmapped} virtual reactions have no Reactome " + f"reaction id and were registered per-variant (sharing is inconsistent for them)") + return {"shared": shared, "distinct": len(variant_cache), "unmapped": unmapped} + + def _build_entity_producer_count(vr_entities: Dict[str, tuple]) -> Dict[str, int]: """Count how many VRs produce each entity as output. @@ -1468,6 +1540,93 @@ def _emit_diagram_set_member_edges( return emitted +def _emit_composition_edges( + pathway_logic_network_data: List[Dict[str, Any]], + reactome_id_to_uuid: Dict[str, str], +) -> None: + """Emit ``composition`` edges: complex node -> node of a complex that CONTAINS it. + + Reactome joins some entities by composition alone. Cytosolic ISGF3 is + consumed by no reaction; it is a component of ISGF3:KPNA1, which is a + component of ISGF3:KPNA1:KPNB1, and only then does translocation appear. A + curator crosses that in their head. A reaction-only network cannot, and the + result on Interferon alpha/beta is that the live ISGF3 node's only outgoing + edges are four dissociation sinks: the entire nuclear branch is severed at + one node, and IFNAR2 or JAK1 knockouts reach none of the ISGF3 targets. + + This is NOT the leaf-subunit bridge that lost four times (-15pp; silo + bridges -73/-77). Those reconnected a released SUBUNIT to its functional + node, whose fan-out is ~112: a broadcast. This connects a COMPLEX to the + complexes it sits inside, and a complex sits inside few complexes + (measured median 1, max 2). Two hops, because an intermediate complex may + exist only as a set member and have no node (ISGF3:KPNA1 has none). + + Simulated at uuid level before this was written: severed curator routes + 392 -> 43 in Interferon, 462 -> 165 in Mitotic G1 (specs/016 in + deltasignal). Reachability is necessary, not sufficient -- the held-out A/B + decides it, on ONE catalog via DS_SKIP_EDGE_TYPES=composition. + + Off by default (LNG_COMPOSITION_EDGES=1 enables). Dissociation sinks are + never sources: they are readout handles, not species that flow. + """ + from collections import defaultdict + from src.neo4j_connector import get_labels, get_containing_complexes + + sink_uuids: Set[str] = set() + existing: Set[tuple] = set() + for e in pathway_logic_network_data: + s_, t_ = str(e.get("source_id")), str(e.get("target_id")) + existing.add((s_, t_)) + if e.get("edge_type") == "dissociation": + sink_uuids.add(t_) + + base_to_uuids: Dict[str, List[str]] = defaultdict(list) + for u, stid in reactome_id_to_uuid.items(): + if u in sink_uuids or not str(stid).startswith("R-"): + continue + base_to_uuids[str(stid).split("::variant::")[0]].append(str(u)) + + n_edges = 0 + fan: List[int] = [] + for base, uuids in base_to_uuids.items(): + try: + labels = get_labels(base) + except (IndexError, KeyError): + # unknown stId (get_labels does .data()[0]); a connection error must + # NOT be swallowed here -- it would silently drop this complex's + # edges while the next lookup aborts the pathway anyway. + labels = [] + if "Complex" not in labels: + continue + containers = get_containing_complexes(base, 2) + targets = [(y, uy) for y in containers for uy in base_to_uuids.get(y, [])] + if not targets: + continue + fan.append(len({y for y, _ in targets})) + for ux in uuids: + for _y, uy in targets: + if ux == uy or (ux, uy) in existing: + continue + existing.add((ux, uy)) + pathway_logic_network_data.append({ + "source_id": ux, + "target_id": uy, + "pos_neg": "pos", + "and_or": "and", + "edge_type": "composition", + "stoichiometry": 1, + }) + n_edges += 1 + if fan: + fan_sorted = sorted(fan) + logger.info( + f"Composition edges: {n_edges} edges from {len(fan)} complexes; " + f"containing-complex fan-out median {fan_sorted[len(fan_sorted)//2]}, max {fan_sorted[-1]}" + ) + else: + logger.info("Composition edges: none (no complex node sits inside another complex node here)") + + def _emit_boundary_decomposition_edges( pathway_logic_network_data: List[Dict[str, Any]], reactome_id_to_uuid: Dict[str, str], @@ -1541,16 +1700,53 @@ def _emit_boundary_decomposition_edges( # stId → existing UUID, so a member reuses the node it already has elsewhere # (free protein, regulator, catalyst) rather than becoming a disconnected dup. - stid_to_existing_uuid: Dict[str, str] = {} + # A root complex's subunit leaf is a boundary INPUT. Reusing an existing + # node for it is right whenever that node is NOT downstream of the root + # complex: a free protein that is itself a root, a catalyst, a regulator, + # or a copy produced by an unrelated reaction (that last case is a real + # feed-forward link Reactome joins only by hasComponent -- severing it cost + # Mitotic G1 28 cases). Reusing a node the root complex can REACH welds a + # cycle Neo4j never had: root complex -> reactions -> ... -> produced + # protein -> (assembly) -> root complex. On the v97 catalog 1,994 of the + # 2,077 cycle-carrying assembly edges were exactly that shape; removing them + # takes TP53's strongly connected component from 836 nodes to 56 and DSB's + # from 1,127 to ~290 (deltasignal specs/018). Reachability is computed over + # every edge emitted so far (bridges and depletion edges included). + # LNG_BOUNDARY_LEAF_REUSE=any restores the old behaviour. + from collections import defaultdict + reuse_mode = os.environ.get("LNG_BOUNDARY_LEAF_REUSE", "downstream_free") + if reuse_mode not in ("downstream_free", "any"): + raise ValueError(f"LNG_BOUNDARY_LEAF_REUSE={reuse_mode!r}: expected 'downstream_free' or 'any'") + _succ: Dict[str, List[str]] = defaultdict(list) + for e in pathway_logic_network_data: + _succ[str(e.get("source_id"))].append(str(e.get("target_id"))) + _reach_cache: Dict[str, Set[str]] = {} + def _downstream_of(root_uuid: str) -> Set[str]: + if root_uuid not in _reach_cache: + seen: Set[str] = set(); stack = [root_uuid] + while stack: + u = stack.pop() + for v in _succ.get(u, ()): + if v not in seen: + seen.add(v); stack.append(v) + _reach_cache[root_uuid] = seen + return _reach_cache[root_uuid] + # stId -> every existing node carrying it, in registry order (deterministic) + stid_to_existing_uuids: Dict[str, List[str]] = defaultdict(list) for existing_uuid, stid in reactome_id_to_uuid.items(): - if stid not in stid_to_existing_uuid: - stid_to_existing_uuid[stid] = existing_uuid - + stid_to_existing_uuids[stid].append(existing_uuid) leaf_uuid_registry: Dict[str, str] = {} - def _leaf_uuid(leaf_stid: str) -> str: - if leaf_stid in stid_to_existing_uuid: - return stid_to_existing_uuid[leaf_stid] + def _leaf_uuid(leaf_stid: str, root_uuid: str) -> str: + candidates = stid_to_existing_uuids.get(leaf_stid, []) + if reuse_mode == "any": + if candidates: + return candidates[0] + else: + downstream = _downstream_of(root_uuid) + for cand in candidates: + if cand not in downstream: + return cand if leaf_stid not in leaf_uuid_registry: leaf_uuid_registry[leaf_stid] = str(uuid.uuid4()) reactome_id_to_uuid[leaf_uuid_registry[leaf_stid]] = leaf_stid @@ -1593,7 +1789,7 @@ def _is_complex(entity_id: str) -> bool: if leaves == {str(stid)}: # nothing below the complex to expose continue for leaf in leaves: - leaf_uuid = _leaf_uuid(leaf) + leaf_uuid = _leaf_uuid(leaf, complex_uuid) if (leaf_uuid, complex_uuid) in seen_edges: continue seen_edges.add((leaf_uuid, complex_uuid)) @@ -1637,6 +1833,11 @@ def _is_complex(entity_id: str) -> bool: f"(separate readout sinks), {len(leaf_uuid_registry)} new assembly leaves" ) + # Composition hierarchy (LNG_COMPOSITION_EDGES). Runs after boundary + # expansion so dissociation sinks are known and excluded as sources. + if os.environ.get("LNG_COMPOSITION_EDGES", "0") == "1": + _emit_composition_edges(pathway_logic_network_data, reactome_id_to_uuid) + def append_regulators( catalyst_map: pd.DataFrame, @@ -1946,13 +2147,16 @@ def create_pathway_logic_network( # Each entity gets a unique UUID per (entity, reaction, role) triple. # No cross-role keys are created (unlike the old self-loop approach). # Boundary entities (root inputs / terminal outputs) share one UUID per stId. - for vr_uid, (input_ids, output_ids, *_) in vr_entities.items(): - for eid in input_ids: - _register_entity_uuid(eid, vr_uid, "input", entity_uuid_registry, - root_input_eids, root_input_uuid_cache) - for eid in output_ids: - _register_entity_uuid(eid, vr_uid, "output", entity_uuid_registry, - terminal_output_eids, terminal_output_uuid_cache) + # Under LNG_SHARE_VARIANT_NODES the same entity in the same role across the + # VARIANTS of one Reactome reaction shares one UUID (see _register_phase1). + _register_phase1( + vr_entities, entity_uuid_registry, + root_input_eids, root_input_uuid_cache, + terminal_output_eids, terminal_output_uuid_cache, + vr_to_reaction=dict(zip(reaction_id_map["uid"].astype(str), + reaction_id_map["reactome_id"].astype(str))), + share_variants=os.environ.get("LNG_SHARE_VARIANT_NODES", "0") == "1", + ) logger.debug(f"Phase 1 complete: {len(entity_uuid_registry)} registry entries") diff --git a/src/neo4j_connector.py b/src/neo4j_connector.py index 69f19e4..3c08384 100755 --- a/src/neo4j_connector.py +++ b/src/neo4j_connector.py @@ -588,6 +588,44 @@ def get_complex_components(entity_id: str) -> Dict[str, int]: raise +_containing_cache: Dict[str, Dict[str, int]] = {} + + +def get_containing_complexes(entity_id: str, max_hops: int = 2) -> Dict[str, int]: + """Complexes that CONTAIN `entity_id`, within `max_hops` hasComponent steps. + + The upward counterpart of :func:`get_complex_components`. Returns + ``{containing_stId: hops}``. Used to emit ``composition`` edges from a + complex node to the complexes it is part of -- a route a curator follows + (ISGF3 -> ISGF3:KPNA1 -> ISGF3:KPNA1:KPNB1) that a reaction-only network + cannot. Two hops, not one, because an intermediate complex may exist only + as a set member and have no node of its own. Fan-out is inherently small + (a complex sits inside few complexes: measured median 1, max 2), which is + what distinguishes this from leaf-subunit bridging. + """ + # `entity` carries the PhysicalEntity label so the planner can use the stId + # index: unlabelled, this was a full-store scan (660-800 ms per call, once + # per Complex per pathway -- 30-90 min across the catalog); labelled, 1-21 ms. + key = f"{entity_id}|{max_hops}" + if key in _containing_cache: + return _containing_cache[key] + try: + data = get_graph().run( + """ + MATCH path=(container:Complex)-[:hasComponent*1..%d]->(entity:PhysicalEntity) + WHERE entity.stId = $entity_id AND container.stId IS NOT NULL + RETURN container.stId AS container_id, min(length(path)) AS hops + """ % max_hops, + entity_id=entity_id, + ).data() + result = {row["container_id"]: int(row["hops"]) for row in data} + _containing_cache[key] = result + return result + except Exception: + logger.error("Error in get_containing_complexes", **_traceback_kwargs()) + raise + + def get_set_members(entity_id: str) -> Set[str]: if entity_id in _members_cache: return _members_cache[entity_id] diff --git a/tests/test_boundary_leaf_reuse.py b/tests/test_boundary_leaf_reuse.py new file mode 100644 index 0000000..47a64cd --- /dev/null +++ b/tests/test_boundary_leaf_reuse.py @@ -0,0 +1,133 @@ +"""A root complex's subunit leaf must not reuse a node that is DOWNSTREAM of that complex. + +Reusing it welds root complex -> reactions -> produced protein -> (assembly) -> root +complex, a cycle Neo4j never had (1,994 of 2,077 cycle-carrying assembly edges on the +v97 catalog; deltasignal specs/018). A produced node that is NOT downstream is a real +feed-forward link (hasComponent with no reaction) and must still be reused -- severing +those cost Mitotic G1 28 cases in the first version of this fix.""" +import os +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 + +C, P, Q = "R-HSA-100", "R-HSA-200", "R-HSA-300" # C: root complex of P and Q +LABELS = {C: ["Complex"], P: ["EntityWithAccessionedSequence"], Q: ["EntityWithAccessionedSequence"]} + + +@pytest.fixture +def stub(monkeypatch): + monkeypatch.setattr(nc, "get_labels", lambda stid: LABELS.get(stid, [])) + 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) + + +def network(): + # C (root complex) -> r1 -> P_out (a PRODUCED copy of P); Q_free is a root protein with no producer + data = [ + {"source_id": "u_C", "target_id": "r1", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r1", "target_id": "u_Pout", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + {"source_id": "u_Qfree", "target_id": "r2", "pos_neg": "pos", "and_or": "and", "edge_type": "catalyst", "stoichiometry": 1}, + {"source_id": "r2", "target_id": "u_X", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + ] + r2u = {"u_C": C, "u_Pout": P, "u_Qfree": Q, "u_X": "R-HSA-400"} + return data, r2u + + +def assembly_sources(data): + return {e["source_id"]: e for e in data if e["edge_type"] == "assembly"} + + +def test_produced_node_is_not_reused_as_a_boundary_leaf(stub, monkeypatch): + monkeypatch.delenv("LNG_BOUNDARY_LEAF_REUSE", raising=False) + data, r2u = network() + _emit_boundary_decomposition_edges(data, r2u) + asm = assembly_sources(data) + targets = {e["target_id"] for e in asm.values()} + assert targets == {"u_C"} + # P's leaf is a FRESH node (u_Pout is produced by r1), Q's leaf reuses the unproduced free node + assert "u_Pout" not in asm + assert "u_Qfree" in asm + fresh = [u for u in asm if u not in ("u_Qfree",)] + assert len(fresh) == 1 and r2u[fresh[0]] == P + # and therefore no cycle: nothing produced feeds back into the root complex + sources_into_C = {e["source_id"] for e in data if e["target_id"] == "u_C"} + produced = {e["target_id"] for e in data if e["edge_type"] == "output"} + assert not (sources_into_C & produced) + + +def test_a_copy_reached_only_through_a_bridge_is_downstream_too(stub, monkeypatch): + # P has two copies: u_Pout (reaction output, downstream of C) and u_Pin (an input copy whose + # only incoming edge is a diagram bridge from u_Pout -- also downstream). Neither may be reused. + monkeypatch.delenv("LNG_BOUNDARY_LEAF_REUSE", raising=False) + data, r2u = network() + data.append({"source_id": "u_Pout", "target_id": "u_Pin", "pos_neg": "pos", "and_or": "and", "edge_type": "diagram_bridge", "stoichiometry": 1}) + data.append({"source_id": "u_Pin", "target_id": "r3", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}) + data.append({"source_id": "r3", "target_id": "u_Y", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}) + r2u.update({"u_Pin": P, "u_Y": "R-HSA-500"}) + _emit_boundary_decomposition_edges(data, r2u) + asm = assembly_sources(data) + assert "u_Pout" not in asm and "u_Pin" not in asm + + +def test_a_produced_copy_that_is_not_downstream_is_reused(stub, monkeypatch): + # The Mitotic G1 shape: P's only copy is produced by a reaction the root complex C does NOT + # reach (r2, fed by Q). Reusing it adds no cycle and keeps the feed-forward link + # P -> (assembly) -> C that Reactome joins only by hasComponent. + monkeypatch.delenv("LNG_BOUNDARY_LEAF_REUSE", raising=False) + data = [ + {"source_id": "u_C", "target_id": "r1", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r1", "target_id": "u_W", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + {"source_id": "u_Qfree", "target_id": "r2", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}, + {"source_id": "r2", "target_id": "u_Pelse", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}, + ] + r2u = {"u_C": C, "u_W": "R-HSA-700", "u_Qfree": Q, "u_Pelse": P} + _emit_boundary_decomposition_edges(data, r2u) + asm = assembly_sources(data) + assert "u_Pelse" in asm and asm["u_Pelse"]["target_id"] == "u_C" # reused: a real link, no cycle + assert "u_Qfree" in asm + assert sum(1 for e in data if e["edge_type"] == "assembly") == 2 + + +def test_two_root_complexes_sharing_a_subunit_each_get_an_acyclic_leaf(stub, monkeypatch): + # C reaches P's produced copy (u_Pout) -> C gets a fresh leaf. C2 does NOT reach it -> + # C2 reuses u_Pout (a real feed-forward link). Neither assembly edge closes a cycle. + monkeypatch.delenv("LNG_BOUNDARY_LEAF_REUSE", raising=False) + data, r2u = network() + data.append({"source_id": "u_C2", "target_id": "r4", "pos_neg": "pos", "and_or": "and", "edge_type": "input", "stoichiometry": 1}) + data.append({"source_id": "r4", "target_id": "u_Z", "pos_neg": "pos", "and_or": "or", "edge_type": "output", "stoichiometry": 1}) + r2u.update({"u_C2": "R-HSA-101", "u_Z": "R-HSA-600"}) + LABELS["R-HSA-101"] = ["Complex"] + monkeypatch.setattr(lng, "get_terminal_components", lambda stid: {P, Q} if stid in (C, "R-HSA-101") else {stid}, raising=False) + _emit_boundary_decomposition_edges(data, r2u) + into = {t_: {e["source_id"] for e in data if e["edge_type"] == "assembly" and e["target_id"] == t_} for t_ in ("u_C", "u_C2")} + assert "u_Pout" not in into["u_C"] and "u_Pout" in into["u_C2"] + assert any(r2u.get(u) == P for u in into["u_C"]) + # acyclicity: no assembly source into a complex is reachable from that complex + succ = {} + for e in data: succ.setdefault(e["source_id"], []).append(e["target_id"]) + def reach(u): + seen, st = set(), [u] + while st: + x = st.pop() + for v in succ.get(x, []): + if v not in seen: seen.add(v); st.append(v) + return seen + for cpx, srcs in into.items(): + assert not (srcs & reach(cpx)) + + +def test_legacy_reuse_is_available_and_welds_the_cycle(stub, monkeypatch): + monkeypatch.setenv("LNG_BOUNDARY_LEAF_REUSE", "any") + data, r2u = network() + _emit_boundary_decomposition_edges(data, r2u) + asm = assembly_sources(data) + assert "u_Pout" in asm # the old behaviour: the produced copy is reused -> C -> r1 -> u_Pout -> C + + +def test_bad_mode_is_an_error(stub, monkeypatch): + monkeypatch.setenv("LNG_BOUNDARY_LEAF_REUSE", "unproduced") + data, r2u = network() + with pytest.raises(ValueError): + _emit_boundary_decomposition_edges(data, r2u) diff --git a/tests/test_composition_edges.py b/tests/test_composition_edges.py new file mode 100644 index 0000000..c98046d --- /dev/null +++ b/tests/test_composition_edges.py @@ -0,0 +1,104 @@ +"""Unit tests for ``_emit_composition_edges`` (LNG_COMPOSITION_EDGES). + +A complex node gets an edge to the node of every complex that CONTAINS it, +within two hasComponent hops. This is the non-broadcasting repair for the +composition-only gap: it connects a complex to the few complexes it sits +inside (fan-out median 1, max 2), never a released subunit to its hub. Neo4j is +stubbed; the emission logic is what is under test. +""" +import os + +import pytest + +import src.neo4j_connector as nc +from src.logic_network_generator import _emit_composition_edges + +# stable ids in a toy pathway +ISGF3, ISGF3_KPNA1, IMPORTIN = "R-HSA-909698", "R-HSA-9710958", "R-HSA-9710965" +STAT1 = "R-HSA-873791" # a protein (EWAS), never a composition source +LABELS = {ISGF3: ["Complex"], ISGF3_KPNA1: ["Complex"], IMPORTIN: ["Complex"], STAT1: ["EntityWithAccessionedSequence"]} +# containment: importin contains ISGF3:KPNA1 (1 hop) which contains ISGF3, so importin contains ISGF3 at 2 hops +CONTAINERS = {ISGF3: {ISGF3_KPNA1: 1, IMPORTIN: 2}, ISGF3_KPNA1: {IMPORTIN: 1}, IMPORTIN: {}, STAT1: {ISGF3: 1}} + + +@pytest.fixture +def stub_neo4j(monkeypatch): + monkeypatch.setattr(nc, "get_labels", lambda stid: LABELS.get(stid, [])) + monkeypatch.setattr(nc, "get_containing_complexes", lambda stid, max_hops=2: CONTAINERS.get(stid, {})) + + +def edge(src, tgt, ty="input"): + return {"source_id": src, "target_id": tgt, "pos_neg": "pos", "and_or": "and", "edge_type": ty, "stoichiometry": 1} + + +def composition_edges(data): + return {(e["source_id"], e["target_id"]) for e in data if e["edge_type"] == "composition"} + + +def test_complex_connects_to_the_complexes_that_contain_it(stub_neo4j): + # nodes: ISGF3 (u1), importin (u3). ISGF3:KPNA1 has NO node (set member) -- the 2-hop case. + uuid_to_stid = {"u1": ISGF3, "u3": IMPORTIN} + data = [edge("rX", "u1", "output")] + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == {("u1", "u3")} + new = [e for e in data if e["edge_type"] == "composition"][0] + assert (new["pos_neg"], new["and_or"], new["stoichiometry"]) == ("pos", "and", 1) + + +def test_one_hop_and_two_hop_both_emitted_when_both_have_nodes(stub_neo4j): + uuid_to_stid = {"u1": ISGF3, "u2": ISGF3_KPNA1, "u3": IMPORTIN} + data = [] + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == {("u1", "u2"), ("u1", "u3"), ("u2", "u3")} + + +def test_a_dissociation_sink_is_never_a_source(stub_neo4j): + # u_sink carries ISGF3's stable id but is a readout sink (target of a dissociation edge). + uuid_to_stid = {"u1": ISGF3, "u_sink": ISGF3, "u3": IMPORTIN} + data = [edge("cplx", "u_sink", "dissociation")] + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == {("u1", "u3")} + assert all(e["source_id"] != "u_sink" for e in data if e["edge_type"] == "composition") + + +def test_proteins_are_not_sources_even_if_something_contains_them(stub_neo4j): + # STAT1 is contained by ISGF3, but STAT1 is a protein: that is the leaf-subunit + # bridge that lost four times, and it must NOT be emitted here. + uuid_to_stid = {"s1": STAT1, "u1": ISGF3} + data = [] + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == set() + + +def test_no_self_edge_and_no_duplicate_of_an_existing_edge(stub_neo4j): + uuid_to_stid = {"u1": ISGF3, "u3": IMPORTIN} + data = [edge("u1", "u3", "assembly")] # an edge already exists u1 -> u3 + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == set() # deduped against the existing pair + assert len(data) == 1 + + +def test_container_without_a_node_yields_nothing(stub_neo4j): + uuid_to_stid = {"u1": ISGF3} # neither container is a node here + data = [] + _emit_composition_edges(data, uuid_to_stid) + assert data == [] + + +def test_variant_nodes_resolve_to_their_base_complex(stub_neo4j): + uuid_to_stid = {"v1": f"{ISGF3}::variant::R-HSA-1", "v2": f"{ISGF3}::variant::R-HSA-2", "u3": IMPORTIN} + data = [] + _emit_composition_edges(data, uuid_to_stid) + assert composition_edges(data) == {("v1", "u3"), ("v2", "u3")} + + +def test_flag_default_is_off(monkeypatch): + # The emitter is only called from _emit_boundary_decomposition_edges under + # the flag; with the flag unset the call site must not reach it. + monkeypatch.delenv("LNG_COMPOSITION_EDGES", raising=False) + import src.logic_network_generator as lng + calls = [] + monkeypatch.setattr(lng, "_emit_composition_edges", lambda *a, **k: calls.append(1)) + src = open(lng.__file__).read() + assert 'os.environ.get("LNG_COMPOSITION_EDGES", "0") == "1"' in src + assert calls == [] diff --git a/tests/test_connectivity_ladder.py b/tests/test_connectivity_ladder.py index 63e0b3a..e5b6f76 100644 --- a/tests/test_connectivity_ladder.py +++ b/tests/test_connectivity_ladder.py @@ -260,6 +260,10 @@ def forward(self): def uuids(self): return _stid_to_uuids(BUNDLE) + @pytest.fixture(scope="class") + def network(self): + return pd.read_csv(BUNDLE / "logic_network.csv", dtype=str) + def test_every_entity_in_the_chain_is_a_node_here(self, uuids): """Guards the xfail below. @@ -290,14 +294,34 @@ def test_nuclear_isgf3_reaches_the_readout(self, forward, uuids): "the break has moved downstream of nuclear import." ) - @pytest.mark.xfail( - strict=True, - reason="The cytosolic half is severed at the boundary bridge, so no " - "perturbation upstream of ISGF3 reaches the readout. This is the " - "end-to-end symptom of the Tier 2 and Tier 3 failures.", - ) - def test_cytosolic_isgf3_reaches_the_readout(self, forward, uuids): - assert self._reaches(forward, uuids[ISGF3_CYTOSOL], uuids[EXPRESSION]) + def test_cytosolic_isgf3_reaches_the_readout(self, forward, uuids, network): + """The end-to-end symptom, and the fix that clears it. + + With a reaction-only network the cytosolic half is severed at the + boundary bridge and no perturbation upstream of ISGF3 reaches the + readout. `LNG_COMPOSITION_EDGES=1` emits complex -> containing-complex + edges along Reactome's hasComponent hierarchy (ISGF3 -> ISGF3:KPNA1 -> + ISGF3:KPNA1:KPNB1); 15 such edges reconnect this branch, and the oracle + in deltasignal specs/016 measured severed curator routes 392 -> 52. + + So the expectation depends on how the bundle was built: a bundle with + no composition edges is EXPECTED to fail here (that is the defect); + one with them must pass. Tiers 2-3 stay xfail -- they describe the + leaf-subunit bridge, which is the repair measured to broadcast and is + deliberately not made. + """ + has_composition = (network["edge_type"] == "composition").any() + reaches = self._reaches(forward, uuids[ISGF3_CYTOSOL], uuids[EXPRESSION]) + if not has_composition: + assert not reaches, ( + "A bundle WITHOUT composition edges reached the readout: either " + "the severing was repaired some other way (record how) or the " + "detector is wrong." + ) + pytest.xfail("no composition edges in this bundle; regenerate with " + "LNG_COMPOSITION_EDGES=1 to exercise the repair") + assert reaches, ("composition edges are present but cytosolic ISGF3 still does " + "not reach the readout -- the repair did not join the chain") # --- the invariant itself, on synthetic fixtures ---------------------------- diff --git a/tests/test_variant_node_sharing.py b/tests/test_variant_node_sharing.py new file mode 100644 index 0000000..932bbc9 --- /dev/null +++ b/tests/test_variant_node_sharing.py @@ -0,0 +1,91 @@ +"""LNG_SHARE_VARIANT_NODES: the same entity in the same role across the variants of +ONE Reactome reaction shares one UUID; across different reactions it does not. + +Motivation (deltasignal specs/016): HDR's BCDX2 complex existed as 33 node copies +-- one per variant reaction -- and no solver rule over copies is right: min over +them measured -83 held-out, max -117. One node is the fix. +""" +from src.logic_network_generator import _register_phase1 + + +def phase1(vr_entities, vr_to_reaction, share, roots=frozenset(), terms=frozenset()): + reg = {} + _register_phase1(vr_entities, reg, set(roots), {}, set(terms), {}, + vr_to_reaction=vr_to_reaction, share_variants=share) + return reg + + +# Reaction R1 has two variants (v1, v2) that both produce BCDX2 ("B") and consume "A"; +# reaction R2 has one variant (v3) that also produces "B". +VR = {"v1": (["A", "S1"], ["B", "P1"]), "v2": (["A", "S2"], ["B", "P2"]), "v3": (["C"], ["B"])} +V2R = {"v1": "R1", "v2": "R1", "v3": "R2"} + + +def test_off_every_variant_gets_its_own_uuid(): + reg = phase1(VR, V2R, share=False) + assert reg[("B", "v1", "output")] != reg[("B", "v2", "output")] + assert reg[("A", "v1", "input")] != reg[("A", "v2", "input")] + + +def test_on_same_entity_same_reaction_same_role_shares_one_uuid(): + reg = phase1(VR, V2R, share=True) + assert reg[("B", "v1", "output")] == reg[("B", "v2", "output")] + assert reg[("A", "v1", "input")] == reg[("A", "v2", "input")] + + +def test_on_different_reactions_stay_distinct(): + reg = phase1(VR, V2R, share=True) + assert reg[("B", "v1", "output")] != reg[("B", "v3", "output")] # positional across reactions, as designed + + +def test_on_roles_are_not_conflated(): + vr = {"v1": (["B"], ["B"])} # B both consumed and produced by one variant + reg = phase1(vr, {"v1": "R1"}, share=True) + assert reg[("B", "v1", "input")] != reg[("B", "v1", "output")] + + +def test_on_variant_specific_members_are_not_conflated(): + reg = phase1(VR, V2R, share=True) + assert reg[("S1", "v1", "input")] != reg[("S2", "v2", "input")] # different entities + assert reg[("P1", "v1", "output")] != reg[("P2", "v2", "output")] + + +def test_boundary_entities_are_left_to_their_own_caches(): + # "A" is a root input: the boundary cache already shares it per stId; sharing + # must not create a second uuid for it. + reg = phase1(VR, V2R, share=True, roots={"A"}) + assert reg[("A", "v1", "input")] == reg[("A", "v2", "input")] + assert len({reg[("A", v, "input")] for v in ("v1", "v2")}) == 1 + + +def test_terminal_output_shared_across_reactions_by_its_own_cache_not_the_variant_cache(): + # "B" is an output in R1 (v1, v2) AND in R2 (v3), and a terminal output. The + # boundary cache shares it per stId across ALL three; an implementation that + # ignored the boundary exclusion and shared via the variant cache would give + # v3 a different uuid from v1/v2. All three must be one uuid. + reg = phase1(VR, V2R, share=True, terms={"B"}) + assert len({reg[("B", v, "output")] for v in ("v1", "v2", "v3")}) == 1 + + +def test_missing_reaction_id_read_back_as_nan_string_is_not_a_shared_key(): + # reactome_id read from a cached CSV with no dtype turns NaN into "nan"; + # every variant of every reaction with a missing id must NOT collapse onto + # one uuid through the key (eid, "nan", role). + reg = phase1(VR, {"v1": "nan", "v2": "nan", "v3": "nan"}, share=True) + assert reg[("B", "v1", "output")] != reg[("B", "v2", "output")] + assert reg[("B", "v1", "output")] != reg[("B", "v3", "output")] + + +def test_partial_reaction_map_is_reported(): + from src.logic_network_generator import _register_phase1 + reg = {} + stats = _register_phase1(VR, reg, set(), {}, set(), {}, + vr_to_reaction={"v1": "R1", "v3": "R2"}, share_variants=True) + assert stats["unmapped"] == 1 + # the unmapped variant is registered per-variant, the mapped ones still share among themselves + assert reg[("B", "v1", "output")] != reg[("B", "v2", "output")] + + +def test_unknown_reaction_falls_back_to_per_variant(): + reg = phase1(VR, {}, share=True) # no vr -> reaction map + assert reg[("B", "v1", "output")] != reg[("B", "v2", "output")]