Regulator wiring: set-pool nodes (default on), cap pools (off), drugs.csv

schema/logic_network.linkml.yaml

@@ -58,6 +58,12 @@ enums:
           an assembly edge.
       regulator_variant:
         description: A set-variant produced by catalyst / negative-regulator decomposition.
+      set_pool:
+        description: >-
+          An EntitySet catalyst or regulator kept as ONE node (LNG_SET_POOL,
+          deltasignal specs/033): its members feed it by set_member edges and it
+          feeds each reaction it serves once. diagram_entity_id is the set;
+          member_leaves are its members.
       reaction:
         description: A virtual-reaction node (edges route input->reaction->output).
       other:

src/neo4j_connector.py

@@ -1,7 +1,7 @@
 import os
 
 from src.credential_redaction import scrub
-from typing import Any, Dict, List, Optional, Set, Union
+from typing import Any, Dict, List, Optional, Set, Tuple, Union
 
 import pandas as pd
 from py2neo import Graph  # type: ignore
@@ -861,6 +861,77 @@ def get_reference_entity_id(entity_id: str) -> Union[str, None]:
 
 
 
+_drug_structure_cache: Dict[str, Tuple[bool, List[Tuple[str, str]]]] = {}
+_drug_meta_cache: Dict[str, Tuple[str, str]] = {}
+
+
+def get_drug_entities(stable_ids) -> Dict[str, Dict[str, str]]:
+    """The entities among ``stable_ids`` that are DRUG-DERIVED, with their
+    schema class and display name (deltasignal specs/032).
+
+    Drug-derived is decided structurally, from the release itself:
+      - a Reactome ``Drug`` (ChemicalDrug, ProteinDrug, RNADrug) is;
+      - a complex is if ANY component is (a MAP2K bound to trametinib);
+      - a set is only if EVERY member is, so a ligand set that merely includes
+        a therapeutic antibody does not flag its physiological members.
+
+    A benchmark describes a cell without the drug, and a consumer that holds
+    these at baseline needs the list to travel with the network rather than
+    keep its own copy, which is how the cofactor lists once diverged.
+    """
+    ids = sorted({s for s in stable_ids if s})
+    missing = [s for s in ids if s not in _drug_structure_cache]
+    if missing:
+        query = """
+            UNWIND $ids AS s
+            MATCH (n:PhysicalEntity {stId: s})-[:hasComponent|hasMember|hasCandidate*0..10]->(x)
+            WITH DISTINCT x
+            OPTIONAL MATCH (x)-[r:hasComponent|hasMember|hasCandidate]->(y)
+            RETURN x.stId AS x, x:Drug AS is_drug, x.schemaClass AS c, x.displayName AS d,
+                   collect(CASE WHEN y IS NULL THEN NULL ELSE [type(r), y.stId] END) AS kids
+        """
+        try:
+            rows = get_graph().run(query, ids=missing).data()
+        except Exception:
+            logger.error("Error in get_drug_entities", **_traceback_kwargs())
+            raise
+        for r in rows:
+            if r.get("x"):
+                _drug_structure_cache[r["x"]] = (
+                    bool(r["is_drug"]), [tuple(k) for k in r["kids"] if k and k[1]])
+                _drug_meta_cache[r["x"]] = (r.get("c") or "", r.get("d") or "")
+        for s in missing:
+            _drug_structure_cache.setdefault(s, (False, []))
+
+    memo: Dict[str, bool] = {}
+
+    def derived(s: str, depth: int = 0) -> bool:
+        if s in memo:
+            return memo[s]
+        memo[s] = False                      # a malformed self-containing entity cannot recurse
+        is_drug, kids = _drug_structure_cache.get(s, (False, []))
+        comps = [y for t, y in kids if t == "hasComponent"]
+        members = [y for t, y in kids if t != "hasComponent"]
+        if is_drug:
+            out = True
+        elif depth >= 10:
+            out = False
+        elif comps:
+            out = any(derived(y, depth + 1) for y in comps)
+        elif members:
+            out = all(derived(y, depth + 1) for y in members)
+        else:
+            out = False
+        memo[s] = out
+        return out
+
+    # Class and name come back with the structure query: a second, unlabelled
+    # MATCH (n {stId: s}) scanned every node (review of PR #98: ~30 s a pathway).
+    return {s: {"schema_class": _drug_meta_cache.get(s, ("", ""))[0],
+                "name": _drug_meta_cache.get(s, ("", ""))[1]}
+            for s in ids if derived(s)}
+
+
 def get_reactome_release() -> Optional[int]:
     """The Reactome release number of the connected graph, or None.
 

src/pathway_generator.py

@@ -11,6 +11,7 @@
 from src.logic_network_generator import (
     create_pathway_logic_network,
     export_cofactors,
+    export_drugs,
     export_containment,
     export_entity_reaction_proxy_mapping,
     export_node_reaction_context,
@@ -35,6 +36,8 @@
     "LNG_HANDOFF_EDGES",
     "LNG_HANDOFF_HUB_MAX",
     "LNG_SET_MEMBERS_OR",
+    "LNG_SET_POOL",
+    "LNG_CAP_POOLS",
     "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
@@ -494,6 +497,19 @@ def generate_pathway_file(
             logger.error(f"Failed to write node provenance files: {e}", exc_info=True)
             # Don't raise - supplementary
 
+        # Drug-derived entities (deltasignal specs/032), in their OWN block: a
+        # failed drug query must not cost containment.csv, which the default
+        # solver needs (review of PR #98). A missing drugs.csv is reported by
+        # the consumer as "no drug table", not mistaken for "no drugs".
+        try:
+            export_drugs(
+                result.logic_network,
+                result.uuid_mapping,
+                str(pathway_output_dir / "drugs.csv"),
+            )
+        except Exception as e:
+            logger.error(f"Failed to write drugs.csv: {e}", exc_info=True)
+
         logger.info(f"Output directory: {pathway_output_dir}")
 
     except (ConnectionError, ValueError) as e:

src/reaction_generator.py

@@ -95,6 +95,10 @@ def _int_env(name: str, default: str) -> int:
 
 MAX_VARIANTS = _int_env("LNG_MAX_VARIANTS", "512")
 
+# Ids whose set combinations exceeded MAX_VARIANTS and were bundled into one
+# variant (deltasignal specs/036): LNG_CAP_POOLS pools only these.
+CAPPED_IDS: Set[str] = set()
+
 class _DecompositionStore:
     """Append-mostly buffer for decomposition rows with O(1) lookups.
 
