specs/047: the experimental gap to MP-BioPath, case by case; adapter fix

bench/analysis/experimental_gap.py

@@ -0,0 +1,105 @@
+#!/usr/bin/env python3
+"""Where DeltaSignal loses to MP-BioPath on the experimental axis, case by case.
+
+MP-BioPath's per-case calls are published in
+mp-biopath-pathways/analysis_results/MP_BioPathReactomePathwayAccuracy-Tests-updated.tsv
+(`mp_biopath_state` against `experimental_value`, 849 cases, 643 correct).
+This joins them with an arm's experimental_cases.tsv and splits the cases into
+both right / only MP-BioPath / only DeltaSignal / neither, by pathway,
+perturbation, our failure category and (truth, our call).
+
+A case that does not join, or whose experimental value disagrees with our
+expected label, is an error, not a skip: a silent mis-join would make the gap
+look like whatever the join happened to keep.
+
+Usage: experimental_gap.py ARM_RESULTS_DIR [--mpb TSV] [--out cases.tsv]
+"""
+from __future__ import annotations
+
+import argparse
+import collections
+import csv
+import sys
+from pathlib import Path
+
+DEFAULT_MPB = (Path.home() / "gitroot" / "mp-biopath-pathways" / "analysis_results"
+               / "MP_BioPathReactomePathwayAccuracy-Tests-updated.tsv")
+LABEL = {"0": "DOWN", "1": "NORMAL", "2": "UP"}
+
+
+def load_mpb(path: Path) -> dict[tuple, tuple[str, str]]:
+    """(pathway, gene, direction, key_output) -> (experimental, mp-biopath state)."""
+    out: dict[tuple, tuple[str, str]] = {}
+    with path.open(newline="") as handle:
+        for r in csv.DictReader(handle, delimiter="\t"):
+            if r["experimental_value"] in ("NA", ""):
+                continue
+            gene, direction = r["scenario"].rsplit("_", 1)
+            key = (r["pathway_name"], gene, direction, r["keyoutput_id"])
+            if key in out:
+                raise SystemExit(f"duplicate MP-BioPath row {key}")
+            out[key] = (r["experimental_value"], r["mp_biopath_state"])
+    return out
+
+
+def split(mpb: dict, ours_rows: list[dict]) -> list[dict]:
+    rows = []
+    for o in ours_rows:
+        key = (o["pathway"], o["gene"], o["direction"], o["key_output"])
+        if key not in mpb:
+            raise SystemExit(f"case not in MP-BioPath's table: {key}")
+        truth, them = mpb[key]
+        if truth != o["expected"]:
+            raise SystemExit(f"experimental value {truth} != our expected {o['expected']} for {key}")
+        us = o["predicted"] if o["valid"] == "1" else "invalid"
+        rows.append({"pathway": key[0], "gene": key[1], "direction": key[2], "key_output": key[3],
+                     "truth": truth, "mpbiopath": them, "deltasignal": us,
+                     "quadrant": ("MPB" if them == truth else "-") + "/" + ("DS" if us == truth else "-"),
+                     "category": o.get("category", "")})
+    missing = set(mpb) - {(r["pathway"], r["gene"], r["direction"], r["key_output"]) for r in rows}
+    if missing:
+        raise SystemExit(f"{len(missing)} MP-BioPath cases have no row in the arm, e.g. {sorted(missing)[0]}")
+    return rows
+
+
+def report(rows: list[dict]) -> str:
+    out = []
+    q = collections.Counter(r["quadrant"] for r in rows)
+    mpb_right = q["MPB/DS"] + q["MPB/-"]; ds_right = q["MPB/DS"] + q["-/DS"]
+    out.append(f"cases {len(rows)}: MP-BioPath {mpb_right}, DeltaSignal {ds_right}, gap {mpb_right - ds_right:+d}")
+    out.append(f"  both right {q['MPB/DS']}, only MP-BioPath {q['MPB/-']}, "
+               f"only DeltaSignal {q['-/DS']}, neither {q['-/-']}; either {len(rows) - q['-/-']}")
+    pw = collections.defaultdict(collections.Counter)
+    for r in rows:
+        pw[r["pathway"]][r["quadrant"]] += 1
+    out.append("per pathway: only-MPB | only-DS | net")
+    for p, c in sorted(pw.items(), key=lambda x: -(x[1]["MPB/-"] - x[1]["-/DS"])):
+        out.append(f"  {p[:40]:40s} {c['MPB/-']:4d} | {c['-/DS']:4d} | {c['MPB/-'] - c['-/DS']:+4d}  (n={sum(c.values())})")
+    lost = [r for r in rows if r["quadrant"] == "MPB/-"]
+    out.append("only-MPB by our category: " + str(collections.Counter(r["category"] for r in lost).most_common()))
+    out.append("only-MPB by (truth, our call): " + str(collections.Counter(
+        (LABEL.get(r["truth"]), LABEL.get(r["deltasignal"], r["deltasignal"])) for r in lost).most_common()))
+    pert = collections.Counter((r["pathway"][:24], r["gene"], "KO" if r["direction"] == "0" else "UP") for r in lost)
+    out.append("only-MPB by perturbation: " + str(pert.most_common(15)))
+    return "\n".join(out)
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("arm", type=Path)
+    ap.add_argument("--mpb", type=Path, default=DEFAULT_MPB)
+    ap.add_argument("--out", type=Path)
+    a = ap.parse_args()
+    with (a.arm / "experimental_cases.tsv").open(newline="") as handle:
+        ours = list(csv.DictReader(handle, delimiter="\t"))
+    rows = split(load_mpb(a.mpb), ours)
+    print(report(rows))
+    if a.out:
+        with a.out.open("w", newline="") as handle:
+            w = csv.DictWriter(handle, fieldnames=list(rows[0]), delimiter="\t")
+            w.writeheader(); w.writerows(rows)
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

