Faithful set-variant networks + LinkML output contract + diagram connectivity + determinism

.gitignore

@@ -47,3 +47,4 @@ uuid_mapping_*.csv
 reaction_connections_*.csv
 decomposed_uid_mapping_*.csv
 best_matches_*.csv
+graphify-out/

bin/backfill-proxy-mapping.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python3
+"""Backfill entity_reaction_proxy_mapping.csv for already-generated pathways.
+
+Regenerating a network from scratch is expensive; this reuses the existing
+logic_network.csv + stid_to_uuid_mapping.csv and only adds the new proxy file
+(which just needs Neo4j + those two artifacts). Going forward the generator
+emits the file itself; this is a one-time catch-up for the existing catalog.
+"""
+import os
+import sys
+from pathlib import Path
+
+import pandas as pd
+
+# Allow running as `bin/backfill-proxy-mapping.py` — put the repo root on the
+# path so `src` imports resolve regardless of cwd.
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
+
+from src.neo4j_connector import prefetch_entity_data  # noqa: E402
+from src.logic_network_generator import export_entity_reaction_proxy_mapping  # noqa: E402
+
+OUTPUT = Path(__file__).resolve().parent.parent / "output"
+
+
+def backfill(pathway_dir: Path) -> None:
+    name = pathway_dir.name
+    if not name.endswith(tuple(f"R-HSA-{n}" for n in [""])) and "R-HSA-" not in name:
+        return
+    pid = "R-HSA-" + name.rsplit("R-HSA-", 1)[-1]
+    net_f = pathway_dir / "logic_network.csv"
+    map_f = pathway_dir / "stid_to_uuid_mapping.csv"
+    rc_f = pathway_dir / "cache" / "reaction_connections.csv"
+    if not (net_f.exists() and map_f.exists()):
+        print(f"  skip {name}: missing artifacts")
+        return
+
+    network = pd.read_csv(net_f)
+    mapping = pd.read_csv(map_f)
+    reactome_id_to_uuid = dict(zip(mapping["uuid"].astype(str),
+                                   mapping["stable_id"].astype(str)))
+    # Approximate reaction_id_map from the mapping (reaction stIds carry UUIDs
+    # too); only reaction stIds get queried, so entity rows are harmless.
+    rid_map = mapping.rename(columns={"uuid": "uid", "stable_id": "reactome_id"})
+
+    # Warm caches with the pathway's reactions for fast terminal-component walks.
+    if rc_f.exists():
+        rc = pd.read_csv(rc_f, dtype=str)
+        rxn_ids = set(rc["preceding_reaction_id"].dropna()) | \
+                  set(rc["following_reaction_id"].dropna())
+        if rxn_ids:
+            prefetch_entity_data(list(rxn_ids))
+
+    out_f = pathway_dir / "entity_reaction_proxy_mapping.csv"
+    export_entity_reaction_proxy_mapping(
+        network, rid_map, reactome_id_to_uuid, pid, str(out_f))
+    df = pd.read_csv(out_f)
+    print(f"  {name}: {len(df)} rows, "
+          f"{df['entity_stable_id'].nunique() if len(df) else 0} entities")
+
+
+if __name__ == "__main__":
+    targets = sys.argv[1:]
+    if targets:
+        dirs = [OUTPUT / t for t in targets]
+    else:
+        dirs = [d for d in sorted(OUTPUT.iterdir()) if d.is_dir()]
+    for d in dirs:
+        try:
+            backfill(d)
+        except Exception as e:
+            print(f"  ERROR {d.name}: {e}")

bin/create-pathways.py

@@ -3,6 +3,17 @@
 
 import os
 import sys
+
+# Determinism: the generator's output depends on set/dict iteration order, which
+# Python randomizes per-process via hash seeding — so regenerating a pathway
+# yields structurally different networks run-to-run (~5.8% of edges for TP53).
+# Pin the hash seed so generation is reproducible. PYTHONHASHSEED must be set
+# before the interpreter starts, so re-exec once if it isn't already fixed.
+# Override with LNG_ALLOW_NONDETERMINISM=1. See reactome/logic-network-generator#42.
+if os.environ.get("PYTHONHASHSEED") != "0" and os.environ.get("LNG_ALLOW_NONDETERMINISM") != "1":
+    os.environ["PYTHONHASHSEED"] = "0"
+    os.execv(sys.executable, [sys.executable] + sys.argv)
+
 from typing import List, Tuple
 
 import pandas as pd

