Derive the cofactor set, and carry the realisation links the diagram draws

src/diagram_connectivity.py

@@ -129,6 +129,81 @@ def diagram_shared_product_pairs(pathway_id: str) -> Set[Tuple[str, str]]:
     return pairs
 
 
+def diagram_set_member_pairs(pathway_id: str) -> Set[Tuple[str, str]]:
+    """(member_stId, set_stId) pairs the DIAGRAM draws and Neo4j does not have.
+
+    A diagram layout carries a `links` array alongside its `edges`, and nothing
+    in this pipeline read it. Most of those links restate set membership that
+    Neo4j already holds as hasMember, so they are redundant — but not all.
+
+    Worked example, verified against Release97. `SMAD7:SMURF/NEDD4L`
+    (R-HSA-2169026) and `SMAD7:SMURF2` (R-HSA-2167883) are both curated as
+    COMPLEXES with no containment relation between them; the generic one holds
+    a `SMURF/NEDD4L` DefinedSet where the specific holds SMURF2. Semantically
+    the specific realises the generic, the diagram says so with an
+    EntitySetAndMemberLink, and the generated network had them as two
+    unconnected nodes. Same shape for the EPH-ephrin oligomer complexes.
+
+    Only `EntitySetAndMemberLink` is returned. `EntitySetAndEntitySetLink` is a
+    set-to-set overlap rather than a realisation, and `Interaction` / `FlowLine`
+    point at entities that usually have no node at all — consuming those means
+    ADDING nodes, which is a different and much larger change (issue #41).
+
+    The link's `inputs` are the member glyph and its `outputs` the set glyph,
+    so the pair is returned member-first: the member realises the set.
+    """
+    ddir = _diagram_dir()
+    diagram_stid = _covering_diagram_stid(pathway_id)
+    if not diagram_stid:
+        return set()
+    try:
+        layout = json.loads((ddir / f"{diagram_stid}.json").read_text())
+        graph = json.loads((ddir / f"{diagram_stid}.graph.json").read_text())
+    except Exception:
+        logger.warning(f"{pathway_id}: could not read diagram links")
+        return set()
+
+    glyph_to_dbid = {n["id"]: n.get("reactomeId")
+                     for n in layout.get("nodes", []) if n.get("reactomeId")}
+    dbid_to_stid = {n["dbId"]: n["stId"]
+                    for n in graph.get("nodes", []) if n.get("stId")}
+
+    # A pathway with no diagram of its own borrows an ancestor's, which also
+    # carries its SIBLINGS' links. `diagram_shared_product_pairs` filters those
+    # out by reaction; do the equivalent here by entity, so a sub-pathway does
+    # not import a realisation drawn in another sub-pathway's context.
+    own_entities: Set[str] = set()
+    if diagram_stid != pathway_id:
+        try:
+            from src.neo4j_connector import get_pathway_participating_entities
+            own_entities = set(get_pathway_participating_entities(pathway_id))
+        except Exception:
+            logger.warning(
+                f"{pathway_id}: borrowing diagram {diagram_stid} but could not "
+                "scope its links to this pathway; skipping them rather than "
+                "importing a sibling's."
+            )
+            return set()
+
+    pairs: Set[Tuple[str, str]] = set()
+    for link in layout.get("links", []) or []:
+        if link.get("renderableClass") != "EntitySetAndMemberLink":
+            continue
+        members = {dbid_to_stid.get(glyph_to_dbid.get(x.get("id")))
+                   for x in (link.get("inputs") or [])}
+        sets = {dbid_to_stid.get(glyph_to_dbid.get(x.get("id")))
+                for x in (link.get("outputs") or [])}
+        for member in filter(None, members):
+            for parent in filter(None, sets):
+                if member == parent:
+                    continue
+                if own_entities and (member not in own_entities
+                                     and parent not in own_entities):
+                    continue
+                pairs.add((member, parent))
+    return pairs
+
+
 def augment_reaction_connections(pathway_id: str,
                                  reaction_connections: pd.DataFrame) -> pd.DataFrame:
     """Union diagram-drawn product->substrate pairs into reaction_connections.

src/logic_network_generator.py

@@ -936,22 +936,75 @@ def _expand_complex_variants(complex_id: str) -> List[tuple]:
     return variants if variants else [(complex_id, 1)]
 
 
-_COFACTOR_STIDS: frozenset = frozenset({
-    "R-ALL-113592",  # ATP
-    "R-ALL-29358",   # ATP variant
-    "R-ALL-113582",  # ADP
-    "R-ALL-29370",   # ADP variant
-    "R-ALL-29360",   # ADP variant
-    "R-ALL-29356",   # H2O
-    "R-ALL-29372",   # Pi
-    "R-ALL-29390",   # Pi variant
-    "R-ALL-29438",   # PPi
-    "R-ALL-217093",  # NADP+
-    "R-ALL-110114",  # NADPH
-    "R-ALL-29986",   # NAD+
-    "R-ALL-73473",   # NADH
+# Offline seed for the cofactor set. Every id here was verified against
+# Release97; the list this replaced had SIX of thirteen wrong, which is why
+# the live set is now derived rather than typed:
+#
+#   R-ALL-29438 was commented "PPi" and is GTP     <- still a cofactor
+#   R-ALL-29360 was commented "ADP variant" and is NAD+  <- still a cofactor
+#   R-ALL-29390 was commented "Pi variant" and is PXLP (pyridoxal 5'-phosphate)
+#   R-ALL-217093 / R-ALL-110114 / R-ALL-29986 do not exist in Release97 at all
+#
+# Note the distinction, which an earlier version of this seed got wrong: a
+# WRONG COMMENT is not a wrong ENTRY. GTP and NAD+ are both in
+# `_COFACTOR_CHEBI`, so dropping them because their labels were wrong would
+# make the offline path exclude FEWER real cofactors than the list it replaced.
+# Only PXLP was a genuine false positive; the three absent ids were inert.
+#
+# Used only when Neo4j is unreachable. `get_cofactor_species()` derives the
+# real set by ChEBI identity against the connected release — 253 species at
+# Release97, every compartment variant included.
+_COFACTOR_STIDS_SEED: frozenset = frozenset({
+    "R-ALL-113592",  # ATP [cytosol]
+    "R-ALL-29358",   # ATP [nucleoplasm]
+    "R-ALL-113582",  # ADP [nucleoplasm]
+    "R-ALL-29370",   # ADP [cytosol]
+    "R-ALL-29356",   # H2O [cytosol]
+    "R-ALL-29372",   # Pi [cytosol]
+    "R-ALL-73473",   # NADH [cytosol]
+    "R-ALL-29438",   # GTP [cytosol]   (was mis-commented "PPi")
+    "R-ALL-29360",   # NAD+ [cytosol]  (was mis-commented "ADP variant")
 })
 
+_cofactor_stids_cache: Optional[frozenset] = None
+
+
+def _cofactor_stids() -> frozenset:
+    """Stable ids treated as metabolic cofactors, derived from the release.
+
+    A cofactor is too shared for a depletion edge, a handoff leaf match or a
+    diagram bridge to mean anything — every reaction in a pathway touches ATP.
+    Deriving the set keeps it honest across releases; the hardcoded seed is
+    only a fallback for when Neo4j is unavailable, and it is deliberately
+    minimal rather than a stale snapshot.
+    """
+    global _cofactor_stids_cache
+    if _cofactor_stids_cache is not None:
+        return _cofactor_stids_cache
+    try:
+        from src.neo4j_connector import get_cofactor_species
+        derived = frozenset(e["stable_id"] for e in get_cofactor_species())
+    except Exception:
+        # Deliberately NOT cached. Memoising the failure would let one
+        # transient connection reset on pathway 1 silently build pathways
+        # 2..N with 9 cofactors instead of 253 — while `export_cofactors`
+        # queries Neo4j separately, succeeds, and ships a cofactors.csv
+        # listing all 253 beside a network built without them. Retrying per
+        # call is cheap next to a whole catalog built wrong.
+        logger.warning(
+            "Could not derive the cofactor set from Neo4j; using the offline "
+            "seed for THIS call only. Depletion edges, handoff leaves and "
+            "diagram bridges will treat fewer species as cofactors than they "
+            "should. Networks built now are not comparable with ones built "
+            "while the database was reachable."
+        )
+        return _COFACTOR_STIDS_SEED
+    if not derived:
+        logger.warning("Neo4j returned no cofactor species; using the offline seed.")
+        return _COFACTOR_STIDS_SEED
+    _cofactor_stids_cache = derived
+    return _cofactor_stids_cache
+
 # Ubiquitin entity stIds (human + cross-species variants). A reaction that
 # takes one of these as INPUT is a ubiquitination reaction (Ub is consumed
 # and attached to a target protein). Reactions whose OUTPUT is Ub are
@@ -971,7 +1024,14 @@ def _expand_complex_variants(complex_id: str) -> List[tuple]:
 # caller excluded cofactor hubs, but nothing did. Curators sometimes draw a
 # single shared Ub glyph, which is how two glyphs became 134 of S Phase's 158
 # bridges. See issue #61.
-_BRIDGE_EXCLUDED_STIDS: frozenset = _COFACTOR_STIDS | _UBIQUITIN_STIDS
+def _bridge_excluded_stids() -> frozenset:
+    """Species a diagram bridge must never be drawn across.
+
+    A function rather than a module constant because the cofactor half is
+    derived from the connected release; binding it at import time would freeze
+    whatever the seed happened to be.
+    """
+    return _cofactor_stids() | _UBIQUITIN_STIDS
 
 
 def _emit_substrate_depletion_edges(
@@ -1151,7 +1211,7 @@ def _emit_substrate_depletion_edges(
                 inp_stid = reactome_id_to_uuid.get(inp_uuid, "")
                 if inp_stid == cat_stid and inp_stid:
                     continue  # same biological entity at different positions
-                if inp_stid in _COFACTOR_STIDS:
+                if inp_stid in _cofactor_stids():
                     continue
                 key = (cat_uuid, inp_uuid)
                 if key in seen_edges:
@@ -1183,7 +1243,7 @@ def _emit_substrate_depletion_edges(
             cat_stid = reactome_id_to_uuid.get(cat_uuid, "")
             for subst_stid in subst_stids:
                 if subst_stid == cat_stid: continue
-                if subst_stid in _COFACTOR_STIDS: continue
+                if subst_stid in _cofactor_stids(): continue
                 target_uuids = stid_to_uuids_in_net.get(subst_stid, [])
                 for tgt_uuid in target_uuids:
                     if tgt_uuid == cat_uuid: continue
@@ -1227,7 +1287,7 @@ def _node_leaves(node_id: str) -> frozenset:
             s = set(get_terminal_components(node_id)) if "Complex" in get_labels(node_id) else {node_id}
         except Exception:
             s = {node_id}
-    leaves = frozenset(s - _COFACTOR_STIDS - _UBIQUITIN_STIDS)
+    leaves = frozenset(s - _cofactor_stids() - _UBIQUITIN_STIDS)
     _handoff_leaf_cache[node_id] = leaves
     return leaves
 
@@ -1329,6 +1389,85 @@ def _emit_precedingevent_handoff_edges(
     )
 
 
+def _emit_diagram_set_member_edges(
+    pathway_logic_network_data: List[Dict[str, Any]],
+    reactome_id_to_uuid: Dict[str, str],
+    set_member_pairs: Optional[Set[Tuple[str, str]]],
+) -> int:
+    """Connect a specific complex to the generic one the DIAGRAM says it realises.
+
+    Curators sometimes draw a realisation relationship between two entities
+    that Reactome stores with no containment between them — a generic complex
+    holding a DefinedSet component, and the specific complexes that instantiate
+    it. `SMAD7:SMURF2` -> `SMAD7:SMURF/NEDD4L` is the worked example. Without
+    this the two are unconnected nodes and a perturbation of the specific form
+    never reaches the generic one.
+
+    Direction is member -> set: more of the specific means more of the generic
+    pool. `or` because any member realises it, so the edge never imposes
+    AND-completeness on the target.
+
+    OFF BY DEFAULT, and the reason is worth stating. The target already carries
+    `and` assembly edges to every constituent protein, and DeltaSignal combines
+    an AND cluster with an OR cluster as `max(and, or)`, so a DECREASE arriving
+    on this edge is discarded while an increase is not: measured on TGF-beta,
+    driving the members down leaves the generic at 1.0, driving them up reaches
+    80.0. An edge that propagates one direction only is worse than no edge,
+    because it looks like it works. Enable together with DS_OR_COMBINE=gate,
+    never alone.
+
+    Both endpoints must already be nodes. A link whose other end is absent
+    means the entity participates in no curated reaction here, and inventing a
+    node for it is a much larger change (issue #41).
+    """
+    if not set_member_pairs:
+        return 0
+    by_stid: Dict[str, List[str]] = {}
+    for node_uuid, stid in reactome_id_to_uuid.items():
+        by_stid.setdefault(str(stid), []).append(str(node_uuid))
+    seen = {(e["source_id"], e["target_id"]) for e in pathway_logic_network_data}
+
+    # ONE edge per stable-id pair, not the cartesian product of occurrences.
+    # Positional decomposition gives an entity many uuids, so all-pairs turns
+    # two curated relationships into 48 edges in EPH-Ephrin alone — the same
+    # blow-up that made the all-pairs silo bridge unusable. Connect the
+    # best-connected occurrence on each side, which is where follow-on signal
+    # has somewhere to go.
+    degree: Dict[str, int] = {}
+    for e in pathway_logic_network_data:
+        degree[e["source_id"]] = degree.get(e["source_id"], 0) + 1
+        degree[e["target_id"]] = degree.get(e["target_id"], 0) + 1
+
+    # Ties are common (several occurrences each with degree 1) and a uuid4 is
+    # regenerated every run, so breaking ties on the uuid STRING would attach
+    # the edge to a different occurrence run to run — defeating the
+    # reproducibility the pinned PYTHONHASHSEED exists to give. Break on
+    # insertion order instead, which follows the deterministic build.
+    ordinal = {u: i for i, u in enumerate(reactome_id_to_uuid)}
+
+    emitted = 0
+    for member_stid, set_stid in sorted(set_member_pairs):
+        srcs = by_stid.get(member_stid, [])
+        tgts = by_stid.get(set_stid, [])
+        if not srcs or not tgts:
+            continue
+        src = max(srcs, key=lambda u: (degree.get(u, 0), -ordinal.get(u, 0)))
+        tgt = max(tgts, key=lambda u: (degree.get(u, 0), -ordinal.get(u, 0)))
+        if src == tgt or (src, tgt) in seen:
+            continue
+        seen.add((src, tgt))
+        pathway_logic_network_data.append({
+            "source_id": src,
+            "target_id": tgt,
+            "pos_neg": "pos",
+            "and_or": "or",
+            "edge_type": "diagram_set_member",
+            "stoichiometry": 1,
+        })
+        emitted += 1
+    return emitted
+
+
 def _emit_boundary_decomposition_edges(
     pathway_logic_network_data: List[Dict[str, Any]],
     reactome_id_to_uuid: Dict[str, str],
@@ -1638,6 +1777,7 @@ def create_pathway_logic_network(
     reaction_connections: pd.DataFrame,
     best_matches: Any,
     diagram_bridge_pairs: Optional[Set[Tuple[str, str]]] = None,
+    diagram_set_member_pairs: Optional[Set[Tuple[str, str]]] = None,
 ) -> PathwayResult:
     """Create a pathway logic network from decomposed UID mappings and reaction connections.
 
@@ -1870,12 +2010,17 @@ def create_pathway_logic_network(
     # only (cofactor hubs already excluded by the caller). See #39.
     if diagram_bridge_pairs:
         n_bridges = 0
+        # Hoisted: this was a module constant before the set became
+        # derived, and rebuilding a ~256-element union per
+        # (pair x producer_vr x consumer_vr) triple is 10^5-10^6
+        # allocations on a variant-expanded pathway.
+        bridge_excluded = _bridge_excluded_stids()
         for a_rid, b_rid in diagram_bridge_pairs:
             for p_vr in reactome_to_vr.get(a_rid, []):
                 p_outputs = set(vr_entities.get(p_vr, ([], [], {}, {}))[1])
                 for f_vr in reactome_to_vr.get(b_rid, []):
                     f_inputs = set(vr_entities.get(f_vr, ([], [], {}, {}))[0])
-                    for eid in (p_outputs & f_inputs) - _BRIDGE_EXCLUDED_STIDS:
+                    for eid in (p_outputs & f_inputs) - bridge_excluded:
                         src = entity_uuid_registry.get((eid, p_vr, "output"))
                         tgt = entity_uuid_registry.get((eid, f_vr, "input"))
                         # Skip if missing or already the same node (already
@@ -1951,6 +2096,17 @@ def create_pathway_logic_network(
         reactome_id_to_uuid=reactome_id_to_uuid,
     )
 
+    # Realisation links the DIAGRAM draws between a specific complex and the
+    # generic one it instantiates, where Reactome stores no containment. See
+    # _emit_diagram_set_member_edges.
+    n_set_member = _emit_diagram_set_member_edges(
+        pathway_logic_network_data=pathway_logic_network_data,
+        reactome_id_to_uuid=reactome_id_to_uuid,
+        set_member_pairs=diagram_set_member_pairs,
+    )
+    if n_set_member:
+        logger.info(f"Diagram set-member links: +{n_set_member} edges")
+
     # Restore curator-intended connectivity that complex-bundling drops: two
     # precedingEvent-linked reactions that hand off a shared COMPONENT (bound in
     # a complex on one side, free/other-complex on the other) share no whole

src/pathway_generator.py

@@ -35,6 +35,7 @@
     "LNG_HANDOFF_EDGES",
     "LNG_HANDOFF_HUB_MAX",
     "LNG_SET_MEMBERS_OR",
+    "LNG_DIAGRAM_SET_MEMBER",
     # Determinism controls: node ids are uuid4 and several selections iterate
     # sets, so hash seeding changes emitted content (~5.8% of TP53 edges per
     # bin/create-pathways.py). A cache built unseeded is not comparable to a
@@ -327,6 +328,7 @@ def generate_pathway_file(
         # LNG_DIAGRAM_CONNECTIVITY=0 to disable diagram connectivity entirely.
         from src.diagram_connectivity import (
             augment_reaction_connections,
+            diagram_set_member_pairs,
             diagram_shared_product_pairs,
         )
         diagram_bridge_pairs = None
@@ -338,9 +340,35 @@ def generate_pathway_file(
 
         # Generate logic network
         logger.info("Creating pathway logic network...")
+        # Realisation links the diagram draws between a specific complex and
+        # the generic one it instantiates. Has its own flag AND honours
+        # LNG_DIAGRAM_CONNECTIVITY=0, so the documented "disable diagram
+        # connectivity entirely" kill switch above stays true.
+        # DEFAULT OFF. The edge is correct in principle and inert in practice:
+        # the target complex already carries `and` assembly edges to every
+        # constituent protein, and DeltaSignal combines an AND cluster with an
+        # OR cluster as max(and, or), so a DECREASE through this edge is
+        # discarded — measured on TGF-beta, member DOWN leaves the generic at
+        # 1.0 while member UP reaches 80.0. One-directional propagation is
+        # worse than none, because it looks like it works.
+        #
+        # It becomes load-bearing under DS_OR_COMBINE=gate, which is itself
+        # default-off and measured inert precisely because nothing in the
+        # catalog produced mixed and/or clusters. This produces them. The two
+        # therefore have to be enabled and measured TOGETHER, and neither alone.
+        set_member_pairs = None
+        if (os.environ.get("LNG_DIAGRAM_SET_MEMBER", "0") == "1"
+                and os.environ.get("LNG_DIAGRAM_CONNECTIVITY", "1") != "0"):
+            try:
+                set_member_pairs = diagram_set_member_pairs(pathway_id)
+            except Exception:
+                logger.warning("Could not read diagram set-member links",
+                               exc_info=True)
+
         result = create_pathway_logic_network(
             decomposed_uid_mapping, connectivity, best_matches,
             diagram_bridge_pairs=diagram_bridge_pairs,
+            diagram_set_member_pairs=set_member_pairs,
         )
 
         # Save logic network (main output file users need)