bench/analysis/mpbiopath_network_adapter.py

@@ -44,6 +44,7 @@
 
 import argparse
 import csv
+import re
 import sys
 from pathlib import Path
 
@@ -80,6 +81,10 @@ def catalog_ids(catalog: Path) -> dict[str, str]:
     mapping: dict[str, str] = {}
     for directory in sorted(p for p in catalog.iterdir() if p.is_dir()):
         name = directory.name
+        if re.fullmatch(r"R-HSA-\d+", name):
+            # Current generator naming: the directory IS the stable id.
+            mapping[name] = name
+            continue
         if "_R-HSA-" not in name:
             continue
         prefix, _, numeric = name.rpartition("_R-HSA-")
@@ -100,8 +105,11 @@ def case_pathway_names(cases_tsv: Path) -> dict[str, str]:
     """
     mapping: dict[str, str] = {}
     with cases_tsv.open(newline="") as handle:
-        for row in csv.DictReader(handle, delimiter="\t"):
-            pid = (row.get("pathway_id") or "").strip()
+        lines = (ln for ln in handle if not ln.startswith("#"))
+        for row in csv.DictReader(lines, delimiter="\t"):
+            # A cases dump (pathway_id, pathway_name) or the versioned
+            # pathway list bench/catalog_pathways.tsv (id, pathway_name).
+            pid = (row.get("pathway_id") or row.get("id") or "").strip()
             name = (row.get("pathway_name") or "").strip()
             if pid and name:
                 mapping[pid] = name
@@ -165,7 +173,7 @@ def main() -> int:
             missing_mpb.append(f"{prefix} ({stable})")
             continue
         edges = parse_network(source)
-        target = args.out / f"{prefix}_{stable}"
+        target = args.out / (stable if prefix == stable else f"{prefix}_{stable}")
         write_catalog_dir(target, edges)
         written.append((target.name, len(edges)))
 

bench/analysis/test_experimental_gap.py

@@ -0,0 +1,40 @@
+"""experimental_gap.py: the join must be exact and the quadrants right."""
+import sys
+from pathlib import Path
+
+import pytest
+
+sys.path.insert(0, str(Path(__file__).parent))
+from experimental_gap import report, split  # noqa: E402
+
+
+def case(p, g, d, k, pred, exp, valid="1", cat="pass"):
+    return {"pathway": p, "gene": g, "direction": d, "key_output": k, "predicted": pred,
+            "expected": exp, "valid": valid, "category": cat}
+
+
+MPB = {("P", "A", "0", "1"): ("0", "0"), ("P", "A", "0", "2"): ("2", "1"),
+       ("P", "B", "2", "1"): ("2", "2"), ("Q", "C", "0", "9"): ("1", "2")}
+
+
+def test_quadrants_and_invalid_counts_as_wrong():
+    rows = split(MPB, [case("P", "A", "0", "1", "2", "0"),            # only MPB
+                       case("P", "A", "0", "2", "2", "2"),            # only DS
+                       case("P", "B", "2", "1", "1", "2", valid="0"),  # invalid -> only MPB
+                       case("Q", "C", "0", "9", "0", "1")])           # neither
+    assert [r["quadrant"] for r in rows] == ["MPB/-", "-/DS", "MPB/-", "-/-"]
+    text = report(rows)
+    assert "MP-BioPath 2, DeltaSignal 1, gap +1" in text
+    assert "either 3" in text
+
+
+def test_a_label_disagreement_is_an_error():
+    with pytest.raises(SystemExit):
+        split(MPB, [case("P", "A", "0", "1", "0", "2")])
+
+
+def test_a_case_missing_on_either_side_is_an_error():
+    with pytest.raises(SystemExit):
+        split(MPB, [case("P", "A", "0", "1", "0", "0")])               # 3 MPB cases unmatched
+    with pytest.raises(SystemExit):
+        split({}, [case("P", "A", "0", "1", "0", "0")])

bench/analysis/test_mpbiopath_adapter.py

@@ -15,7 +15,7 @@
 
 sys.path.insert(0, str(Path(__file__).parent))
 from mpbiopath_network_adapter import (  # noqa: E402
-    CONJUNCTION, POLARITY, UUID_PREFIX, case_pathway_names, parse_network,
+    CONJUNCTION, POLARITY, UUID_PREFIX, case_pathway_names, catalog_ids, parse_network,
     write_catalog_dir,
 )
 
@@ -89,9 +89,22 @@ def check(label, actual, expected):
     check("case name bridge", case_pathway_names(cases),
           {"R-HSA-453279": "Mitotic_G1-G1_S_phases"})
 
+    # Current generator naming (bare stable-id directories) and the versioned
+    # pathway list as the name bridge; both were unreadable until 2026-10-09,
+    # so the adapter silently adapted nothing.
+    cat = tmp / "catalog"
+    (cat / "R-HSA-453279").mkdir(parents=True)
+    (cat / "Old_Name_R-HSA-69620").mkdir()
+    check("bare and legacy catalog dirs", catalog_ids(cat),
+          {"R-HSA-453279": "R-HSA-453279", "Old_Name": "R-HSA-69620"})
+    plist = tmp / "catalog_pathways.tsv"
+    plist.write_text("# comment\nid\tpathway_name\nR-HSA-453279\tMitotic_G1-G1_S_phases\n")
+    check("pathway list as name bridge", case_pathway_names(plist),
+          {"R-HSA-453279": "Mitotic_G1-G1_S_phases"})
+
 if failures:
     print("FAILED:")
     for f in failures:
         print("  -", f)
     sys.exit(1)
-print(f"all adapter assertions passed ({12} checks)")
+print(f"all adapter assertions passed ({14} checks)")