bin/validate-against-mpbiopath.py

@@ -44,6 +44,27 @@
 DOWN, NORMAL, UP = 0, 1, 2
 MAX_ITERATIONS = 50
 
+# AND-combination mode for the propagator. "min" is Boolean MIN (rate-limiting
+# reagent wins; can't carry upregulations through gates whose other inputs
+# are unperturbed). "signed" sums each input's displacement from NORMAL and
+# uses the sign: net up → UP, net down → DOWN, balanced or unperturbed →
+# NORMAL. Switch via --propagator on the CLI.
+PROPAGATOR_MODE = "min"
+
+
+def signed_and(contribs: list[int]) -> int:
+    """Sum-of-displacements AND. (state - NORMAL) per input; sign decides output."""
+    displacement = sum(c - NORMAL for c in contribs)
+    if displacement > 0:
+        return UP
+    if displacement < 0:
+        return DOWN
+    return NORMAL
+
+
+def combine_and(contribs: list[int]) -> int:
+    return signed_and(contribs) if PROPAGATOR_MODE == "signed" else min(contribs)
+
 
 def invert(state: int) -> int:
     return {DOWN: UP, NORMAL: NORMAL, UP: DOWN}[state]
@@ -71,6 +92,26 @@ def build_stid_to_uuids(pathway_dir: Path) -> dict[str, list[str]]:
     return out
 
 