@@ -288,6 +292,7 @@ def get_broken_apart_ids(
             merged: Set[str] = set()
             for member in new_broken_apart_members:
                 merged |= member
+            CAPPED_IDS.add(str(reactome_id))
             logger.warning(
                 f"variant cap hit for {reactome_id}: {cartesian_size} "
                 f"combinations > {MAX_VARIANTS}; bundling "

tests/test_cap_pools.py

@@ -0,0 +1,117 @@
+"""LNG_CAP_POOLS (deltasignal specs/036): in a reaction the variant cap BUNDLED
+(every alternative of every set merged into one variant), a set input becomes
+one pool: each alternative is an AND unit of its leaves, the pool ORs the
+alternatives, and the reaction reads the pool once.
+
+Shape: RAF "MAP2Ks and MAPKs bind to the activated RAF complex". Its input set S
+("RAF/MAPK scaffolds") has alternatives A (a single protein) and K (a complex,
+leaves K1, K2); over the cap all three leaves were separate required inputs.
+
+Review of #combined-037: the first version pooled ANY set resolving to several
+nodes, which also fires in UNCAPPED reactions where the one chosen alternative
+is a complex (miR-93 RISC in PIP3: its four subunits were averaged, so a miR-93
+knockout read 0.75x instead of 0)."""
+import pandas as pd
+import pytest
+
+import src.logic_network_generator as m
+from src import neo4j_connector
+
+S, C, A, K, K1, K2 = "R-HSA-S", "R-HSA-C", "R-HSA-A", "R-HSA-K", "R-HSA-K1", "R-HSA-K2"
+LEAVES = {S: {A, K1, K2}, A: {A}, K: {K1, K2}, C: {C}}
+
+
+@pytest.fixture
+def stub(monkeypatch):
+    monkeypatch.setattr(neo4j_connector, "get_reaction_input_output_ids",
+                        lambda rid, io: {S, C} if io == "input" else {"R-HSA-OUT"})
+    monkeypatch.setattr(neo4j_connector, "get_reaction_io_stoichiometry",
+                        lambda rid, io: {S: 2} if io == "input" else {})
+    monkeypatch.setattr(neo4j_connector, "get_labels",
+                        lambda e: ["CandidateSet"] if e == S else ["EntityWithAccessionedSequence"])
+    monkeypatch.setattr(neo4j_connector, "get_set_members", lambda e: [A, K] if e == S else [])
+    monkeypatch.setattr(m, "_matching_leaves", lambda e: LEAVES.get(e, {e}))
+    monkeypatch.setattr(m, "modifier_isoform_set_ids", lambda: set())
+    monkeypatch.setattr(m, "CAPPED_IDS", {"1"})
+
+
+def resolve(present, rid=1):
+    uid_index = {"in": ([], set(present) | {C}, {}), "out": ([], {"R-HSA-OUT"}, {})}
+    rmap = pd.DataFrame({"uid": ["vr1"], "input_hash": ["in"], "output_hash": ["out"], "reactome_id": [rid]})
+    return m._resolve_vr_entities(rmap, uid_index)
+
+
+def test_a_capped_bundle_becomes_a_pool_of_alternatives(stub, monkeypatch):
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    ve = resolve({A, K1, K2})
+    assert m._vr_input_pools == {"vr1": {S: {"alts": {A: {A}, K: {K1, K2}}, "stoich": 2}}}
+    assert set(ve["vr1"][0]) == {A, K1, K2, C}      # members stay VR inputs (Phase 2 joins)
+
+
+def test_an_uncapped_reaction_is_never_pooled(stub, monkeypatch):
+    # the miR-93 case: one chosen alternative that is a complex resolves to its
+    # subunits, which ARE co-required
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    resolve({K1, K2}, rid=2)
+    assert m._vr_input_pools == {}
+
+
+def test_the_cap_gate_itself(stub, monkeypatch):
+    # two alternatives present in a reaction the cap did NOT bundle: not pooled
+    # (only a capped reaction merged its alternatives)
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    resolve({A, K1, K2}, rid=2)
+    assert m._vr_input_pools == {}
+
+
+def test_one_alternative_present_is_not_a_pool(stub, monkeypatch):
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    resolve({K1, K2})
+    assert m._vr_input_pools == {}
+
+
+def test_a_leaf_another_input_also_maps_to_stays_direct(stub, monkeypatch):
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    monkeypatch.setattr(m, "_matching_leaves", lambda e: {A, K1, K2} if e == S else ({K1, K2, C} if e == C else LEAVES.get(e, {e})))
+    monkeypatch.setattr(neo4j_connector, "get_labels",
+                        lambda e: ["CandidateSet"] if e in (S,) else (["Complex"] if e == C else ["EntityWithAccessionedSequence"]))
+    monkeypatch.setattr(m, "_map_annotated_entity_to_nodes",
+                        lambda e, mem: {A, K1, K2} & mem if e == S else ({K1} if e == C else {e}))
+    resolve({A, K1, K2})
+    alts = m._vr_input_pools["vr1"][S]["alts"]
+    assert alts == {A: {A}, K: {K2}}                       # K1 is also C's node: not pooled
+
+
+def test_off_and_no_leak_between_pathways(stub, monkeypatch):
+    monkeypatch.delenv("LNG_CAP_POOLS", raising=False)
+    resolve({A, K1, K2})
+    assert m._vr_input_pools == {}
+    monkeypatch.setenv("LNG_CAP_POOLS", "1")
+    resolve({A, K1, K2})
+    assert m._vr_input_pools
+    resolve({A}, rid=2)
+    assert m._vr_input_pools == {}
+
+
+def test_emission_wires_alternatives_as_and_units_into_one_pool():
+    reg = {(x, "vr1", "input"): f"u-{x}" for x in (A, K1, K2, C)}
+    r2u, data = {}, []
+    pools = {S: {"alts": {A: {A}, K: {K1, K2}}, "stoich": 2}}
+    m._emit_vr_inputs("vr1", [A, K1, K2, C], {C: 1}, pools, reg, r2u, data, "R-HSA-RX")
+    into_vr = [(e["source_id"], e["edge_type"], e["stoichiometry"]) for e in data if e["target_id"] == "vr1"]
+    pool = next(u for u, s in r2u.items() if s == S)
+    alt = next(u for u, s in r2u.items() if s == K)
+    assert sorted(into_vr) == sorted([(pool, "input", 2), ("u-" + C, "input", 1)])  # one input per component
+    assert {(e["source_id"], e["and_or"], e["edge_type"]) for e in data if e["target_id"] == pool} == \
+        {("u-" + A, "or", "set_member"), (alt, "or", "set_member")}
+    assert {(e["source_id"], e["and_or"], e["edge_type"]) for e in data if e["target_id"] == alt} == \
+        {("u-" + K1, "and", "assembly"), ("u-" + K2, "and", "assembly")}
+    assert not any(e["target_id"] == "vr1" and e["source_id"] in ("u-" + A, "u-" + K1, "u-" + K2) for e in data)
+
+
+def test_without_pools_emission_is_the_plain_input_fan_in():
+    reg = {(x, "vr1", "input"): f"u-{x}" for x in (A, C)}
+    data = []
+    m._emit_vr_inputs("vr1", [A, C], {}, {}, reg, {}, data, "R-HSA-RX")
+    assert sorted(e["source_id"] for e in data) == ["u-" + A, "u-" + C]
+    assert {e["edge_type"] for e in data} == {"input"}

tests/test_provenance_exports.py

@@ -285,3 +285,94 @@ def test_export_cofactors_finds_an_entity_split_across_uuids(tmp_path, monkeypat
 
     df = pd.read_csv(out)
     assert int(df.in_network.sum()) == 1, "split entity missed when only one uuid is used"
+
+
+# --- drugs.csv (deltasignal specs/032) ---------------------------------------
+
+def test_export_drugs_lists_present_drug_entities_only(tmp_path, monkeypatch):
+    seen = {}
+    def fake(stids):
+        seen["asked"] = set(stids)
+        return {"R-HSA-D": {"schema_class": "ChemicalDrug", "name": "trametinib [cytosol]"}}
+    monkeypatch.setattr(neo4j_connector, "get_drug_entities", fake)
+    monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97)
+    edges = pd.DataFrame([{"source_id": "u-d", "target_id": "u-x"}])
+    out = tmp_path / "drugs.csv"
+    m.export_drugs(edges, {"R-HSA-D": "u-d", "R-HSA-X::variant::R-HSA-A": "u-x"}, str(out))
+    df = pd.read_csv(out)
+    assert list(df.columns) == ["stable_id", "schema_class", "name", "reactome_release"]
+    assert list(df.stable_id) == ["R-HSA-D"] and set(df.reactome_release) == {97}
+    # A variant is judged by its parent AND the members it chose.
+    assert seen["asked"] == {"R-HSA-D", "R-HSA-X", "R-HSA-A"}
+
+
+def test_export_drugs_lists_a_variant_by_its_exact_node_id(tmp_path, monkeypatch):
+    # Set S mixes a drug D and a protein P, so S is not drug-derived; the variant
+    # that CHOSE D is, and it must be listed under the exact id a consumer reads
+    # (stid_to_uuid_mapping.csv), not the parent's.
+    monkeypatch.setattr(neo4j_connector, "get_drug_entities",
+                        lambda stids: {"R-HSA-D": {"schema_class": "ChemicalDrug", "name": "d"}})
+    monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97)
+    vd, vp = "R-HSA-S::variant::R-HSA-D", "R-HSA-S::variant::R-HSA-P"
+    edges = pd.DataFrame([{"source_id": "u-1", "target_id": "u-2"}])
+    out = tmp_path / "drugs.csv"
+    m.export_drugs(edges, {vd: "u-1", vp: "u-2"}, str(out))
+    assert list(pd.read_csv(out).stable_id) == [vd]
+
+
+def test_export_drugs_writes_a_header_only_file_when_none(tmp_path, monkeypatch):
+    monkeypatch.setattr(neo4j_connector, "get_drug_entities", lambda s: {})
+    monkeypatch.setattr(neo4j_connector, "get_reactome_release", lambda: 97)
+    out = tmp_path / "drugs.csv"
+    m.export_drugs(pd.DataFrame([{"source_id": "u-x", "target_id": "u-y"}]), {"R-HSA-X": "u-x"}, str(out))
+    assert out.exists() and len(pd.read_csv(out)) == 0
+
+
+def test_drug_rule_complex_any_component_set_every_member(monkeypatch):
+    # D is a drug. C1 = complex(P, D) is drug-derived; S_all = set{D, D2} is;
+    # S_mixed = set{D, P} is NOT (its physiological member must not be held);
+    # C2 = complex(P, S_mixed) is NOT either; P is not.
+    struct = {
+        "D": (True, []), "D2": (True, []), "P": (False, []),
+        "C1": (False, [("hasComponent", "P"), ("hasComponent", "D")]),
+        # candidate-only, so dropping hasCandidate from the rule is caught
+        "S_all": (False, [("hasCandidate", "D"), ("hasCandidate", "D2")]),
+        "S_mixed": (False, [("hasMember", "D"), ("hasMember", "P")]),
+        "C2": (False, [("hasComponent", "P"), ("hasComponent", "S_mixed")]),
+    }
+    monkeypatch.setattr(neo4j_connector, "_drug_structure_cache", dict(struct))
+
+    monkeypatch.setattr(neo4j_connector, "_drug_meta_cache", {})
+
+    def boom():
+        raise AssertionError("fully cached: no query expected")
+    monkeypatch.setattr(neo4j_connector, "get_graph", boom)
+    got = neo4j_connector.get_drug_entities(struct)
+    assert set(got) == {"D", "D2", "C1", "S_all"}
+
+
+def test_get_drug_entities_queries_uncached_ids_and_keeps_misses(monkeypatch):
+    # The Cypher path: the structure query returns one row per reachable entity
+    # (with class and name), and an id it does not return (a reaction stid) is
+    # cached as not-a-drug instead of being re-queried.
+    monkeypatch.setattr(neo4j_connector, "_drug_structure_cache", {})
+    monkeypatch.setattr(neo4j_connector, "_drug_meta_cache", {})
+    calls = []
+
+    class G:
+        def run(self, q, **kw):
+            calls.append(sorted(kw["ids"]))
+
+            class R:
+                def data(_):
+                    return [{"x": "C", "is_drug": False, "c": "Complex", "d": "P:D",
+                             "kids": [["hasComponent", "P"], ["hasComponent", "D"]]},
+                            {"x": "P", "is_drug": False, "c": "EWAS", "d": "P", "kids": [None]},
+                            {"x": "D", "is_drug": True, "c": "ChemicalDrug", "d": "d", "kids": []}]
+            return R()
+    monkeypatch.setattr(neo4j_connector, "get_graph", lambda: G())
+    got = neo4j_connector.get_drug_entities({"C", "R-RXN"})
+    assert got == {"C": {"schema_class": "Complex", "name": "P:D"}}
+    assert calls == [["C", "R-RXN"]]
+    assert neo4j_connector.get_drug_entities({"C", "R-RXN"}) == got
+    assert len(calls) == 1                         # second call served from the cache

tests/test_regulators_and_catalysts.py

@@ -8,6 +8,7 @@
 5. Regulatory relationships are properly created
 """
 
+import os
 import pytest
 import pandas as pd
 from typing import Dict, List, Any
@@ -29,6 +30,7 @@ def _mock_decompose(entity_id, variant_decomposition=False, bundle_complex=False
     return [(entity_id, 1)]
 
 
+@patch.dict(os.environ, {"LNG_SET_POOL": "0"})   # these pin the member fan-out; pooling: tests/test_set_pool.py
 class TestRegulatorsAndCatalysts:
     """Test regulatory and catalytic relationships in logic networks."""