+def load_entity_reaction_proxies(pathway_dir: Path) -> dict[str, list[str]]:
+    """entity stable_id → list of proxy reaction UUIDs.
+
+    Curated species (often Complexes containing an EntitySet) get expanded into
+    virtual variants during generation, so the parent's stId isn't in
+    stid_to_uuid_mapping.csv. The generator's entity_reaction_proxy_mapping.csv
+    points each such species at the UUIDs of the reaction that produces it, so
+    we can read reaction flux as a proxy for the species' state. Absent file →
+    no proxies (pathway was generated by an older generator).
+    """
+    out: dict[str, list[str]] = defaultdict(list)
+    f = pathway_dir / "entity_reaction_proxy_mapping.csv"
+    if not f.exists():
+        return out
+    proxies = pd.read_csv(f)
+    for _, row in proxies.iterrows():
+        out[str(row["entity_stable_id"])].append(str(row["proxy_uuid"]))
+    return out
+
+
 def gene_name_to_stids(graph, gene_names: list[str]) -> dict[str, list[str]]:
     """Map gene names to all PhysicalEntity stable IDs whose reference entity has that gene name."""
     rows = graph.run(
@@ -129,9 +170,9 @@ def propagate(
                     or_contribs.append(src_state)
 
             if and_contribs and or_contribs:
-                new_state[uuid_] = max(min(and_contribs), max(or_contribs))
+                new_state[uuid_] = max(combine_and(and_contribs), max(or_contribs))
             elif and_contribs:
-                new_state[uuid_] = min(and_contribs)
+                new_state[uuid_] = combine_and(and_contribs)
             elif or_contribs:
                 new_state[uuid_] = max(or_contribs)
         if new_state == state:
@@ -254,11 +295,19 @@ def validate_one_pathway(
         return {"status": "no_curator_file", "name": pathway_name}
     network = load_network(pathway_dir)
     stid_to_uuids = build_stid_to_uuids(pathway_dir)
+    entity_reaction_proxies = load_entity_reaction_proxies(pathway_dir)
 
-    # Resolve key outputs (numeric dbId → list of UUIDs)
+    # Resolve key outputs (numeric dbId → list of UUIDs). Fall back to the
+    # generator's proxy mapping when the curated species was decomposed into
+    # virtual variants and isn't directly addressable — read the producing
+    # reaction's flux as a proxy for the species' state.
     key_output_uuids: dict[str, list[str]] = {}
     for ko in curator["key_output"].astype(str):
-        key_output_uuids[ko] = stid_to_uuids.get(f"R-HSA-{ko}", [])
+        sid = f"R-HSA-{ko}"
+        uuids = stid_to_uuids.get(sid, [])
+        if not uuids:
+            uuids = entity_reaction_proxies.get(sid, [])
+        key_output_uuids[ko] = uuids
 
     # Resolve perturbation genes
     perturbations = parse_perturbation_columns(curator)
@@ -366,7 +415,14 @@ def main():
              "Experimental is only available for 10 pathways and contains -999 (not measured) "
              "cells which are skipped.",
     )
+    ap.add_argument(
+        "--propagator", choices=["min", "signed"], default="min",
+        help="AND-combination rule. 'min' = Boolean MIN (default); 'signed' = sum of "
+             "displacements from NORMAL, sign decides direction (handles upregulation symmetrically).",
+    )
     args = ap.parse_args()
+    global PROPAGATOR_MODE
+    PROPAGATOR_MODE = args.propagator
 
     output_root = Path(args.output_dir)
     pathways = pd.read_csv("/tmp/mpbio_pathways.tsv", sep="\t")

bin/validate-logic-network.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+"""Validate a generated logic-network output dir against the LinkML schema.
+
+The CSVs are the flat serialization; multivalued slots (member_leaves,
+source_sets, chosen_members) are pipe-delimited cells. This loads them back
+into the native LinkML data model (real arrays) and validates the assembled
+``LogicNetwork`` instance against schema/logic_network.linkml.yaml.
+
+Usage:
+    poetry run python bin/validate-logic-network.py <pathway_output_dir>
+"""
+import sys
+from pathlib import Path
+
+import pandas as pd
+from linkml.validator import validate
+
+SCHEMA = Path(__file__).resolve().parent.parent / "schema" / "logic_network.linkml.yaml"
+MULTIVALUED = ("member_leaves", "source_sets", "chosen_members")
+
+
+def _split(v):
+    if v is None or (isinstance(v, float) and pd.isna(v)) or v == "":
+        return []
+    return [p for p in str(v).split("|") if p]
+
+
+def _clean(v):
+    return None if (v is None or (isinstance(v, float) and pd.isna(v)) or v == "") else v
+
+
+def load_instance(d: Path) -> dict:
+    nodes = []
+    for r in pd.read_csv(d / "nodes.csv", dtype=str, keep_default_na=False).to_dict("records"):
+        node = {k: _clean(v) for k, v in r.items() if k not in MULTIVALUED}
+        for k in MULTIVALUED:
+            node[k] = _split(r.get(k))
+        nodes.append({k: v for k, v in node.items() if v is not None or k in MULTIVALUED})
+
+    contexts = pd.read_csv(d / "node_reaction_context.csv", dtype=str,
+                           keep_default_na=False).to_dict("records")
+
+    edges = []
+    for r in pd.read_csv(d / "logic_network.csv", keep_default_na=False).to_dict("records"):
+        e = {k: _clean(v) for k, v in r.items()}
+        if e.get("stoichiometry") not in (None, ""):
+            e["stoichiometry"] = int(float(e["stoichiometry"]))
+        edges.append({k: v for k, v in e.items() if v is not None})
+
+    return {"pathway_id": d.name, "nodes": nodes,
+            "node_reaction_contexts": contexts, "edges": edges}
+
+
+def main() -> int:
+    if len(sys.argv) != 2:
+        print(__doc__)
+        return 2
+    d = Path(sys.argv[1])
+    report = validate(load_instance(d), str(SCHEMA), "LogicNetwork")
+    if not report.results:
+        print(f"VALID: {d} conforms to {SCHEMA.name}")
+        return 0
+    print(f"INVALID: {len(report.results)} problem(s) in {d}")
+    for res in report.results[:25]:
+        print(f"  [{res.severity}] {res.message}")
+    return 1
+
+
+if __name__ == "__main__":
+    sys.exit(main())