diff --git a/CLAUDE.md b/CLAUDE.md index 89681c3..7fc15ad 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -332,6 +332,19 @@ env vars only override for benchmark sweeps: `containment.csv`; the solve reports `self_inhibitor_rule` (on / off / inert) and how many inhibitor slots were flagged. Measured numbers are in specs/023. specs/022 has the rule; specs/023 has the adoption. +- `DS_DRUG_MODE=propagate` (default) | `inert` (specs/032): drug-derived + nodes (a Reactome `Drug`, a complex containing one, a set that is all drugs), + listed by the generator's `drugs.csv`, held at baseline like cofactors, for a + cell without the drug. The solve reports `drug_rule` and `drugs_held`. +- `DS_SET_POOL_MODE=product` (default) | `extreme` | `geomean` | `mean` | `max` + (specs/033): how a set POOL node (the generator's `LNG_SET_POOL`: one node + per set-valued catalyst or regulator, fed by `set_member` edges) combines its + members. It has no effect on networks without `set_member` edges. +- `DS_CONSERVED_MODE=off` (default) | `inert` (specs/035): a loop input released + UNCHANGED by its own loop (every producer consumes a complex containing it: + RAS GAPs, scaffold partners, F-actin) is held at baseline like a cofactor. + Signal carriers transformed inside the loop are not touched. The solve reports + `conserved_rule` and `conserved_held`. - Export aggregation default: `stoichiometry_weighted`. ### Where design decisions live @@ -413,6 +426,27 @@ for the last feature that touched it rather than re-deriving from the code. reaction forms (a curation-review list). 96 entity → reaction breaks are copy multiplicity; the targeted join is untried, near a family that failed three times. About 41 cases for an any-copy pool. Analysis only. +- `specs/032-inert-drugs/` — `DS_DRUG_MODE`: drug-derived nodes held at + baseline. Not adopted: held-out +30, but experimental −5, RAF worse. +- `specs/033-set-pool-nodes/`, `specs/038-pool-max/` — one OR node per + set-valued catalyst or regulator (`LNG_SET_POOL`), feeding the reaction once. + Adoptable on faithfulness under `DS_SET_POOL_MODE=product`: held-out +6 / −2, + experimental +3 / +8 (the +8 concentrated in 4 perturbations) on two + measurements. `max` (any member suffices) is refuted (held-out −120, + experimental −57); `mean` does not gain. +- `specs/034-regulator-wiring/` — inventory of how catalysts and regulators are + wired; the RAF collapse traced by iteration; MP-BioPath's network and its + 1,800 hand edits compared with ours. UUIDs already resolve 86% of their + loop artefacts. Our loops are mostly curated interconversion cycles, and the + 219 inhibitions they added are mostly the reversal enzyme of such cycles + (187 absent from ours). The design input for what comes next. +- `specs/035-conserved-participants/` — `DS_CONSERVED_MODE`: loop inputs + released unchanged held at baseline. Not adopted: it held PIP3 itself, and + experimental was −122. +- `specs/036-capped-bundles-as-pools/`, `specs/037-combined-raf-fixes/` — the + variant cap bundles a set's alternatives into one all-required node + (`LNG_CAP_POOLS` fixes the structure, default off). Combined with de-duplication + it was held-out −162; de-duplication is the damage. Not adopted. - `specs/009-solver-defaults/` — the one-variable-at-a-time re-measurement behind the `DS_*` defaults above (cited in that section too). - `.specify/memory/constitution.md` — project principles the specs are @@ -426,12 +460,12 @@ behaviour, not the whole directory; `ls specs/` is the complete list. Per-feature numbers belong in that feature's `research.md`, not here. ### Testing Strategy -**Most of the test suite cannot fail.** Thirteen files contain assertions; the +**Most of the test suite cannot fail.** Sixteen files contain assertions; the other seven execute code and print output. | file | assertions | note | |---|---|---| -| `test/test_config_validation.jl` | 196 | | +| `test/test_config_validation.jl` | 201 | | | `test/test_loop_elasticity.jl` | 173 | + 1 `@test_broken` | | `test/test_propagator_invariants.jl` | 121 | | | `test/test_loop_pool.jl` | 81 | | @@ -439,6 +473,9 @@ other seven execute code and print output. | `test/test_and_curves.jl` | 71 | | | `test/test_scc_break_roles.jl` | 56 | | | `test/test_self_inhibition.jl` | 80 | specs/022, added 2026-09-25 | +| `test/test_drugs.jl` | 59 | specs/032, added 2026-09-26 | +| `test/test_set_pool.jl` | 43 | specs/033, added 2026-09-26 | +| `test/test_conserved.jl` | 15 | specs/035, added 2026-09-26 | | `test/test_cli_observations.jl` | 39 | | | `test/test_api_errors.jl` | 44 | | | `test/test_cycle_handling.jl` | 34 | + 2 `@test_broken` | diff --git a/bench/analysis/regulator_inventory.py b/bench/analysis/regulator_inventory.py new file mode 100644 index 0000000..74d57d0 --- /dev/null +++ b/bench/analysis/regulator_inventory.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 +"""Inventory of how catalysts and regulators are wired into a catalog build +(specs/034). Every catalyst/regulator EDGE into a reaction node is classified on +four axes: + + role catalyst | positive | negative | requirement | pos_gene_expr | neg_gene_expr + (from the Regulation subclass in Neo4j; the generator only knows + catalyst / PositiveRegulation / NegativeRegulation) + shape of the curated entity: protein | complex | complex_with_set | set | + set_of_complexes | small_molecule | drug-derived (drugs.csv when present) + copy of the SOURCE node that is wired in: + root no incoming edge (what the benchmark pins) + upstream produced, and the reaction does not reach it + downstream produced, and the reaction reaches it (a loop closes) + pool an LNG_SET_POOL node + combine the edge's and_or flag + +With --cases, each curator/experimental case is attributed to the regulator +edges its perturbed gene's pinned nodes can reach through a regulator edge, and +accuracy is tallied per (role, shape, copy) cell. + +Usage: regulator_inventory.py --catalog BUILD_DIR [--cases ARM_DIR] +""" +from __future__ import annotations + +import argparse +import collections +import csv +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + +ROLE_Q = """ +UNWIND $rids AS rid MATCH (r:ReactionLikeEvent {stId: rid}) +CALL { + WITH r MATCH (r)-[:catalystActivity]->(:CatalystActivity)-[:physicalEntity]->(e) + RETURN 'catalyst' AS role, e.stId AS e + UNION WITH r MATCH (r)-[:regulatedBy]->(rg)-[:regulator]->(e) + RETURN CASE WHEN rg:Requirement THEN 'requirement' + WHEN rg:PositiveGeneExpressionRegulation THEN 'pos_gene_expr' + WHEN rg:NegativeGeneExpressionRegulation THEN 'neg_gene_expr' + WHEN rg:PositiveRegulation THEN 'positive' + ELSE 'negative' END AS role, e.stId AS e +} +RETURN rid AS r, role, e +""" + +SHAPE_Q = """ +UNWIND $ids AS s MATCH (e:PhysicalEntity {stId: s}) +RETURN s, + CASE WHEN e:SimpleEntity THEN 'small_molecule' + WHEN e:EntitySet AND size([(e)-[:hasMember|hasCandidate]->(m:Complex) | m]) > 0 THEN 'set_of_complexes' + WHEN e:EntitySet THEN 'set' + WHEN e:Complex AND size([(e)-[:hasComponent*1..6]->(x:EntitySet) | x]) > 0 THEN 'complex_with_set' + WHEN e:Complex THEN 'complex' + ELSE 'protein' END AS shape, + [(e)-[:hasComponent|hasMember|hasCandidate*0..8]->(l:PhysicalEntity) | l.stId] AS leaves +""" + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--catalog", type=Path, required=True) + ap.add_argument("--cases", type=Path, help="arm dir with {curator,experimental}_cases.tsv") + ap.add_argument("--cell", help="role/shape/copy: list its edges and the FAILED cases routed through them") + a = ap.parse_args() + import benchmark_vs_mpbiopath as BM + from py2neo import Graph + g = Graph(BM.NEO4J_URL, auth=(BM.NEO4J_USER, BM.NEO4J_PASSWORD)) + + names = {} + for line in open(Path(__file__).resolve().parents[1] / "catalog_pathways.tsv"): + if line.startswith(("#", "id\t")): + continue + sid, name = line.rstrip("\n").split("\t") + names[sid] = name + + cells = collections.Counter() # (role, shape, copy, and_or) -> edges + want = tuple(a.cell.split("/")) if a.cell else None + cell_edges = collections.Counter() # (pathway, regulator stid, reaction stid) -> failed cases + case_cells = collections.defaultdict(lambda: [0, 0]) # (role, shape, copy) -> [n, correct] + cases = {} + if a.cases: + for ax in ("curator", "experimental"): + p = a.cases / f"{ax}_cases.tsv" + if p.exists(): + for r in csv.DictReader(open(p), delimiter="\t"): + if r["valid"] == "1" and not r.get("exclusion"): + cases.setdefault(r["pathway"], []).append((ax, r)) + + for d in sorted(a.catalog.glob("R-HSA-*")): + pid = d.name + lab = {x["uuid"]: x["stable_id"] for x in csv.DictReader(open(d / "stid_to_uuid_mapping.csv"))} + kinds = {x["uuid"]: x["node_kind"] for x in csv.DictReader(open(d / "nodes.csv"))} + drugs = set() + if (d / "drugs.csv").exists(): + drugs = {x["stable_id"] for x in csv.DictReader(open(d / "drugs.csv"))} + E = list(csv.DictReader(open(d / "logic_network.csv"))) + fwd = collections.defaultdict(set) + indeg = collections.Counter() + for e in E: + fwd[e["source_id"]].add(e["target_id"]) + indeg[e["target_id"]] += 1 + reach_cache: dict = {} + + def reach(u): + if u not in reach_cache: + seen, st = {u}, [u] + while st: + for t in fwd[st.pop()]: + if t not in seen: + seen.add(t) + st.append(t) + reach_cache[u] = seen + return reach_cache[u] + + reg_edges = [e for e in E if e["edge_type"] in ("catalyst", "regulator")] + rx_stids = sorted({lab.get(e["target_id"], "") for e in reg_edges} - {""}) + roles = collections.defaultdict(list) # (reaction stid) -> [(role, entity stid)] + for r in g.run(ROLE_Q, rids=rx_stids).data(): + roles[r["r"]].append((r["role"], r["e"])) + ent = sorted({e for v in roles.values() for _, e in v}) + shape, leaves = {}, {} + for r in g.run(SHAPE_Q, ids=ent).data(): + shape[r["s"]], leaves[r["s"]] = r["shape"], set(r["leaves"]) + edge_cell = {} + for e in reg_edges: + src, tgt = e["source_id"], e["target_id"] + rstid, sstid = lab.get(tgt, ""), lab.get(src, "") + want_neg = e["pos_neg"] == "neg" + # the curated annotation this edge comes from: the entity equal to, + # or containing, the source node's entity, with a matching sign + cands = [(ro, en) for ro, en in roles.get(rstid, []) + if (ro in ("negative", "neg_gene_expr")) == want_neg + and (en == sstid or sstid in leaves.get(en, set()))] + if e["edge_type"] == "catalyst": + cands = [c for c in cands if c[0] == "catalyst"] or cands + else: + cands = [c for c in cands if c[0] != "catalyst"] or cands + role, en = cands[0] if cands else ("unmatched", sstid) + shp = "drug" if (en in drugs or sstid in drugs) else shape.get(en, "?") + if kinds.get(src) == "set_pool": + copy = "pool" + elif indeg[src] == 0: + copy = "root" + elif src in reach(tgt): + copy = "downstream" + else: + copy = "upstream" + cells[(role, shp, copy, e["and_or"])] += 1 + edge_cell[(src, tgt)] = (role, shp, copy) + # attribute cases: a case uses a regulator edge if its pins reach the source + for ax, r in cases.get(names.get(pid, ""), []): + pins = [u for u in r["gene_uuids"].split("|") if u] + if not pins: + continue + got = set() + seen, st = set(pins), list(pins) + while st: + for t in fwd[st.pop()]: + if t not in seen: + seen.add(t) + st.append(t) + outs = set(u for u in r["output_uuids"].split("|") if u) + ok = r["predicted"] == r["expected"] + for (s, t), c in edge_cell.items(): + if s in seen and outs & reach(t): + got.add(c) + if want and c == want and not ok: + cell_edges[(ax, names.get(pid, pid)[:30], r["gene"], r["direction"], + lab.get(s, s), lab.get(t, t))] += 1 + for c in got: + key = (ax,) + c + case_cells[key][0] += 1 + case_cells[key][1] += ok + + print("role / shape / copy / and_or -> edges") + for k, v in sorted(cells.items(), key=lambda kv: -kv[1]): + print(f" {v:7d} {' / '.join(k)}") + if want: + print(f"\nFAILED cases through {a.cell} (axis, pathway, gene, dir, regulator node, reaction):") + for k, v in cell_edges.most_common(25): + print(f" {v:4d} {k}") + return 0 + if case_cells: + print("\ncases whose route uses a regulator edge of the cell (a case may use many cells):") + for k, (n, ok) in sorted(case_cells.items(), key=lambda kv: -kv[1][0]): + if n >= 20: + print(f" {k[0]:12} {' / '.join(k[1:]):45} n={n:6d} acc={ok / n:.1%}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/bench/benchmark_mpbiopath_cases.py b/bench/benchmark_mpbiopath_cases.py index 9fe879a..873576e 100644 --- a/bench/benchmark_mpbiopath_cases.py +++ b/bench/benchmark_mpbiopath_cases.py @@ -260,6 +260,10 @@ def load_dbid_to_uuids(pathway_dir: Path) -> dict[str, list[str]]: if nodes_file.exists(): with nodes_file.open(newline="") as handle: for row in csv.DictReader(handle): + # specs/033: a set pool is not a readout or pin target (see + # benchmark_vs_mpbiopath.load_stid_to_uuids). + if normalize_cell(row.get("node_kind")) == "set_pool": + continue stable_ids = {normalize_cell(row.get("diagram_entity_id"))} stable_ids.update(normalize_cell(item) for item in normalize_cell(row.get("member_leaves")).split("|")) for stable_id in stable_ids - {""}: diff --git a/bench/benchmark_vs_mpbiopath.py b/bench/benchmark_vs_mpbiopath.py index 7eef9f6..33486fa 100644 --- a/bench/benchmark_vs_mpbiopath.py +++ b/bench/benchmark_vs_mpbiopath.py @@ -330,6 +330,14 @@ def load_stid_to_uuids(pathway_dir: Path): seen = set() with open(nodes_csv) as f: for row in csv.DictReader(f): + # specs/033: a set POOL node is the set as one catalyst/regulator, + # not a readout or a pin target. Indexed, it would capture set- + # valued readouts (p-ERK dimers, p-AKT) and join the mean of every + # member readout, changing what is MEASURED rather than the model + # (review of the set-pool change). Readouts and pins stay as in + # an unpooled build. + if (row.get("node_kind") or "") == "set_pool": + continue uuid = str(row["uuid"]) keys = set() de = (row.get("diagram_entity_id") or "").strip() @@ -889,10 +897,19 @@ def run_pathway(pathway_id: str, pathway_name: str, gene_to_stids_cache=None, # Record what was actually pinned, so a run's protocol is on the record # rather than inferred from its flags (specs/023: 89% of pins had silently # become mid-pathway complexes for two months). + # specs/032: a pin on a drug node would override its inert hold and carry + # the perturbation through the drug, so count them (expected 0). + drug_uuids = set() + drugs_csv = pathway_dir / "drugs.csv" + if drugs_csv.exists(): + with open(drugs_csv) as f: + for r in csv.DictReader(f): + drug_uuids.update(stid_to_uuids.get(r["stable_id"], [])) for us in gene_to_uuids.values(): PIN_TALLY["perturbations"] += bool(us) PIN_TALLY["pinned"] += len(set(us)) PIN_TALLY["pinned_roots"] += sum(1 for u in set(us) if indeg[u] == 0) + PIN_TALLY["pinned_drugs"] += sum(1 for u in set(us) if u in drug_uuids) # Resolve key_output dbIds to their REAL stIds (may be R-ALL-, R-NUL-, not # just R-HSA-) so species that are genuinely in the network are found. @@ -1002,7 +1019,10 @@ def run_pathway(pathway_id: str, pathway_name: str, gene_to_stids_cache=None, # bundle's cofactor list (the solve would report self_inhibitor_rule inert). network_payload = {"nodes": parsed["nodes"], "edges": edges, "pathways": parsed["pathways"], "containment": parsed.get("containment"), - "cofactor_stids": parsed.get("cofactor_stids")} + "cofactor_stids": parsed.get("cofactor_stids"), + # specs/032: without it an edge-drop arm under DS_DRUG_MODE=inert + # reached the solver with no drug table (review of the branch). + "drug_stids": parsed.get("drug_stids")} # Use the server-cached network by id (fast: send only observations per # solve). But if SKIP_EDGE_TYPES modified the edges above, the cached # network is stale, so send the full modified payload instead. @@ -1041,9 +1061,16 @@ def run_pathway(pathway_id: str, pathway_name: str, gene_to_stids_cache=None, raise SystemExit(f"{pathway_name}: the solve reports self_inhibitor_rule=" f"{ds_result['self_inhibitor_rule']!r} -- the network reached the solver " "without its containment table, so this is not the default model.") + if ds_result and ds_result.get("drug_rule") == "inert: no drug table": + raise SystemExit(f"{pathway_name}: DS_DRUG_MODE=inert but the network has no drugs.csv " + "(a build that predates it), so no drug is held (specs/032).") if ds_result: SOLVE_TALLY["solves"] += 1 SOLVE_TALLY["converged"] += bool(ds_result.get("converged")) + SOLVE_TALLY[f"drug_rule={ds_result.get('drug_rule', 'unknown')}"] += 1 + SOLVE_TALLY["drugs_held"] += int(ds_result.get("drugs_held") or 0) + SOLVE_TALLY[f"conserved_rule={ds_result.get('conserved_rule', 'unknown')}"] += 1 + SOLVE_TALLY["conserved_held"] += int(ds_result.get("conserved_held") or 0) activities = ds_result.get("node_activities", {}) if ds_result else {} for r in rows: @@ -1330,8 +1357,18 @@ def main(): f"{r['total']}\t{r['correct']}\t{r['accuracy']:.6f}\t" f"{r['valid_total']}\t{r['valid_correct']}\t{r['valid_accuracy']:.6f}\n") print(f"\nConverged: {SOLVE_TALLY['converged']} of {SOLVE_TALLY['solves']} solves") + rules = sorted(k.split("=", 1)[1] for k in SOLVE_TALLY if k.startswith("drug_rule=")) + print(f"\nDrugs: rule {'/'.join(rules) or 'unknown'}, {SOLVE_TALLY['drugs_held']} node-solves held (specs/032)") + if "inert" in rules and SOLVE_TALLY["drugs_held"] == 0: + # An inert arm that held nothing measured the default model. + raise SystemExit("DS_DRUG_MODE=inert but no solve held a drug node; the arm is the control.") + crules = sorted(k.split("=", 1)[1] for k in SOLVE_TALLY if k.startswith("conserved_rule=")) + print(f"\nConserved: rule {'/'.join(crules) or 'unknown'}, {SOLVE_TALLY['conserved_held']} node-solves held (specs/035)") + if "inert" in crules and SOLVE_TALLY["conserved_held"] == 0: + raise SystemExit("DS_CONSERVED_MODE=inert but no solve held a conserved node; the arm is the control.") print(f"\nPinned: {PIN_TALLY['pinned']} nodes over {PIN_TALLY['perturbations']} " - f"perturbations, {PIN_TALLY['pinned_roots']} of them roots (DS_PIN_SCOPE={PIN_SCOPE})") + f"perturbations, {PIN_TALLY['pinned_roots']} of them roots, " + f"{PIN_TALLY['pinned_drugs']} drug nodes (DS_PIN_SCOPE={PIN_SCOPE})") print(f"\nPer-pathway report: {args.report}") if args.dump_cases: diff --git a/bench/test_analysis_numbers.py b/bench/test_analysis_numbers.py index 5c94334..29f8196 100644 --- a/bench/test_analysis_numbers.py +++ b/bench/test_analysis_numbers.py @@ -438,3 +438,18 @@ def test_drop_carriers_removes_small_molecule_links_both_ways(self): pruned2 = self.co.drop_carriers(adj, lambda n: n.startswith("R-HSA-") and n != "R-HSA-rx2") assert "R-HSA-B" not in self.co.signed_reach(pruned2, "R-HSA-A") + + +def test_set_pool_nodes_are_not_readouts_or_pin_targets(tmp_path): + # specs/033: indexing a pool under its set id and member leaves would move + # set-valued readouts onto the pool and add it to member readouts. + (tmp_path / "nodes.csv").write_text( + "uuid,node_kind,diagram_entity_id,compartment,member_leaves,source_sets,chosen_members\n" + "u-a,simple_entity,R-HSA-A,,R-HSA-A,,\n" + "u-pool,set_pool,R-HSA-SET,,R-HSA-A|R-HSA-B,,\n") + m = bench.load_stid_to_uuids(tmp_path) + assert m["R-HSA-A"] == ["u-a"] + assert "R-HSA-SET" not in m and "R-HSA-B" not in m + import benchmark_mpbiopath_cases as cases + d = cases.load_dbid_to_uuids(tmp_path) + assert all("u-pool" not in v for v in d.values()) diff --git a/cli/deltasignal.jl b/cli/deltasignal.jl index dfacf6b..80e97f5 100755 --- a/cli/deltasignal.jl +++ b/cli/deltasignal.jl @@ -203,6 +203,8 @@ function execute_parse_command(args) "cofactor_stids" => collect(network.cofactor_stids), # Same reason, for the self-inhibitor rule (specs/022). "containment" => network_containment_json(network), + # And for DS_DRUG_MODE (specs/032); null when the bundle has no list. + "drug_stids" => drug_stids_json(network), "set_mappings" => Dict(set_id => Dict( "original_set_id" => mapping.original_set_id, "original_name" => mapping.original_name, @@ -379,7 +381,8 @@ function execute_solve_command(args) cofactor_stids = haskey(network_json, "cofactor_stids") ? Set(String.(network_json["cofactor_stids"])) : Set{String}() containment = containment_from_json(get(network_json, "containment", nothing)) - network = ReactionNetwork(nodes, edges, set_mappings, cofactor_stids, containment) + drug_stids = drug_stids_from_json(get(network_json, "drug_stids", nothing)) + network = ReactionNetwork(nodes, edges, set_mappings, cofactor_stids, containment, drug_stids) println("✓ Network loaded: $(length(nodes)) nodes, $(length(edges)) edges") # Load observations CSV diff --git a/docs/API.md b/docs/API.md index cad0370..1495810 100644 --- a/docs/API.md +++ b/docs/API.md @@ -101,7 +101,9 @@ Response: ``` `containment` (restricted to stable ids in this network) and `cofactor_stids` are what the solver needs to run the default model; send them back with an -inline network. +inline network. `drug_stids` (specs/032, additive) lists the bundle's +drug-derived stable ids, or is `null` when the bundle has no `drugs.csv`; send +it back too, or `DS_DRUG_MODE=inert` reports `"inert: no drug table"`. Node display names are enriched from the Reactome ContentService when the generator only provided stable ids (degrades gracefully if that service is @@ -149,6 +151,10 @@ Response: - **`node_activities`**: `uuid → activity` on the **0–1 internal scale** (`× 100` for display — see Conventions). - **`self_inhibitors`** (additive, specs/022): the number of inhibitor slots FLAGGED as self-contained (the inhibitor contains one of its own reaction's inputs AND is reachable from it in the network). A flagged slot is damped only when the inhibitor moves with that input. On by default (`DS_SELF_INHIBITOR_WEIGHT=0.1`; `off` disables it). +- **`drug_rule`** / **`drugs_held`** (additive, specs/032): `"propagate"` + (default), `"inert"`, or `"inert: no drug table"`, and how many drug-derived + nodes this solve held at baseline. Under `inert` a drug participates at fold + 1.0 and never carries a perturbation; an explicit observation still wins. - **`self_inhibitor_rule`** (additive): `"on"`, `"off"`, or `"inert: no containment table"`. The last means the solve did NOT run the default model, because the network came without a containment table. `/api/parse` returns `containment`; send it back with a POSTed network to keep the rule active. - **`scc`** (additive, specs/017): how the cyclic components were resolved — `method` (the `DS_SCC_METHOD` in force), `pooled`, `iterated`, diff --git a/scripts/run_arm.sh b/scripts/run_arm.sh index 369f2b0..8f6c9b7 100755 --- a/scripts/run_arm.sh +++ b/scripts/run_arm.sh @@ -156,6 +156,8 @@ for gt in ("curator", "experimental"): m = re.search(r"^Protocol: .*$", text, re.M); rec[f"{gt}_protocol"] = m.group(0) if m else None m = re.search(r"^Pinned: .*$", text, re.M); rec[f"{gt}_pinned"] = m.group(0) if m else None m = re.search(r"^Converged: .*$", text, re.M); rec[f"{gt}_converged"] = m.group(0) if m else None + m = re.search(r"^Drugs: .*$", text, re.M); rec[f"{gt}_drugs"] = m.group(0) if m else None + m = re.search(r"^Conserved: .*$", text, re.M); rec[f"{gt}_conserved"] = m.group(0) if m else None json.dump(rec, open(os.path.join(out, "ARM.json"), "w"), indent=2) print(json.dumps(rec, indent=2)) PY diff --git a/specs/032-inert-drugs/research.md b/specs/032-inert-drugs/research.md new file mode 100644 index 0000000..0f3a6b0 --- /dev/null +++ b/specs/032-inert-drugs/research.md @@ -0,0 +1,146 @@ +# specs/032 — Drug entities as inert participants + +**Status:** pre-registered 2026-09-26, before any arm has run. + +## Why + +The experimental axis is where we trail MP-BioPath's published result: +- ours 67.14% (570/849); +- theirs 75.74%; +- paired against their supplementary table, 67.4% against 73.8%, −54 cases. + +**One pathway carries 32 of the 54:** RAF/MAP kinase cascade (R-HSA-5673001), +where we score 28.6% and they score 93.9%. On MP-BioPath's own networks +DeltaSignal scores RAF at 100%, so the cause is our network, not the solver. + +**Traced (canonical build `20260926-1221_590301c`, HRAS up):** +- The baseline is healthy: 1,571 of 1,571 nodes at 1.0, converged in 5 + iterations. +- Under HRAS up the solve does not converge (535 iterations), and 127 of the + 128 nodes of the second cyclic component read 0. +- That component's only negative edges are 4 regulator edges from + **drug-bound complexes**: "activated RAF:scaffold:MAP2K:MAPK complex:dual + (single) mechanism MAP2K inhibitors". They inhibit "RAF phosphorylates + MAP2K dimer" and "MAP2Ks phosphorylate MAPKs". +- The drug is a root held at baseline, so raising RAS or RAF raises the + drug-bound complexes, which shut the kinase steps off. +- Every RAS/RAF perturbation reads DOWN, including the NF1 knockdown, which + should read UP. +- Reactome 97 curates these drug actions inside the pathway. MP-BioPath's 2019 + networks predate them. + +**Scope in the catalog:** +- 202 drug entities, or entities that contain a drug: 134 ChemicalDrug, 29 + ProteinDrug, 39 complexes. +- They occur in 24 of 92 pathways: RAF 37, DNA Repair 28, Fanconi 24, Mitotic + G1 24, MET 23, PDGF 19, KIT 18, IFN α/β 14, ERBB2 14, and others. +- No benchmark pin and no readout is a drug entity (checked on both axes). + +## The rule + +A benchmark case describes a cell **without the drug**. So an entity that is, +or contains, a Reactome `Drug` is held at baseline: +- it stays a participant (an AND input of fold 1.0; an inhibitor at baseline + divides by 1); +- it can never carry or amplify a perturbation. + +This is the cofactor rule (specs/007) applied to drugs. The curated edges stay +in the network, which keeps faithfulness to Neo4j. Pinning at 0 is rejected: +under `divide` inhibition, a regulator at 0 de-represses its target tenfold, +which models drug *withdrawal*, not drug absence. An explicit observation +still wins. + +**Implementation:** +- **Generator:** ships `drugs.csv` beside `cofactors.csv`, listing the + pathway's drug entities from the release's own data. It is an additive file; + the network is byte-identical. +- **Solver:** `DS_DRUG_MODE` = `propagate` (default, byte-identical) or + `inert`. The solve reports how many nodes were held, and the benchmark + refuses an `inert` arm whose solves hold none. + +## Pre-registration + +**Arms:** +- `drugctrl`: code defaults. +- `druginert`: `DS_DRUG_MODE=inert`. +- Both on the same build, the first one that ships `drugs.csv`, through + `scripts/run_arm.sh`. + +**Adopt only if all three hold:** +- experimental net ≥ +10 with p < 0.05; +- curator held-out net ≥ −15 (the regeneration noise floor); +- no pathway loses more than 10 curator held-out cases. + +The rule has no parameter, and nothing is tuned on either axis. + +**Predictions, written before the arm:** +- RAF experimental rises by at least 20 of its 49 cases; +- RAF curator improves; +- the other drug pathways move by less than 10 each; +- pathways without drugs are byte-identical. + +**Also reported:** +- curator tuning; +- per-pathway net on both axes; +- concentration: pathways moved, and distinct perturbations; +- the convergence count for RAF under each arm; +- the paired gap against MP-BioPath's supplementary table. + +## Result — NOT ADOPTED (2026-09-26) + +**Setup:** +- Build `20260926-1359_f2ba9cc_drugs`; both arms at solver 093282f, through + `scripts/run_arm.sh`. +- `drugctrl`: code defaults, which report drug rule `propagate` with 0 held. +- `druginert`: `DS_DRUG_MODE=inert`; it held 6,940 drug node-solves on curator + and 3,016 on experimental. +- Pins are identical between the arms (1,269 nodes over 864 perturbations), and + so is convergence (1,659 of 1,725 on curator, 212 of 244 on experimental). + +| split | net | fixed / broke | p | pathways | +|---|---|---|---|---| +| curator held-out | **+30** | 30 / 0 | 1.9e-09 | 4 (+4 / −0): MET +10, ROCKs +8, KIT +6, VEGF +6 | +| curator tuning | −16 | 5 / 21 | 0.0025 | RAF −17, ERBB2 +1 | +| **experimental** | **−5** | 3 / 8 | 0.23 | 1 (RAF) | + +**Gates:** +- experimental must be at least +10 with p < 0.05. It is −5: **FAIL.** +- held-out must be no worse than −15. It is +30: pass. +- no pathway may lose more than 10 held-out cases: pass. + +**The predictions failed.** RAF experimental did not rise by at least 20: it +fell. RAF curator did not improve: it is −17. **Not adopted.** `DS_DRUG_MODE` +stays, default `propagate`, as the record. + +### Why the mechanism in this spec was incomplete + +HRAS up in RAF (R-HSA-5673001), solved on the drugs build: +- `propagate`: 537 iterations, not converged, 287 nodes at 0. +- `inert` (40 drug nodes held): still not converged, 235 at 0. +- In the 128-node component, 127 nodes are at 0 under `propagate` and 124 under + `inert`. + +**The drug-bound complexes' inhibitory edges were real but not the cause.** +The component is a **positive loop closed through catalyst edges**: +- "Phosphorylation of RAF" (80 variant copies) takes the phosphorylated MAP2K + dimers (p-S218,S222 MAP2K1 dimer; p-S222,S226 MAP2K2 dimer; the p-2S + MAP2K1:MAP2K2 heterodimer) as **catalysts**. +- Inside this component those dimers are produced only by "Dissociation of + RAS:RAF complex". +- That reaction consumes the activated RAF:scaffold:p-2S MAP2K:p-2T MAPK + complex, which is made downstream of RAF phosphorylation. +- "MAP2Ks and MAPKs bind to the activated RAF complex" consumes F-actin, CNKSR2 + and Ca2+, which in this component come back only from the same dissociation. +- The external inputs are healthy: p21 RAS:GTP:RAF complex is at 100x, and ATP + and the kinases are at 1x. +- Pushed hard, the loop does not settle, and it ends in the all-zero state. + +This is the loop knife-edge (specs/013, 014) on a component built by +**catalyst and dissociation recycling**. Holding the drugs removes 3 of the +127 zeros. + +**Open, not traced:** +- whether the catalyst edges (p-MAP2K → Phosphorylation of RAF) match what + Reactome curates for that reaction, or were derived; +- why the held-out +30 (MET, ROCKs, KIT, VEGF) moved; the drug rule is sound + there, but it was not gated separately. diff --git a/specs/033-set-pool-nodes/research.md b/specs/033-set-pool-nodes/research.md new file mode 100644 index 0000000..09e0497 --- /dev/null +++ b/specs/033-set-pool-nodes/research.md @@ -0,0 +1,230 @@ +# specs/033 — One node per set-valued catalyst or regulator + +**Status:** pre-registered 2026-09-26, before any code or arm. + +## The defect + +A Reactome catalyst or regulator can be an **EntitySet**, meaning any one of +these members plays the role. An example is "RAF activating kinases" +(R-HSA-5672734), the catalyst of "Phosphorylation of RAF": JAK2, SRC, BRAF, +RAF1, ARAF, CaMKII and three phosphorylated MAP2K dimers. + +**What the generator does now** (`append_regulators`, subunit/member +decomposition, a design chosen 2026-07-15 because bundling lost single-gene +propagation): +- It attaches **every member as its own edge on every reaction copy**, marked + AND for catalysts and positive regulators and OR for negative regulators. +- Nothing is copied per member. + +**So the reaction node computes** +`inputs × JAK2 × SRC × … × MAP2K dimer`: +- two members up multiply; +- one member at 0 kills the reaction; +- a member made downstream of the reaction makes the reaction require its own + product. + +**Traced in RAF/MAP kinase (specs/032):** the three MAP2K dimers are made only +by "Dissociation of RAS:RAF complex", downstream of RAF phosphorylation. The +component self-requires, does not converge, and ends at 0 under every RAS/RAF +perturbation (experimental 28.6%, against MP-BioPath's 93.9%). + +**Inventory (build `20260926-1221_590301c`, groups = reaction copy × set):** + +| role | groups | current member flag | +|---|---|---| +| catalyst | 1,881 | AND | +| positive regulator | 56 | AND | +| negative regulator | 146 | OR | + +- 63 pathways have at least one catalyst group. +- These edges close only 201 of 8,064 cyclic nodes (RAF 84). + +**What the ground truth says** when the perturbed gene is itself a member of a +catalyst set: + +| axis | cases | accuracy | expected change | +|---|---|---|---| +| curator | 1,030 | 70.7% (KO) / 63.9% (UP) | 68% | +| experimental | 53 | 57.7% (KO) / 29.6% (UP) | 96–100% | + +For comparison, other genes score 85–87% on curator, with 38% expected change. +**So a single member is expected to matter.** A plain mean over members (which +makes one member of nine nearly invisible) contradicts the ground truth. The +failures are mostly missed changes (curator 125 DOWN read as NORMAL, 133 UP +read as NORMAL), 208 of them in Mitotic G1. + +## The design + +**Generator (`LNG_SET_POOL=1`):** for a bare EntitySet catalyst or regulator +(not a modifier-isoform set): +- emit **one pool node per set per pathway**, with the set's stable id, which + is the entity the diagram draws; +- add a `set_member` edge (pos, OR) from each member to the pool, emitted once; +- add **one** edge from the pool to each reaction copy, in the original role + and flag. + +Every member reaches every reaction the set serves, exactly as before, so +reachability and cyclic structure are unchanged. Complexes are untouched: +under complex-as-node they are already one node. + +**Solver (`DS_SET_POOL_MODE`):** a node whose incoming edges are all +`set_member` is a pool, and its fold is combined from its members' folds. + +| mode | combination | one member KO | one member up 80x | several members up | +|---|---|---|---|---| +| `product` | ∏ folds, capped as hill_sat is | 0 | 80x | multiply | +| `extreme` | the member with the largest \|log fold\| (a tie goes to the lower fold) | 0 | 80x | the largest | +| `geomean` | (∏ folds)^(1/n) | 0 | 80^(1/n) | damped | +| `mean` | arithmetic mean (the existing OR rule) | (n−1)/n | ≈80/n | diluted | + +- `product` on pooled networks is close to today's arithmetic for catalysts + and positive regulators. +- For negative regulators the pool is a single inhibitor term, where today + each member divides separately. +- Default: `product` until adopted. + +## Pre-registration + +**Arms,** all through `scripts/run_arm.sh`, on one pooled build plus the +canonical build: +- `base`: the canonical build `20260926-1221_590301c`, code defaults; +- `pool_product`, `pool_extreme`, `pool_geomean`, `pool_mean`: the pooled + build with each mode. + +**Choosing the mode (no held-out data used):** +- the winner is the mode with the best **tuning-split** curator macro-F1 plus + experimental macro-F1 (each on its own case set; summed); +- ties go to the simpler mode, in the order product < extreme < geomean < mean. + +**Adopt** the pooled generator plus the chosen mode against `base` only if all +of these hold: +- curator held-out net ≥ −15 (the regeneration noise floor) and not + significantly negative (p ≥ 0.05, or net ≥ 0); +- experimental net ≥ 0; +- tuning curator net ≥ −15. + +It is adopted on faithfulness (one node per curated participant) as long as +it is not harmful. That is the same basis as specs/020's variant sharing. + +**Predictions:** +- RAF experimental rises under `extreme` and `geomean` (the loop stops + requiring its own product); +- `pool_product` moves few predictions against `base`, apart from the TP53 + uuid draw; +- `pool_mean` loses set-member cases on both axes. + +**Also reported:** +- per-pathway net on both axes; +- concentration; +- convergence counts; +- the set-member-case accuracy table above for each arm; +- cyclic node count (it must be unchanged); +- node and edge counts. + +## Pre-registration amendment (after review, before any arm) + +An adversarial review of both code changes found three things the +pre-registration did not state. They are fixed or specified here, **before +any arm has run**: + +1. **Readouts and pins are unchanged.** A pool node carries the set's stId and + its members, and the benchmark indexed every node under both. + - A set-valued readout that is also a catalyst set (p-T,Y MAPK dimers, + p-AKT, GSK3, p-MAPK8/9/10) would have resolved to the pool instead of + the members' mean. + - Every member readout would have gained the pool in its mean. + - Both benchmark loaders now skip `set_pool` rows (test in + `bench/test_analysis_numbers.py`), so the arms measure the model, not a + change of readout. +2. **Depletion.** + - A pooled phosphatase or ubiquitin-ligase set now depletes its substrate + through its pool: one edge instead of one per member. Under `product` that + is the same fold; under `extreme` or `geomean` the depletion strength + follows the mode. This is part of the design and is measured with it. + - The "a catalyst does not deplete itself" guards compared stIds, and a + pool's stId is the set's, so a pool would have depleted its own member. + This happens in 4 E3-ligase sets and would have closed new cycles. + - Fixed in LNG 81499ef: a pool never depletes one of its members. +3. **The cycle gate** is restated as **cyclic nodes excluding pool nodes must be + unchanged**. A pool on an existing loop (RAF's is the motivating one) is + itself cyclic. + - On the first pooled build: 8,078 cyclic nodes, of which 14 are pools, so + 8,064 non-pool, equal to the base. + - It is re-checked on the rebuilt catalog. + +**Known approximations of `product` against today:** +- **Caps:** the pool caps the member product at 100x before the reaction + multiplies. Two members at 20x read 100x, not 400x. The difference appears + only above the cap. +- **Negative-regulator sets become one divide term.** It is equal to the old + per-member divides below the cap. Under `DS_INHIBITOR_OR=1` (not the + default) min over members becomes 1/product. +- **The self-inhibitor rule** (specs/022) can now flag a pooled inhibitor set + whose containment includes the reaction's input, where before it flagged + only the containing member. + +**Build:** the arms run on a pooled build regenerated at LNG 81499ef, not on +the first pooled build. + +## Result — the pooled structure ADOPTABLE with `product` (2026-09-26) + +**Setup:** +- Base: canonical `20260926-1221_590301c`. Pooled: `20260926-2229_81499ef_setpool2` + (299 pools, 1,911 set_member edges, 64 depletion edges from pools). +- Non-pool cyclic nodes: 8,064, equal to the base. **Gate 3 passes.** +- All arms at solver 039ce74, through `run_arm.sh`. Pins are identical across + arms (1,269 nodes over 864 perturbations). + +**Choosing the mode** (tuning curator macro-F1 plus experimental macro-F1, as +pre-registered): + +| mode | tuning curator mF1 | experimental mF1 | sum | +|---|---|---|---| +| **product** | 0.7558 | 0.5933 | **1.3491** | +| extreme | 0.7557 | 0.5933 | 1.3490 | +| mean | 0.7401 | 0.5952 | 1.3353 | +| geomean | 0.7341 | 0.5933 | 1.3274 | + +**`product` against base:** + +| split | net | fixed / broke | p | gate | +|---|---|---|---|---| +| curator held-out | +6 | 10 / 4 | 0.18 | ≥ −15, not significantly negative: **pass** | +| curator tuning | +10 | 10 / 0 | 0.002 | ≥ −15: **pass** | +| experimental | +3 | 3 / 0 | 0.25 | ≥ 0: **pass** | + +Held-out for the other modes is reported only; none of it was used to +choose: extreme +9 (p = 0.035); geomean −20 (TP53 −79, IFN-γ −17); mean −6 +(TP53 −68). + +**Verdict:** the pooled structure plus `product` passes every gate. It is +**adoptable on faithfulness**: one node per curated catalyst or regulator, +the set node the diagram draws, and 10,950 fewer edges, with no measurable +harm. + +**Predictions:** +- `pool_product` moves few predictions: **held** (24 of 23,511 curator cases). +- RAF experimental rises under extreme or geomean: **only partly**: + +| arm | RAF experimental (of 49) | RAF curator (of 84) | experimental predictions DOWN / NORMAL / UP | +|---|---|---|---| +| base | 14 | 40 | 34 / 12 / 3 | +| product | 17 | 43 | 31 / 12 / 6 | +| extreme | 17 | 43 | 31 / 12 / 6 | +| geomean | 19 | 46 | 29 / 12 / 8 | +| mean | 22 | 48 | 26 / 12 / 11 | + +MP-BioPath scores RAF experimental 46 of 49. **The collapse has a second +cause.** "MAP2Ks and MAPKs bind to the activated RAF complex" takes F-actin, +CNKSR2 and Ca2+ as inputs, and in this component they are produced only by +"Dissociation of RAS:RAF complex", downstream. The input is wired to the +recycled, downstream copy. This is the same wrong-copy pattern as specs/030's +joins, and it is traced next. + +## Note added at merge (2026-09-27) + +Under the hierarchy (specs/030), a root complex with a set component now joins +the set's POOL node (one assembly edge: pool → complex) instead of the member +leaves. So `DS_SET_POOL_MODE` also shapes that assembly. Under `product` the +result is numerically close to the member join. It is included in every pooled +build measured here and in specs/038. diff --git a/specs/034-regulator-wiring/research.md b/specs/034-regulator-wiring/research.md new file mode 100644 index 0000000..be6f04c --- /dev/null +++ b/specs/034-regulator-wiring/research.md @@ -0,0 +1,408 @@ +# specs/034 — How catalysts and regulators are wired: inventory + +**Status:** analysis, a design input. Nothing is proposed for adoption here. + +**Why (Adam, 2026-09-26):** "we might need to put a lot more thought to how +these regulators are connected ... there are a lot of potential +configurations." Four findings in a row traced to regulator wiring: +- set members wired as separate required inputs (specs/033); +- an arbitrary copy chosen for joins (specs/030, STAT1); +- drug complexes as inhibitors (specs/032); +- recycled catalysts closing loops (RAF). + +Before designing, this inventory counts which configurations occur and which +cases route through them. + +**Tool:** `bench/analysis/regulator_inventory.py --catalog BUILD --cases ARM`. +- Build `20260926-1221_590301c` (canonical); cases from arm + `results/039ce74/poolbase`. +- The role comes from the Regulation subclass in Neo4j. + +## 1. Annotations (Neo4j, 92 pathways) + +| role | protein | complex | complex with set | set | set of complexes | small molecule | total | +|---|---|---|---|---|---|---|---| +| catalyst | 263 | 406 | 427 | 156 | 108 | 0 | 1,361 | +| positive gene-expression regulation | 42 | 141 | 120 | 4 | 4 | 1 | 312 | +| negative regulation | 46 | 59 | 108 | 13 | 9 | 2 | 237 | +| positive regulation | 28 | 50 | 38 | 15 | 4 | 12 | 147 | +| negative gene-expression regulation | 4 | 53 | 64 | 1 | 3 | 0 | 125 | +| requirement | 3 | 1 | 10 | 1 | 1 | 9 | 25 | + +**The generator has three roles** (`_CATALYST_CYPHER`, `_POS_REG_CYPHER`, +`_NEG_REG_CYPHER`), and Neo4j's labels fold the rest into them: +`Requirement` and `PositiveGeneExpressionRegulation` both carry +`PositiveRegulation`, and `NegativeGeneExpressionRegulation` carries +`NegativeRegulation`. So a Requirement ("needed, not sufficient") is wired +exactly like a positive regulator: an AND input whose rise multiplies the +reaction. + +## 2. Edges in the network, by role and wired copy + +`copy` is the SOURCE node of the edge: +- **root:** no incoming edge; +- **upstream:** produced, and not reached from the reaction; +- **downstream:** produced by something the reaction itself reaches, so the edge + closes a loop. + +| role | edges | root | upstream | downstream | +|---|---|---|---|---| +| catalyst | 26,457 | 20% | 48% | **32%** (8,481) | +| negative | 1,954 | 16% | 82% | 2% | +| positive | 1,948 | 16% | 83% | 1% | +| positive gene expression | 409 | 8% | 88% | 5% | +| negative gene expression | 158 | 4% | 95% | 1% | +| requirement | 118 | 29% | 42% | 29% | + +**A third of all catalyst edges are fed by a copy the reaction regenerates +downstream.** That is catalytic recycling (the loop taxonomy of 2026-06-12), +and it is rare for regulators. By shape: + +| catalyst shape | edges from a downstream copy | +|---|---| +| set | 6,575 | +| complex with set | 1,435 | +| set of complexes | 369 | +| complex | 89 | +| protein | 13 | + +The cells with the most edges are catalyst / complex with set / upstream +(7,210), catalyst / set / downstream (6,575) and catalyst / set / root (3,124). + +## 3. Cases routed through each configuration + +A case counts toward a cell if its pinned nodes reach the edge's source and the +edge's reaction reaches its readout. **A case can count in many cells, so these +accuracies are associations, not causes.** They point at where to trace. + +**Low cells (n ≥ 40):** + +| axis | configuration | n | accuracy | +|---|---|---|---| +| curator | positive / complex with set / downstream | 240 | **42.9%** | +| curator | catalyst / set of complexes / root | 30 | 40.0% | +| curator | negative / protein / downstream | 42 | 35.7% | +| curator | catalyst / set of complexes / downstream | 316 | 55.7% | +| curator | positive gene expression / complex with set / downstream | 280 | 58.2% | +| curator | requirement / complex with set / upstream | 46 | 63.0% | +| curator | negative / complex / upstream | 312 | 65.1% | +| curator | negative / complex with set / upstream | 948 | 67.8% | +| experimental | catalyst / set of complexes / downstream | 43 | **20.9%** | +| experimental | negative gene expression / complex with set / upstream | 56 | 37.5% | +| experimental | negative / complex / upstream | 62 | 46.8% | +| experimental | negative / complex with set / upstream | 87 | 48.3% | +| experimental | positive gene expression / complex with set / upstream | 97 | 56.7% | + +**High cells**, for contrast: catalyst / set / upstream (curator 91.3%, n = 738), +catalyst / set of complexes / upstream (experimental 96.4%, n = 165), +negative / set / root (curator 100%, n = 60). + +## 4. Mismatches between wiring and Reactome's meaning + +| # | mismatch | size | status | +|---|---|---|---| +| M1 | Requirement wired as a multiplying positive regulator | 25 annotations, 118 edges | not tested | +| M2 | set-valued catalysts and regulators as separate AND/OR member terms | 1,881 + 56 + 146 groups | specs/033: pooling is neutral-to-positive; draft LNG #99 | +| M3 | catalyst fed by a downstream (recycled) copy | 8,481 edges, 32% of catalysts | `DS_SCC_BREAK_CATALYST` (TP53 +104, held-out −65) and specs/018 roles: not adopted | +| M4 | inputs fed by a downstream copy (RAF F-actin, CNKSR2, Ca2+) | not yet counted | specs/033 result: traced, open | +| M5 | drug-derived regulators propagate | 202 entities, 24 pathways | specs/032: inert not adopted | +| M6 | negative regulators that are complexes (with or without sets) | the low cells in §3 | not traced | +| M7 | gene-expression regulation wired like reaction regulation | 437 annotations | semantics plausible; low experimental cells, not traced | + +**Next:** trace one case from each low cell in §3 (M6 and M7 first, because they +are the least explained), and count M4. Then propose a rule set for review. + +## 5. Traced: M6 (negative regulators that are complexes) + +`regulator_inventory.py --cell negative/complex/upstream`. The failed cases +concentrate on five regulator → reaction pairs: + +| pathway, perturbation | inhibitor | inhibited reaction | mechanism | +|---|---|---|---| +| RAF (KRAS/HRAS/NRAS up, NF1 KO) | BRAP:KSR1:MARK3 | "RAS:GTP:'activator' RAF homo/heterodimerizes …" | inside the RAF loop (see §6) | +| PIP3, PDPK1 KO | **N:M:PDPK1**: PDPK1 bound by **SARS-CoV-2** N and M | "PDPK1 phosphorylates AKT at T308" | self-contained, so the KO also removes the inhibitor and cancels itself; the virus is exogenous | +| ERBB2 up | **ERBB2:trastuzumab**:ERBIN:HSP90:CDC37 | "Trans-autophosphorylation of ERBB2 heterodimers" | self-contained and a **drug** | +| Activin (INHBA, ACVR2A) | INHIBIN-A:TGFBR3:ACVR2A | "Activin … binds Activin Receptor ACVR2A,B:ACVR1B" | self-contained sequestration (specs/022) | +| WNT (WNT1, WNT5A) | pT298-NLK dimer | "TCF/LEF:CTNNB1 bind canonical WNT target promoters" | not traced | + +**Exogenous entities** are not in the cell a benchmark case describes: + +| kind | entities | pathways | +|---|---|---| +| drug-derived | 202 | 24 | +| contains a non-human-species protein | 66 | 9 | + +The pathogens are Influenza A, RSV, Rotavirus, HCV, Measles, SARS-CoV-2, +N. meningitidis and C. trachomatis. DDX58/IFIH1 has 45 of the 66, and PIP3 has +the SARS-CoV-2 complex above. + +specs/032 held drugs inert: held-out **+30 / 0 broken**. It failed its gate +only on RAF's experimental cases, where the loop in §6 dominates. + +## 6. Counted: M4 (inputs supplied only by the reaction's own downstream) + +| where an input edge's source comes from | edges | share | +|---|---|---| +| a root | 38,599 | 43.7% | +| upstream | 41,047 | 46.4% | +| **downstream only** (every producer is reached from the reaction) | **7,989** | **9.0%** | +| downstream, plus an outside supply | 746 | 0.8% | + +- This happens in 70 pathways: Class I MHC 1,098, DNA Repair 525, DNA Damage + Bypass 422, HOX 421, DSB 407, EGFR 407, Cell Cycle Checkpoints 399, HDR 398. +- **87% (6,948) have no other usable copy** of the entity in the network; 11% + have a root copy and 3% an upstream one. +- The entity exists only inside its own cycle, so nothing outside the cycle + anchors its baseline. That is RAF's F-actin, CNKSR2 and Ca2+ (supplied only by + "Dissociation of RAS:RAF complex"), and the knife-edge loops of specs/013 and + 014. +- MP-BioPath cut such loops by hand, which made the entry entity a root. + `DS_PIN_SCOPE=root_cycle` reproduces that only for the PERTURBED gene. + +## 7. Design candidates this points to (for review, none built) + +- **D1 — anchor cycle-only entities.** An input whose every producer is + downstream of its consumer gets a baseline supply: a root copy (fold 1) + OR-pooled with the recycled copy. This is the curated entity's initial + supply, the pool a cell starts with. + - Mechanism: the any-copy pool (specs/031, 033) with a synthetic root member. + - It removes the all-zero trap without cutting a curated edge. + - Addresses M4 and RAF, and plausibly M3, since recycled catalysts are the + same shape. +- **D2 — exogenous entities inert.** Drugs AND pathogen-derived entities are + held at baseline (specs/032's rule, widened to species). Addresses M5 and + part of M6. +- **D3 — Requirement as a limiter.** A knockout blocks the reaction; a rise + does not raise it, so min(fold, 1). Addresses M1. +- **D4 — set pools (specs/033),** already measured neutral to positive. + +D1 is the largest and the most structural. D2 and D3 are small, principled and +independent. Each is to be pre-registered and measured on its own before any +combination. + +## 8. RAF traced by iteration (2026-09-27) + +HRAS up (80x) on the canonical build. The solve was snapshotted at iteration +budgets 1, 2, 3, 5, 10, 20, 50 and 535. Fixed-point iteration is deterministic +within a process, so each budget continues the same trajectory. + +| iteration | RAS:GTP:RAF complex (input) | activated RAF:scaffold:MAP2K:MAPK | 128-node component below 0.5x | +|---|---|---|---| +| 2 | 1.04x | **0.98x** | 0 | +| 3 | 1.16x | 0.92x | 0 | +| 5 | 1.57x | 0.60x | 4 | +| 10 | 3.76x | 0.04x | 115 | +| 20 | 100x | 6e-5 | 127 | + +**The up-signal turns into a down-signal at the step that forms the activated +scaffold complex.** That step is "MAP2Ks and MAPKs bind to the activated RAF +complex" (R-HSA-5672972). At iteration 3: + +- **The push:** negative regulator R-HSA-5675413, p21 RAS:GTP:activated RAF1 + homo/heterodimer:PEBP1 (RKIP), rises with RAS (1.154x). It contains the + reaction's own input, the self-contained shape of specs/022. +- **The amplifier:** in our network the reaction has **27 AND inputs**, about 20 + of them recycled from the downstream "Dissociation of RAS:RAF complex" (each + 0.994x). Their product is 0.994^20 ≈ 0.89 per pass: a loop gain of about 20 + on any step down, so it collapses instead of settling. + +**Reactome curates this reaction with 4 inputs:** +- p21 RAS:GTP:activated RAF; +- MAPKs (a DefinedSet of 2); +- MAP2K homo/heterodimers (a DefinedSet of 3); +- RAF/MAPK scaffolds (a **CandidateSet of 8**). + +That is at most 2 × 3 × 8 alternatives. The network has **one reaction node +that requires every member of every set at once**, including the leaves of all +8 scaffold candidates (F-actin, TLN1, VWF, fibrin, integrins, CNKSR1/2, KSR2, +IQGAP1, ARRB1/2 …). + +**How common this is:** +- Counted by leaves unique to each member (members sharing a subunit are not + counted twice), **194 reaction copies in 16 pathways** require at least 2 + alternatives of one input set: DDX58/IFIH1 96, Class I MHC 33, Neurexins 13, + RAF 13, RET 12. +- A looser count (any shared leaf) read 3,709. That count was an artefact of + shared subunits and is withdrawn. + +**Hypothesis under test:** the variant cap (`LNG_MAX_VARIANTS`, 512) bundles an +over-cap reaction's alternatives into one node, which is an AND of all of them. +RAF is being regenerated with the cap at 512 and at 10^6 to test it. + +**This also explains the three failed RAF fixes.** Drugs held inert, set +pools, and conserved participants each removed a secondary term. None of them +removed the ~20-way AND of alternatives that turns a small step down into a +collapse. + +## 9. RAF: MP-BioPath's network against ours (2026-09-27) + +MP-BioPath's RAF network (`~/codes_and_results/mp_biopath/V_78_direct/pathways_2/RAF_MAP_kinase_cascade.tsv`, +Reactome v78, loops cut by hand) scores **100%** experimental with DeltaSignal's +own solver (docs/RESULTS §1b). Ours scores 29%. Same solver, so the difference is +the network: + +| | MP-BioPath | ours (v97) | +|---|---|---| +| size | 842 nodes, 1,077 edges | 1,571 nodes, 4,762 edges | +| largest loops | **12 and 6 nodes** | 138 and 128 | +| sets | **one node per set, members OR-ed at the curated-alternative level** ("RAF/MAPK scaffolds" ← its 8 candidates; "RAF activating kinases" ← 7 kinases) | expanded: input sets into reaction copies (and past the variant cap into all-required bundles); catalyst and regulator sets into every member on every copy | +| "MAP2Ks and MAPKs bind to the activated RAF complex" | 4 inputs (the 3 set nodes and RAS:GTP:RAF), inhibitor **free PEBP1** | 28 member-level inputs, inhibitor PEBP1 bound to RAS:GTP:RAF | +| "Phosphorylation of RAF" catalyst | **one** set node | 9 members, all required | +| "RAF phosphorylates MAP2K dimer" | the complex once, plus 2 drug inhibitors | the complex as input AND catalyst (squared) | +| "Dissociation of RAS:RAF complex" outputs | p-MAP2K dimers, MAPKs, 5672712: **no scaffold release** | releases every scaffold leaf (F-actin, CNKSR, talin, …), which recycle into the binding step | +| RAF's loop | the real MAP2K → RAF feedback, through the OR catalyst node | welded through ~20 recycled scaffold leaves | + +**Implication.** Every patch tried here (specs/032, 033, 035, 036, 037) tried to +re-create inside our representation what MP-BioPath's representation has +natively: one node per curated set, OR over its alternatives, and no member-level +fan-out. + +The generator's design is the opposite. Sets "should always decompose to member +species, they should NOT survive as network nodes" (set-expansion fix, July +2026). That choice is what produces copy multiplicity (specs/031), capped +bundles (specs/036) and member fan-out (specs/033). + +**Whether to revisit it is a design decision for Adam, not an experiment to +run unasked.** The scaffold release on dissociation is a curation difference +(v78 against v97), to be checked in Neo4j before it is attributed to the +generator. + +## 10. The MP-BioPath hand edits (2026-09-27) + +**Source:** +- `~/gitroot/mp-biopath-test-neo4j-generated-pathways/PathwayAnalysis/QA/*_changes.tsv` + (87 files; columns source, destination, interaction type, and/or, action, + comment), applied by `bin/add_changes_to_pathway.pl`. +- Unedited: `reactome_exports/`. Edited: `pathways/`. +- Adam: "some of these loops/cycles were in reactome and were intentional, some + others were loops caused by the fact that a single entity in a compartment has + the same id even if it is in multiple places." + +**Totals:** +- **1,800 edits in 71 pathways:** 1,171 deletions, 629 additions. +- Cyclic nodes fell from **1,007 (export) to 413 (edited)**, −59%. +- Comments are rare: "set to member link" ×4, "…with Active AKT", "receptor:ligand + set to member links". + +| action | edge (MP-BioPath node kinds) | n | inside a loop in the export | +|---|---|---|---| +| delete | entity → entity (set-to-member / same-id links) | 607 | 55 | +| delete | entity → reaction (inputs) | 163 | 12 | +| delete | **reaction → entity (outputs)** | 148 | **84** | +| delete | entity → set-decomposition node | 96 | 1 | +| delete | set-decomposition → entity | 76 | 22 | +| delete | **reaction → set-decomposition node** | 63 | **38** | +| delete | negative entity → reaction | 13 | 0 | +| add | positive entity → entity | 182 | – | +| add | **negative entity → entity** | 131 | – | +| add | **negative entity → reaction** | 87 | – | +| add | positive reaction → entity | 78 | – | +| add | positive entity → reaction | 36 | – | +| add | positive, involving an unnamed (synthetic) node | ~116 | – | + +**Three separable kinds of edit:** +1. **Loop cutting.** About 120 of the 208 in-loop deletions remove a reaction's + output into a recycled species or a set's members. RAF: "Dissociation of + RAS:RAF complex" → the MAPK set's member-level AND/OR decomposition was + deleted and replaced by an edge to the set node. +2. **Same-id and set-to-member artefacts** (most of the 607 entity → entity + deletions, which are mostly outside loops). This is the class the + generator's per-position UUIDs are meant to handle. +3. **218 negative edges added** that Reactome does not curate, 131 of them + entity → entity. Our networks have none of them. This is a direct source of + difference that no generator change can recover from Neo4j alone. + +**Next:** +- Map the loop-cutting deletions (class 1) onto our v97 networks: do the same + recycling edges close loops there? That would give a curator-made label set of + "artefact loop edges" to design loop handling against, instead of guessing. +- List the 218 added inhibitions and ask whether each is biology Reactome lacks + or a modelling patch. + +## 11. Our loops are mostly curated reaction cycles, not artefacts (2026-09-27) + +**MP-BioPath's loop-cutting deletions mapped onto our v97 networks** (canonical +`20260926-1221_590301c`): 122 reaction → output edges were deleted INSIDE a loop, +in 25 pathways. + +| where the deleted edge lands in our network | edges | +|---|---| +| present, and **not** in a loop | **85** | +| present, and still closes a loop | 14 (DSB/HDR 5687758 → 68462; Intrinsic Apoptosis BAD / p-S99-BAD, 139904, 141643, 350870; Mitotic G1 69195 → 68377) | +| no such edge | 20 | +| pathway not in our catalog | 3 | + +**So the per-position UUIDs already resolve about 86% (85 of 99) of the loop +artefacts MP-BioPath cut by hand.** + +**What our 8,064 cyclic nodes are made of** (nodes that stop being cyclic when +one edge type is removed): + +| edge type removed | cyclic nodes no longer cyclic | +|---|---| +| output | 8,058 (100%) | +| input | 7,878 (98%) | +| depletion | 1,426 (18%) | +| regulator | 1,033 (13%) | +| catalyst | 487 (6%) | +| dissociation | 0 | +| assembly | 0 | + +**Most cyclic nodes sit in loops made of curated input/output edges alone**: +Reactome's own reaction cycles (interconversions such as RAS GDP ⇄ GTP, +phosphorylation ⇄ dephosphorylation, association ⇄ dissociation). Derived edges +contribute a minority. This supports Adam's position: keep the loops, and +handle them mathematically. + +**The math they need.** A cycle A → B → A driven by opposing enzymes (GEF / GAP, +kinase / phosphatase) has a CONSERVED total whose split is set by the two +activities: B/A ∝ forward / backward. The solver instead multiplies around the +cycle: gain 1 at baseline, the knife-edge of specs/013 and 014, and the RAF +collapse. specs/017's loop pool used the product of entries, which is not +this ratio. **A ratio-based treatment of interconversion cycles is the +untested candidate.** It needs a design pass (identifying interconversion +cycles and their opposing catalysts from the curated reactions) before any +pre-registration. + +## 12. The 219 negative edges MP-BioPath added (2026-09-27) + +Adam: do the work first; the curators can then be asked what was done and +whether it can be replicated automatically in the new networks. + +**Patterns** (by the schema classes of source and target, from Neo4j): + +| pattern | examples | +|---|---| +| **reversal enzyme ⊣ the modified form, or ⊣ the forward reaction** (the largest group) | phosphatases: PPP5C, PTPN6 / PTPN11 ⊣ "Phosphorylation of INFAR1 by TYK2", myosin phosphatase ⊣ p-MRLC; demethylases: KDM1A/B, KDM2A/B, KDM4A, KDM5A-D, KDM6B, JMJD6 ⊣ histone methylation reactions; deubiquitinase OTUD5; **RAS GAPs, SPRED:NF1 ⊣ p21 RAS:GTP** | +| ubiquitin ligase ⊣ its target | ITCH, RNF125:E2, CUL1:SKP1:SKP2:CKS1B | +| sequestering trap or decoy ⊣ ligand | FST, FSTL3 ⊣ Activin; "Ligand Trap" ⊣ BMP2; IL13RA2 ⊣ "IL13 binds IL13RA:TYK2"; I-SMAD ⊣ BMP:BMPR | +| transcriptional repressor ⊣ target-gene expression | RB1 ⊣ E2F1-driven TK1, CCNA1, DHFR, TYMS expression | + +**Which of them our v97 network already has** (canonical build; the source and +target entity with a negative edge between them): + +| status | edges | +|---|---| +| same negative edge exists (a `depletion` edge) | 7 | +| the source is a negative source in ours, but with a different target | 9 | +| **absent** | **187** | +| pathway not in our catalog | 16 | + +**Reading.** The dominant pattern is the **backward enzyme of an +interconversion cycle** (GAP, phosphatase, demethylase, deubiquitinase) used as +a brake on the forward-modified species. It is the hand-made counterpart of §11: +in a cycle A ⇄ B the split B/A follows forward / backward activity, and +MP-BioPath encoded "backward" as an inhibition of B. + +- Our generator derives a small part of it: depletion edges for phosphatases + (Pi as output) and Ub ligases (Ub as input). +- It does not derive it for demethylases, GAPs (GTP hydrolysis), deubiquitinases + or transcriptional repressors. +- The traps and decoys are ordinary sequestration. + +**Automatic equivalent to test:** for each curated reversal reaction B → A with +catalyst E_b, make E_b act against B, or better, give the A ⇄ B cycle its +ratio-based treatment (§11). +- This is derivable from Reactome: pairs of reactions whose input/output + entities are the modified/unmodified forms of one protein. +- It is the pattern to confirm with the curators. diff --git a/specs/035-conserved-participants/research.md b/specs/035-conserved-participants/research.md new file mode 100644 index 0000000..ed53299 --- /dev/null +++ b/specs/035-conserved-participants/research.md @@ -0,0 +1,148 @@ +# specs/035 — Conserved loop participants held at baseline (design D1) + +**Status:** pre-registered 2026-09-26, before any code or arm. + +## Why (specs/034 §6) + +- **7,989 reaction-input edges** (9% of all) in 70 pathways are supplied ONLY by + the reaction's own downstream: every producer of the input node is reached + from the reaction. +- 87% have no other usable copy of the entity, so nothing outside the loop + anchors its baseline. A loop whose AND inputs are all loop-carried has a gain + of 1 at baseline (the knife-edge, specs/013 and 014), and a push can drain it + to the all-zero state. +- **RAF/MAP kinase is traced** (specs/032, 033): the 128-node component still + collapses to 0 under every RAS/RAF perturbation after pooling. Experimental + is 17 of 49, where MP-BioPath gets 46. + +**The nearest failed attempt, and the difference:** +- `DS_SCC_BREAK_CATALYST` froze every recycling catalyst at its entry value: + TP53 +104, held-out −65, not adopted. +- It froze catalysts that also carry real signal. +- This rule is narrower. It touches only entities **released unchanged** by the + loop: every producer of the node consumes a complex that CONTAINS the node's + entity (or is a `dissociation` edge from one). That is a conserved pool, + bound and released, like a cofactor. +- Entities transformed inside the loop (RAS:GTP, the RAF complexes) are not + touched. + +**Split of the 7,989:** + +| class | edges | where | +|---|---|---| +| released unchanged (conserved) | 1,722 (22%) | 36 pathways: Class I MHC 499, DNA Damage Bypass 150, PTK6 150, FGFR2 140, ERBB2 119, ERBB4 112, HDR 81, RAF 67 | +| transformed in the loop | 6,267 (78%) | not touched | + +RAF's conserved set: +- the RAS GAPs: NF1, RASA1–4, RASAL1–3, SYNGAP1, DAB2IP; +- the scaffold partners: F-actin, CNKSR1/2, IQGAP1, KSR2, TLN1, ARRB1/2, + RAP1A/B, APBB1IP; +- the ions. + +## The rule + +`DS_CONSERVED_MODE` = `off` (default, byte-identical) or `inert`. + +A node is **conserved** if all four hold: +1. it is an input of a reaction in a non-trivial strongly connected component; +2. every non-depletion producer of it is reached from that reaction; +3. every producer is a release: a reaction whose input or catalyst CONTAINS the + node's entity (by the bundle's containment table, or the same stId), or a + `dissociation` edge from such a complex; +4. it has at least one producer. + +Under `inert`, conserved nodes are pinned at baseline exactly as cofactors are: +participants at fold 1.0 that never carry a perturbation. An explicit pinning +observation (confidence above the gate) wins, so a perturbed GAP such as NF1 +is still perturbed. The solve reports `conserved_held`, and the benchmark +refuses an `inert` arm that held nothing. + +## Pre-registration + +**Arms:** +- `consctrl`: code defaults; +- `consinert`: `DS_CONSERVED_MODE=inert`; +- both on the canonical build `20260926-1221_590301c`, through + `scripts/run_arm.sh`. + +**Adopt only if all hold:** +- curator held-out net > +15 with p < 0.05; +- experimental net ≥ −15; +- no single pathway loses more than 10 held-out cases; +- the gain is not a single-pathway effect: at least 3 pathways moved, or the + result is reported as such. + +**Predictions:** +- RAF experimental rises (at least +10 of 49); +- RAF convergence improves; +- Class I MHC, DNA Damage Bypass, PTK6 and the FGFR/ERBB pathways move; +- pathways with no conserved nodes are byte-identical. + +**Also reported:** +- tuning curator; +- per-pathway net on both axes; +- concentration; +- convergence counts; +- `conserved_held`. + +## Result — NOT ADOPTED (2026-09-26) + +**Setup:** canonical build `20260926-1221_590301c`; arms `consctrl` and +`consinert` at solver 1fcc1bc. The inert arm held 12,070 curator and 1,256 +experimental node-solves. Convergence: curator 1,652 → 1,658, experimental +217 → 226 of 244. + +| split | net | fixed / broke | p | pathways | +|---|---|---|---|---| +| curator held-out | **+42** | 110 / 68 | 0.002 | 13 (+9 / −3) | +| curator tuning | −265 | 100 / 365 | 1.7e-36 | PIP3 −308 | +| **experimental** | **−122** | 3 / 125 | 2.1e-33 | 3 (all negative) | + +- Best: HDR +54, ERBB2 +26, ERBB4 +18, DNA Damage Bypass +17. +- Worst: PIP3 −308, Mitotic G2/M −42, **TGF-β −30**, **Class I MHC −18**. + +**Gates:** +- held-out passes (+42, p = 0.002); +- **experimental fails** (−122); +- **no-pathway-below −10 fails** (TGF-β −30 and Class I MHC −18, both held-out). + +**Not adopted.** `DS_CONSERVED_MODE` stays in the code as a record, default +`off`. + +**Why PIP3 collapsed:** PIP3 itself (R-ALL-179838) is one of the two nodes held +in that pathway. +- That PIP3 node is produced only by reactions that take PIP3 in and give it + back (it binds its effectors and is released), so it meets "released + unchanged". +- But its level IS the pathway's signal, a second messenger. +- **"Released unchanged" does not separate a conserved scaffold (F-actin, a RAS + GAP) from a signalling molecule that is bound and released.** +- It is also copy multiplicity again: the loop-only PIP3 copy is not the one + PI3K produces. + +**What survives:** +- The held-out gain is real: +42, 13 pathways, 66 perturbations. +- It sits in the pathways the inventory predicted: HDR, ERBB2/4 and DNA Damage + Bypass, the loop-only-input pathways of specs/034 §6. +- A narrower rule (exclude small molecules and second messengers, or require a + protein or complex entity) would be a NEW pre-registration. It is not a tuning + of this one. Not run. + +**Scored per pathway (experimental, ctrl → inert):** + +| pathway | ctrl | inert | +|---|---|---| +| **RAF** | 14 / 49 | **14 / 49** | +| PIP3 | 168 / 200 | **49 / 200** | +| Mitotic G1 | 35 / 88 | 35 / 88 | +| ERBB2 | 42 / 49 | 41 / 49 | + +- **The whole experimental loss is PIP3. RAF, the motivating case, did not + move.** Holding its 29 conserved entities at baseline does not stop the + collapse either. +- Three attempts have now failed to fix RAF: inert drugs (specs/032), set + pools (specs/033) and conserved participants (here). Each targets a + structural suspect found by reading the graph. +- **The next step is to trace the RAF solve's iterations directly** (which + node first leaves baseline, in which direction, and why the loop does not + settle) before proposing another rule. diff --git a/specs/036-capped-bundles-as-pools/research.md b/specs/036-capped-bundles-as-pools/research.md new file mode 100644 index 0000000..b355138 --- /dev/null +++ b/specs/036-capped-bundles-as-pools/research.md @@ -0,0 +1,120 @@ +# specs/036 — Over-cap set bundles as one input per curated component + +**Status:** pre-registered 2026-09-27, before any code or arm. + +## The defect (traced in specs/034 §8) + +The variant cap (`LNG_MAX_VARIANTS`, 512; issue #40) exists because expanding +every set combination is intractable. Uncapped, RAF alone asks for a +52.6 GiB (84,000 × 84,000) array. + +**The problem is what the cap does:** when a reaction's input or output +combinations exceed it, `get_broken_apart_ids` merges every alternative of every +set into ONE member list, a single variant that **requires all of them**. + +**RAF, "MAPs and MAPKs bind to the activated RAF complex" (R-HSA-5672972):** +- Curated with 4 inputs: RAS:GTP:activated RAF; MAPKs (a set of 2); MAP2K dimers + (a set of 3); RAF/MAPK scaffolds (a CandidateSet of 8). +- The generator log: "variant cap hit for R-HSA-5672972: 12000 combinations > + 512; bundling 4 components into one node". +- The network node requires 28 inputs at once. +- The other members of the RAF loop are capped the same way: R-HSA-5672720, + 5672723, 5672724 and 5672980. +- About 20 of those inputs are recycled from the downstream dissociation, so a + 0.6% step down becomes 0.994^20 ≈ 0.89 per pass. The loop collapses to 0 under + every RAS/RAF perturbation (experimental 14 of 49, against MP-BioPath's 46). + +**Scale:** the canonical build's `generate.log` records **96 cap hits in 21 +pathways**: RUNX1 24, HOX 19, WNT 6, CD28 6, RAF 5, Rho GTPases 5, and others. + +## The rule + +**`LNG_CAP_POOLS=1`:** at emission, for a reaction or complex whose inputs come +from a capped bundle, the bundled leaves are grouped by the CURATED input +component they belong to (from Neo4j: the reaction's inputs and their set +membership). +- A component that is a set becomes a **set pool** (specs/033): its leaves feed + the pool by `set_member` edges, and the pool feeds the reaction as ONE AND + input. +- A component that is not a set is wired as before. +- The reaction then has one input per curated component, as Reactome curates it. +- The matching layer (the cap itself, and the content hashes that line inputs up + with outputs) is **unchanged**. + +The solver reads pools with `DS_SET_POOL_MODE` (default `product`). For a pool +of alternatives, `product` would re-create the all-required bundle, so this +change is **measured with `extreme` and `mean`** as well. Which one to use is +chosen on the tuning split, as specs/033 did. + +## Pre-registration + +**Arms,** all through `scripts/run_arm.sh`: +- `base`: the canonical build `20260926-1221_590301c`, code defaults; +- `cap_product`, `cap_extreme`, `cap_mean`: one build with `LNG_CAP_POOLS=1` and + `LNG_SET_POOL=1`, with each pool mode. + +**Choosing the mode:** by tuning curator macro-F1 plus experimental macro-F1; +ties go to product < extreme < mean. + +**Adopt** the chosen arm only if all hold: +- curator held-out net ≥ −15 and not significantly negative; +- experimental net ≥ 0; +- no single pathway loses more than 10 held-out cases. + +It is a faithfulness fix (one input per curated component), so it does not need +a held-out gain to be adopted. + +**Predictions:** +- RAF experimental rises well above 17 of 49 under `extreme` or `mean`, and the + RAF solve converges; +- `cap_product` is close to base; +- the 21 capped pathways are the only ones that move, together with the + specs/033 set-pool pathways, which are measured separately there. + +**Also reported:** +- per-pathway net on both axes; +- concentration; +- convergence counts; +- the RAF iteration trace re-run on the new build; +- cyclic node count excluding pools. + +## Sanity trace on RAF alone (2026-09-27, before any arm) + +**Setup:** +- RAF (R-HSA-5673001) regenerated at LNG feat/set-pool + WIP `LNG_CAP_POOLS=1`, + in the scratchpad (a diagnostic, not a catalog). +- "MAP2Ks and MAPKs bind to the activated RAF complex": **28 inputs → 4**. The + 4 are RAS:GTP:activated RAF plus pools of MAPKs (2), MAP2K dimers (3) and + RAF/MAPK scaffolds (22 leaves of the 8 candidates). +- **Limitation:** that pool is over leaves, not candidates, so larger candidates + weigh more. + +HRAS pinned 80x or 0. Readout: p-T,Y MAPKs in the nucleus (R-HSA-5674340). + +| configuration | HRAS up | HRAS KO | +|---|---|---| +| canonical | 0 (collapse) | – | +| cap pools, any pool mode | 0 (collapse), not converged | 100x (inverted) | +| + drug-derived nodes held at baseline | 0 (collapse) | 100x (inverted) | +| + `DS_DEDUP_ACTIVATORS=1`, `mean` pools | **0.83x**, no collapse (0 zeros), not converged | **0.016x** (correct direction) | +| + `DS_DEDUP_ACTIVATORS=1`, `extreme` pools | 0 (collapse) | 0 | + +**Amplifiers found, in the order they appeared:** +1. Capped bundles: all alternatives required at once. +2. Drug-bound complexes inhibiting "RAF phosphorylates MAP2K dimer" and "MAP2Ks + phosphorylate MAPKs" (R-HSA-9653109 and others rise with RAS). +3. **Catalyst = substrate squaring.** + - "RAF phosphorylates MAP2K dimer" reads the activated RAF:scaffold complex + as both input and catalyst. + - "MAP2Ks phosphorylate MAPKs" does the same with its complex. + - Each squares its own fold (specs/012; de-duplicating was held-out −15 + catalog-wide there). +4. **Open:** with 1–3 addressed and `mean` pools, up-regulation still reads + 0.83x and the component does not converge. Not yet traced. + +**Consequence for the pre-registration:** the RAF prediction ("rises well above +17/49") cannot be met by cap pools alone. Before a combined arm (cap pools + +inert drugs + dedup) is pre-registered, amplifier 4 is traced. Each component +was individually neutral or negative catalog-wide, so the combination is +measured as a combination, and the RAF effect is not generalised from this +trace. diff --git a/specs/037-combined-raf-fixes/research.md b/specs/037-combined-raf-fixes/research.md new file mode 100644 index 0000000..8e3a155 --- /dev/null +++ b/specs/037-combined-raf-fixes/research.md @@ -0,0 +1,154 @@ +# specs/037 — Combined structural fixes (from the RAF trace), measured catalog-wide + +**Status:** pre-registered 2026-09-27, before the build or any arm. + +## Why + +Tracing RAF/MAP kinase iteration by iteration (specs/034 §8, specs/036) found +compounding amplifiers. Each one alone was neutral or harmful catalog-wide, or +did not move RAF: + +| # | amplifier | fix | alone | +|---|---|---|---| +| 1 | variant-capped set inputs, with every alternative required at once | `LNG_CAP_POOLS=1` (specs/036) | RAF still collapses | +| 2 | set-valued catalysts/regulators as separate AND terms | `LNG_SET_POOL=1` (specs/033) | held-out +6, experimental +3 | +| 3 | drug-bound complexes inhibiting the kinase steps | `DS_DRUG_MODE=inert` (specs/032) | held-out +30, experimental −5 | +| 4 | catalyst = substrate squaring | `DS_DEDUP_ACTIVATORS=1` (specs/012) | held-out −15 | + +On RAF alone, 1 + 3 + 4 with `mean` pools turned HRAS KO from "no signal" into +0.016x (correct). HRAS up read 0.83x (still wrong, but no longer a collapse). + +## Arms + +All through `scripts/run_arm.sh` at one solver commit: + +| arm | build | settings | +|---|---|---| +| `base` | canonical `20260926-1221_590301c` | code defaults | +| **`combined`** | a new build from LNG `feat/combined-037` (`LNG_SET_POOL` default on, `LNG_CAP_POOLS=1`, ships `drugs.csv`) | `DS_SET_POOL_MODE=mean`, `DS_DRUG_MODE=inert`, `DS_DEDUP_ACTIVATORS=1` | + +Attribution arms, on the same new build, are **reported, not gated**: + +| arm | settings | +|---|---| +| `pools_mean` | `DS_SET_POOL_MODE=mean` | +| `pools_drugs` | + `DS_DRUG_MODE=inert` | +| `pools_dedup` | + `DS_DEDUP_ACTIVATORS=1` | + +## Pre-registration (`combined` against `base`) + +**Adopt only if all hold:** +- experimental net ≥ +15 with p < 0.05 (the axis this targets: we trail + MP-BioPath's published result by 8.6pp); +- curator held-out net ≥ −15 and not significantly negative; +- no held-out pathway loses more than 10 cases. + +**Predictions:** +- RAF experimental rises (at least +10 of 49); +- PIP3 is not harmed (the specs/035 failure mode); +- the attribution arms show which component carries RAF. + +`mean` is fixed in advance from the RAF trace. It is not chosen on the +evaluation set; `extreme` still collapsed there. + +**Also reported:** +- tuning curator; +- per-pathway net; +- concentration; +- convergence counts; +- `drugs_held`; +- cyclic node count excluding pools; +- the RAF iteration trace on the new build; +- the paired gap against MP-BioPath's supplementary table. + +## Revision after review (2026-09-27, before any arm on the combined build) + +An adversarial review of both integration branches found a blocking generator +bug and weaknesses in this pre-registration. The first combined build was +stopped part-way, and `20260927-0821_d98d2f8_combined037` is incomplete and +unused. + +**Generator (LNG bddfe78):** +- Cap pools fired in reactions the cap never touched: 410 pools in PIP3, which + has 0 cap hits. The miR-93 RISC input's four subunits were averaged, so a + miR-93 knockout read 0.75x instead of 0. +- Now: + - pooling happens only in capped reactions; + - leaves are grouped by alternative, a complex alternative is an AND unit, and + a pool needs at least 2 alternatives; + - shared leaves stay direct; + - the curated stoichiometry is kept; + - depletion onto a pool is retargeted to its members. +- Verified by regeneration: RAF's binding reaction has its 4 curated inputs, + and PIP3 has no input pools (its 11 pools are the specs/033 catalyst and + regulator pools). + +**Pre-registration, revised:** +1. **Control is a fresh regeneration**, not the canonical build. It is built + from the same LNG commit with `LNG_SET_POOL=0` and `LNG_CAP_POOLS=0`. `base` + and `combined` were otherwise different uuid draws, and the gates sit at the + regeneration noise floor. +2. **The experimental gate excludes RAF.** `mean` and dedup were chosen from + the RAF HRAS trace, and RAF is an experimental pathway, so gating on it would + select on the evaluation set. The gate is experimental net ≥ +15, p < 0.05, + **excluding RAF/MAP kinase**. RAF is reported separately as the traced case. +3. **Concentration gate:** the curator held-out gain or loss, and the + experimental gain, must span at least 2 pathways and at least 5 distinct + perturbations. Otherwise McNemar is not valid and the result is reported as + such. +4. **The tuning pathways are gated too:** no pathway, held-out or tuning, may + lose more than 10 cases. PIP3 (specs/035's failure) is covered. +5. Corrections to the table above: + - the specs/033 +6 / +3 was measured under `product`, not `mean`, so it does + not carry over; + - this spec replaces specs/036's plan (its `cap_*` arms and its choice of + mode on the tuning split) with `mean` fixed in advance. +6. **Known limitations, stated:** + - cap pools are minted per reaction copy, not shared; + - pool `member_leaves` in `nodes.csv` lists every leaf of the set; + - de-duplication does not see through a pool. + +## Result — NOT ADOPTED (2026-09-27) + +**Setup:** +- Control: `20260927-0833_bddfe78_ctrl037`, a fresh regeneration with pools off. + Its structure matches the canonical build exactly (72,025 nodes, 225,147 + edges). +- Pooled build: `20260927-0845_bddfe78_comb037` (72,348 nodes, 214,234 edges). +- All arms at solver 6f261e2. + +| arm | curator held-out | curator tuning | experimental (all) | experimental excl. RAF | RAF experimental | +|---|---|---|---|---|---| +| control | — | — | — | 555 / 796 | 14 / 49 | +| **combined** (pools + drugs + dedup, `mean`) | **−162** (53/215, p = 2.6e-24) | −106 | −3 | **−12** | 23 / 49 | +| pools only (`mean`) | −35 (p = 6.9e-7) | −29 | +6 | −2 | 22 / 49 | +| pools + drugs | −4 (p = 0.74) | −36 | +7 | −2 | 23 / 49 | +| pools + dedup | −193 (p = 1.7e-40) | −104 | −6 | −12 | 20 / 49 | + +**Gates for `combined`:** +- experimental excluding RAF ≥ +15: it is −12, **FAIL**; +- held-out: −162, **FAIL**; +- no pathway below −10: **FAIL** (IFN α/β −112, TP53 −80, IFN-γ −34, PIP3 −28). + +**Not adopted.** + +**What the attribution says:** +- **De-duplication is the damage.** pools + dedup reproduces the combined + losses: IFN α/β −112, TP53 −80, IFN-γ −34, PIP3 −28. That confirms and enlarges + specs/012's negative (−15 there). Squaring catalyst-and-substrate is + load-bearing on these networks, even though it amplifies RAF's loop. +- **Pools with `mean`** cost held-out −35 (DSB −13, NOTCH1 −7, TP53 −37 tuning). + They move RAF 14 → 22. +- **Pools + drugs** are neutral on held-out (−4) and experimental (+7; −2 + excluding RAF), with the specs/032 held-out gains (MET +11, ROCKs +8, VEGF +6, + KIT +6) offset by the pool losses. +- **RAF gains in every pooled arm (14 → 20–23 of 49).** The cap pools partly + fix the traced case, but nothing generalises beyond RAF on the experimental + axis. + +**Record:** +- The generator's cap pools stay behind `LNG_CAP_POOLS` (default off). They are + correct, and they fix the structural defect (RAF's binding reaction has its 4 + curated inputs), but they are not adoptable on accuracy with `mean`. +- A `product`-mode cap-pool arm was not run. With cap pools under `product`, the + pool re-creates an AND over alternatives, which is the defect itself. diff --git a/specs/038-pool-max/research.md b/specs/038-pool-max/research.md new file mode 100644 index 0000000..203f78e --- /dev/null +++ b/specs/038-pool-max/research.md @@ -0,0 +1,115 @@ +# specs/038 — Set-valued regulators and catalysts as one OR node, read as `max` + +**Status:** pre-registered 2026-09-27; result recorded below (same day). + +## Why (Adam, 2026-09-27) + +"A reaction should be made unique based on its inputs and outputs, not the +regulator. The regulator should be merged into an OR node, and then have that +regulate the reaction, because any of them fulfil the regulation of the +reaction. We just have to figure out the math." + +The OR node is the set pool of specs/033 (`LNG_SET_POOL`, draft LNG #99). Its +math was measured under `product`, `extreme`, `geomean` and `mean`, but not under +the literal OR: **`max`**, the strongest member fulfils the role. + +- A knocked-out member is covered by its unperturbed peers. +- A raised member lifts the pool. + +The ground truth leans the other way on knockouts (specs/033: a single set +member is expected to change the readout in 68% of curator cases). So this is a +real question, not a formality. + +**Scope:** pools for catalyst and regulator sets only (`LNG_SET_POOL=1`, +`LNG_CAP_POOLS=0`). Inputs keep their virtual-reaction expansion, as Adam +specifies. + +## Arms + +All at one solver commit, through `scripts/run_arm.sh`: + +| arm | build | settings | +|---|---|---| +| `ctrl038` | `20260927-0833_bddfe78_ctrl037` (the fresh control, pools off) | code defaults, re-run at this commit | +| `pmax` | a new build at LNG `bddfe78`, `LNG_SET_POOL=1`, `LNG_CAP_POOLS=0` | `DS_SET_POOL_MODE=max` | +| `pprod` | the same build | `DS_SET_POOL_MODE=product` (the specs/033 reference) | + +## Adopt `max` only if all hold (against `ctrl037`) + +- curator held-out net ≥ −15 and not significantly negative; +- experimental net ≥ 0; +- no pathway, tuning or held-out, loses more than 10 cases; +- any gain spans at least 2 pathways and 5 perturbations (or is reported as + concentrated). + +**Predictions:** +- `max` loses single-member knockout cases (those were 68% "expected change"); +- `max` gains on up-regulation; +- the net depends on which dominates; +- `pprod` reproduces specs/033's +6 / +10 / +3 within the regeneration noise + floor. + +**Also reported:** +- the set-member-case table of specs/033 for each arm (KO and UP separately); +- per-pathway net; +- convergence counts. + +## Result (2026-09-27) + +**Setup:** +- Control `ctrl038`: build `20260927-0833_bddfe78_ctrl037`, pools off, code + defaults. +- Pooled build `20260927-1526_bddfe78_sp038`: `LNG_SET_POOL=1`, `LNG_CAP_POOLS=0`. +- All arms at solver 5a2a927. + +| arm | curator held-out | tuning | experimental | worst pathways | +|---|---|---|---|---| +| **`max`** | **−120** (27 / 147, p = 3.5e-21) | −208 | **−57** (16 / 73, p = 7.2e-10) | PIP3 −140, TP53 −57, IFN-γ −17, IL-4/13 −14 | +| **`product`** | −2 (2 / 4, p = 0.69) | +45 (TP53 +42) | **+8** (8 / 0, p = 0.0078) | Fanconi −3 | + +**`max`: NOT ADOPTED.** It fails every gate. +- "Any one member fulfils the role" covers a knocked-out member with its + unperturbed peers, and the ground truth expects that knockout to matter + (specs/033: 68% of single-member curator cases expect a change). +- PIP3's catalyst sets (46 pool → catalyst edges) carry most of the loss. + +**`product`: passes every gate.** +- Held-out −2, which is within the −15 floor and not significant. +- Experimental +8 (p = 0.0078). +- No pathway below −10. + +The experimental gain spans 2 pathways but only **4 perturbations**, so it is +reported as **concentrated** under the pre-registered concentration rule. +TP53 +42 on tuning is consistent with the uuid-draw effect seen before +(specs/026, 030) and is not attributed to pooling. + +**Conclusion.** One OR node per set-valued regulator or catalyst (Adam's +design) is **adoptable on faithfulness under `product`**, and neutral to +slightly positive, as specs/033 found. The pool's math is the product of member +folds: every member's change counts, which is what the ground truth rewards. +The literal "any one fulfils it" reading (`max`) is refuted on both axes. + +**Scoring the pre-registered predictions:** +- `max` loses single-member knockout cases: **yes**, and on both axes. +- `pprod` reproduces specs/033's +6 / +10 / +3: **partly**. Here it is −2 / +45 / +8. + - Held-out and experimental are within the regeneration noise floor + (held-out ±15, experimental ±15). + - Tuning differs by 35. That is the TP53 uuid draw (+42 here), not pooling. + +**Review of the consolidated PRs** (added before merge): + +*Self-inhibitor flags.* A pool uuid maps to the set's stId, so +`containment.csv` gains the set. This could have changed which inhibitors +specs/022 flags. Measured with `self_inhibitor_setup` over every pathway: + +| build | flagged self-inhibitor slots | +|---|---| +| control | 1,215 | +| pooled | 1,192 | + +Four pathways differ: RAF 9 → 0, DDX58/IFIH1 67 → 61, Cell junction +organization 29 → 24, and R-HSA-450294 3 → 0. The effect is small, runs toward +fewer flags, and is included in the `pprod` numbers above. + +*Edge case.* A pool that also receives a depletion edge falls out of the pool +branch and back to the OR mean. None was found in the catalog. diff --git a/src/DeltaSignal.jl b/src/DeltaSignal.jl index 7be513f..151158b 100644 --- a/src/DeltaSignal.jl +++ b/src/DeltaSignal.jl @@ -20,8 +20,10 @@ include("core/sensitivity.jl") include("core/aggregators.jl") include("core/hill_functions.jl") include("core/cofactors.jl") +include("core/drugs.jl") include("core/silo_bridges.jl") include("core/reaction_model.jl") +include("core/conserved.jl") # Solvers include("solvers/steady_state.jl") @@ -45,6 +47,8 @@ export LogicNetworkEdge, ReactionParams, SolverResult # Parsing functions export parse_logic_network, parse_uuid_mapping, parse_set_mappings, parse_complete_network, create_reaction_network export parse_containment, network_containment_json, containment_from_json +export parse_drug_list, drug_stids_json, drug_stids_from_json, drug_mode, drug_uuids +export conserved_mode, conserved_uuids export convert_to_reaction_network, aggregate_to_pathway_view, validate_network_pathway_mapping export compute_influence_scores, SteadyStateParams diff --git a/src/api/server.jl b/src/api/server.jl index f5e9cd2..a267a40 100644 --- a/src/api/server.jl +++ b/src/api/server.jl @@ -320,7 +320,13 @@ function reaction_network_from_json(data)::DeltaSignal.ReactionNetwork catch throw(ArgumentError("`containment` must be an object mapping a stable id to a list of stable ids.")) end - return DeltaSignal.ReactionNetwork(nodes_dict, edges, set_mappings, cofactor_stids, containment) + drug_stids = try + DeltaSignal.drug_stids_from_json(get(data, :drug_stids, nothing)) + catch + throw(ArgumentError("`drug_stids` must be a list of stable ids or null.")) + end + return DeltaSignal.ReactionNetwork(nodes_dict, edges, set_mappings, cofactor_stids, containment, + drug_stids) end """ @@ -699,6 +705,9 @@ function parse_handler(req) # The bundle's cofactor list, so a parse -> POSTed solve keeps it # (reaction_network_from_json reads it back; it was never sent). "cofactor_stids" => sort(collect(network.cofactor_stids)), + # specs/032: the bundle's drug list (null when it has none), so a + # parse -> POSTed solve runs the same model as a solve by network_id. + "drug_stids" => DeltaSignal.drug_stids_json(network), ) return HTTP.Response(200, JSON_HEADERS, JSON3.write(result)) @@ -836,6 +845,12 @@ function solve_handler(req) # specs/022 (additive): inhibitor slots damped as self-contained. "self_inhibitors" => get(solver_result.diagnostics, "self_inhibitors", 0), "self_inhibitor_rule" => get(solver_result.diagnostics, "self_inhibitor_rule", "unknown"), + # specs/032 (additive): "propagate", "inert", or "inert: no drug table". + "drug_rule" => get(solver_result.diagnostics, "drug_rule", "unknown"), + "drugs_held" => get(solver_result.diagnostics, "drugs_held", 0), + # specs/035 (additive): "off" or "inert", and how many nodes were held. + "conserved_rule" => get(solver_result.diagnostics, "conserved_rule", "unknown"), + "conserved_held" => get(solver_result.diagnostics, "conserved_held", 0), "scc" => Dict( "method" => get(solver_result.diagnostics, "scc_method", "unknown"), "pooled" => get(solver_result.diagnostics, "scc_pooled", 0), diff --git a/src/core/conserved.jl b/src/core/conserved.jl new file mode 100644 index 0000000..c26bf4f --- /dev/null +++ b/src/core/conserved.jl @@ -0,0 +1,81 @@ +""" +Conserved loop participants (specs/035, design D1 of specs/034). + +About 9% of reaction-input edges are supplied ONLY by the reaction's own +downstream: every producer of the input node lies in the same cycle. Nothing +outside the loop anchors such a node, so a loop whose AND inputs are all +loop-carried sits on the gain-1 knife-edge (specs/013, 014), and a push can +drain it to the all-zero state. RAF/MAP kinase is the traced case: F-actin, +CNKSR2 and Ca2+ come back only from "Dissociation of RAS:RAF complex", and the +component reads 0 under every RAS/RAF perturbation. + +Most of those inputs (78%) are TRANSFORMED inside the loop (RAS:GTP, the RAF +complexes) and carry signal; they are not touched. This rule applies only to +nodes RELEASED UNCHANGED by the loop: every producer consumes a complex that +contains the node's entity. That is a conserved pool, bound and released, like +a cofactor, whose level is set by total abundance rather than by loop flux. + + DS_CONSERVED_MODE=off (default) previous behaviour. + DS_CONSERVED_MODE=inert conserved nodes are pinned at baseline, as cofactors + are. An explicit pinning observation still wins. + +Different from DS_SCC_BREAK_CATALYST, which froze every recycling catalyst +(including signal-carrying ones) and lost held-out -65. +""" + +const DS_VALID_CONSERVED_MODES = Set(["off", "inert"]) + +function conserved_mode()::String + value = get(ENV, "DS_CONSERVED_MODE", "off") + value in DS_VALID_CONSERVED_MODES || throw(ArgumentError( + "DS_CONSERVED_MODE must be one of $(join(sort(collect(DS_VALID_CONSERVED_MODES)), ", ")); got $(repr(value))")) + return value +end + +"""uuids of conserved loop participants in `network` (see module docstring).""" +function conserved_uuids(network::ReactionNetwork)::Set{String} + ids = Set{String}(keys(network.nodes)) + for e in network.edges + push!(ids, e.parent_uuid); push!(ids, e.child_uuid) + end + order = sort!(collect(ids)) # label-independent indexing + idx = Dict(u => i for (i, u) in enumerate(order)) + adj = [Set{Int}() for _ in order] + incoming = [LogicNetworkEdge[] for _ in order] + for e in network.edges + push!(adj[idx[e.parent_uuid]], idx[e.child_uuid]) + push!(incoming[idx[e.child_uuid]], e) + end + comp, _ = tarjan_scc_components(adj) + size = Dict{Int, Int}() + for c in comp + size[c] = get(size, c, 0) + 1 + end + stid(u) = (n = get(network.nodes, u, nothing); n === nothing ? nothing : n.reactome_id) + contains(container, s) = container !== nothing && s !== nothing && + (container == s || s in get(network.containment, container, Set{String}())) + + out = Set{String}() + for e in network.edges + e.edge_type == "input" || continue + s, t = idx[e.parent_uuid], idx[e.child_uuid] + (comp[s] == comp[t] && size[comp[s]] > 1) || continue + su = order[s] + su in out && continue + ss = stid(su) + ss === nothing && continue + prods = [p for p in incoming[s] if p.edge_type != "depletion"] + isempty(prods) && continue + all(p -> comp[idx[p.parent_uuid]] == comp[s], prods) || continue + released = all(prods) do p + if p.edge_type == "dissociation" + return contains(stid(p.parent_uuid), ss) + end + # the producing reaction consumes (or is catalysed by) a complex holding ss + any(q -> q.edge_type in ("input", "catalyst") && contains(stid(q.parent_uuid), ss), + incoming[idx[p.parent_uuid]]) + end + released && push!(out, su) + end + return out +end diff --git a/src/core/drugs.jl b/src/core/drugs.jl new file mode 100644 index 0000000..0138560 --- /dev/null +++ b/src/core/drugs.jl @@ -0,0 +1,49 @@ +""" +Drug-derived entities (specs/032). + +Reactome curates drug actions inside signalling pathways: MAP2K and MAPK +inhibitors bind the RAF:scaffold:MAP2K:MAPK complex in the RAF/MAP kinase +cascade, and PARP inhibitors bind PARP in DNA repair. The generator keeps them, +because they are curated, and ships `drugs.csv` naming which nodes are +drug-derived (a Drug; a complex with any drug component; a set whose every +member is one). + +A benchmark case, or any question about an untreated cell, describes a cell +WITHOUT the drug. Left to propagate, a drug is a root held at baseline, so +raising RAS raises the drug-bound complexes, and their inhibitory edges then +shut the kinase cascade off: every RAS/RAF perturbation in RAF/MAP kinase read +DOWN (experimental 28.6% against MP-BioPath's 93.9%, specs/032). + + DS_DRUG_MODE=propagate (default) — previous behaviour. + DS_DRUG_MODE=inert — drug-derived nodes are pinned at baseline, as + cofactors are: participants of fold 1.0 that + never carry a perturbation. Not 0, because under + `divide` an inhibitor at 0 de-represses tenfold, + which would model drug WITHDRAWAL. An explicit + observation still wins. +""" + +const DS_VALID_DRUG_MODES = Set(["propagate", "inert"]) + +function drug_mode()::String + value = get(ENV, "DS_DRUG_MODE", "propagate") + if !(value in DS_VALID_DRUG_MODES) + throw(ArgumentError( + "DS_DRUG_MODE must be one of $(join(sort(collect(DS_VALID_DRUG_MODES)), ", ")); got $(repr(value))")) + end + return value +end + +"""uuids in this network whose Reactome stable id the bundle declares drug-derived.""" +function drug_uuids(network::ReactionNetwork)::Set{String} + out = Set{String}() + stids = network.drug_stids + stids === nothing && return out + for (uuid, node) in network.nodes + rid = node.reactome_id + if rid !== nothing && rid in stids + push!(out, uuid) + end + end + return out +end diff --git a/src/core/reaction_model.jl b/src/core/reaction_model.jl index 4d34f7c..55df6ff 100644 --- a/src/core/reaction_model.jl +++ b/src/core/reaction_model.jl @@ -194,7 +194,9 @@ function create_reaction_from_edges( push!(activator_is_assembly, edge.edge_type in ("assembly", "composition")) # Group composition inputs by the source's base stable id, so node # copies of one entity can be aggregated as alternatives. - g = 0 + # specs/033: a `set_member` edge feeds a set POOL node; -1 marks it + # (composition groups are > 0, ordinary inputs 0). + g = edge.edge_type == "set_member" ? -1 : 0 if edge.edge_type == "composition" node = get(network.nodes, edge.parent_uuid, nothing) rid = node === nothing ? nothing : node.reactome_id @@ -261,7 +263,8 @@ function create_reaction_from_edges( # composition. `max` let the group id (read as "is composition" # under DS_COMPOSITION_MODE=limit*) delete the input role and turn # a 4x input into a <=1 limiter. - d_grp[j] = (d_grp[j] > 0 && activator_group[k] > 0) ? max(d_grp[j], activator_group[k]) : 0 + d_grp[j] = (d_grp[j] > 0 && activator_group[k] > 0) ? max(d_grp[j], activator_group[k]) : + (d_grp[j] == -1 && activator_group[k] == -1) ? -1 : 0 end end activators, activator_is_and = d_act, d_and @@ -951,6 +954,7 @@ struct ReactionEvalConfig depletion_own_product::String self_inhibitor_weight::Float64 # < 0 = off. Default 0.1. specs/022. loop_gain::Float64 # 1.0 = off. specs/028. + set_pool_mode::String # specs/033: how a set pool combines its members. end """ @@ -976,6 +980,8 @@ explicitly so that selecting it stays deliberate): - `DS_COMPOSITION_MODE`: `assembly` is the `else` (composition edges join the assembly-limiting AND cluster); `limit` makes them a pure limiter. - `DS_INHIBITOR_FLOOR_SCOPE`: `none` is the `else` (no floor applied). +- `DS_SET_POOL_MODE` (specs/033): `mean` is the existing OR rule, so it takes + no special branch; the other three combine a set pool's member folds. """ const DS_VALID_MODES = Dict( "DS_INHIBITION_MODE" => Set([ @@ -991,6 +997,7 @@ const DS_VALID_MODES = Dict( "DS_COMPOSITION_MODE" => Set(["assembly", "limit", "limit_novel"]), "DS_DEPLETION_OWN_PRODUCT" => Set(["full", "suppress_only"]), "DS_INHIBITOR_FLOOR_SCOPE" => Set(["loops", "all", "transcription", "none"]), + "DS_SET_POOL_MODE" => Set(["product", "extreme", "geomean", "mean", "max"]), ) const DS_BREAK_ROLES = Set(["catalyst", "assembly", "depletion"]) @@ -1315,6 +1322,7 @@ function resolve_reaction_eval_config()::ReactionEvalConfig # Adam: a loop without inhibitors should still amplify, so this is # meant to be tiny (g close to 1). 1.0 = off, byte-identical. _loop_gain_env(), + _mode_env("DS_SET_POOL_MODE", "product"), ) if cfg.loop_gain < 1.0 && cfg.loop_elasticity < 1.0 throw(ArgumentError("DS_LOOP_GAIN and DS_LOOP_ELASTICITY are two ways to pull a loop " * @@ -1387,6 +1395,56 @@ function _loop_elasticity_env()::Float64 eps end +""" + set_pool_value(vals, bl, mode) + +Value of a set pool node (specs/033) from its members' values `vals`, baseline +`bl`, on the internal 0-1 scale (so the cap is 1.0, i.e. 100x): + +- `product`: the product of member folds, as the AND of the members it replaces; +- `extreme`: the member fold furthest from 1 in log space (a knocked-out member + wins outright; a tie goes to the lower fold, so a knockout beats an equal + rise); +- `geomean`: the geometric mean of member folds (0 if any member is 0); +- `max`: the strongest member (specs/038): "any one of them fulfils the role", + so a knocked-out member is covered by its peers and a raised one lifts the pool. +""" +function set_pool_value(vals::AbstractVector{T}, bl::T, mode::String) where {T<:Real} + bl > zero(T) || return zero(T) + if mode == "product" + f = one(T) + for v in vals + f *= max(v, zero(T)) / bl + f == zero(T) && break + end + return clamp(bl * f, zero(T), one(T)) + elseif mode == "extreme" + best = one(T); dev = T(-1) + for v in vals + fv = max(v, zero(T)) / bl + d = fv <= zero(T) ? T(Inf) : abs(log(fv)) + if d > dev || (d == dev && fv < best) + best, dev = fv, d + end + end + return clamp(bl * best, zero(T), one(T)) + elseif mode == "max" + m = zero(T) + for v in vals + m = max(m, v) + end + return clamp(m, zero(T), one(T)) + elseif mode == "geomean" + s = zero(T) + for v in vals + v <= zero(T) && return zero(T) + s += log(v / bl) + end + return clamp(bl * exp(s / length(vals)), zero(T), one(T)) + end + throw(ArgumentError("set_pool_value: unknown mode $(repr(mode))")) +end + """ Compute reaction output from a flat activity vector. Same math as the dict version (sensitivity transform → geomean activators × Hill-suppression @@ -1528,6 +1586,16 @@ function compute_reaction_output_vec(x::AbstractVector{T}, rxn::IndexedReaction; end end + # specs/033: a set POOL node (every input a `set_member` edge) is one + # curated participant, "any one of these members". Its members arrive as + # OR inputs, so `mean` is the existing OR rule; the other modes combine the + # member folds here and return the pool directly. + if config.set_pool_mode != "mean" && !isempty(rxn.activator_group) && + isempty(rxn.inhibitor_indices) && isempty(rxn.depletion_indices) && + all(==(-1), rxn.activator_group) && length(or_vals) == length(rxn.activator_group) + return set_pool_value(or_vals, bl, config.set_pool_mode) + end + # Inject the assembled-complex limiting value as a single AND input, so it # combines with any catalytic/regulatory inputs via the normal DS_AND_MODE. # Each entity group contributes its strongest copy to the limiting-reactant min. diff --git a/src/io/tsv_parser.jl b/src/io/tsv_parser.jl index c896ba2..3414c3f 100644 --- a/src/io/tsv_parser.jl +++ b/src/io/tsv_parser.jl @@ -42,11 +42,17 @@ struct ReactionNetwork # Empty when the bundle has no such file, or when a network arrives as JSON # (the rule is then inert, which the solve reports). containment::Dict{String, Set{String}} + # Drug-derived stable ids that shipped with this network, from the + # generator's `drugs.csv` (specs/032). `nothing` means the bundle has no such + # file, which is NOT the same as a file listing none: DS_DRUG_MODE=inert + # reports the first as "inert: no drug table" instead of pretending. + drug_stids::Union{Nothing, Set{String}} ReactionNetwork(nodes, edges, set_mappings, cofactor_stids = Set{String}(), - containment = Dict{String, Set{String}}()) = - new(nodes, edges, set_mappings, cofactor_stids, containment) + containment = Dict{String, Set{String}}(), + drug_stids = nothing) = + new(nodes, edges, set_mappings, cofactor_stids, containment, drug_stids) end """ @@ -315,6 +321,7 @@ function parse_complete_network( default_cofactor_path(logic_network_path) : cofactor_path stids = parse_cofactor_list(resolved; required = cofactor_path !== nothing) containment = parse_containment(logic_network_path) + drugs = parse_drug_list(logic_network_path) if isempty(stids) # A file that declares nothing in-network is NOT the same as no file, @@ -335,7 +342,7 @@ function parse_complete_network( end end return ReactionNetwork(network.nodes, network.edges, network.set_mappings, - Set{String}(), containment) + Set{String}(), containment, drugs) end println("Found $(length(stids)) cofactor species declared by the bundle") @@ -364,7 +371,32 @@ function parse_complete_network( end return ReactionNetwork(network.nodes, network.edges, network.set_mappings, stids, - containment) + containment, drugs) +end + +""" +Read the generator's `drugs.csv` beside a logic network (specs/032): the +drug-derived stable ids present in this pathway. `nothing` when there is no +such file (a bundle that predates it), an empty set when the file lists none. +""" +function parse_drug_list(logic_network_path::String)::Union{Nothing, Set{String}} + path = joinpath(dirname(logic_network_path), "drugs.csv") + isfile(path) || return nothing + df = CSV.read(path, DataFrame; types = String) + "stable_id" in names(df) || + throw(ArgumentError("$(path) has no `stable_id` column; it is not a generator drug list")) + return Set(String(x) for x in df.stable_id if !ismissing(x) && !isempty(strip(x))) +end + +"""JSON form of `drug_stids`: `nothing` stays null so "no table" survives a round trip.""" +drug_stids_json(network::ReactionNetwork) = + network.drug_stids === nothing ? nothing : sort(collect(network.drug_stids)) + +"""Inverse of `drug_stids_json`.""" +function drug_stids_from_json(raw)::Union{Nothing, Set{String}} + raw === nothing && return nothing + raw isa AbstractVector || throw(ArgumentError("drug_stids must be a list of stable ids or null")) + return Set(String.(collect(raw))) end """ diff --git a/src/solvers/steady_state.jl b/src/solvers/steady_state.jl index 9d72668..9b845b3 100644 --- a/src/solvers/steady_state.jl +++ b/src/solvers/steady_state.jl @@ -93,6 +93,17 @@ function solve_steady_state( mode = cofactor_mode() cofactors = mode == "inert" ? cofactor_uuids(network) : Set{String}() + # specs/032: drug-derived nodes held at baseline for a cell without the + # drug. Same mechanism as cofactors, own switch and report. + dmode = drug_mode() + drugs = dmode == "inert" ? drug_uuids(network) : Set{String}() + drug_rule = dmode == "propagate" ? "propagate" : + network.drug_stids === nothing ? "inert: no drug table" : "inert" + # specs/035: loop participants released unchanged by their loop, held at + # baseline like cofactors. Resolved on the network as given (before bridges). + cmode = conserved_mode() + conserved = cmode == "inert" ? conserved_uuids(network) : Set{String}() + # Whether to traverse a positional-decomposition silo is a PROCESSING # decision, like the one above: the generator faithfully records that a # curated entity occurs in two places, and this decides whether a signal @@ -103,7 +114,7 @@ function solve_steady_state( @info "Adding $(length(bridges)) silo bridge edge(s)" max_reach=silo_bridge_max_reach() network = ReactionNetwork(network.nodes, vcat(network.edges, bridges), network.set_mappings, network.cofactor_stids, - network.containment) + network.containment, network.drug_stids) end # Convert network to reactions @@ -137,6 +148,32 @@ function solve_steady_state( for u in cofactors if !haskey(observations, u)) observations = merge(observations, pins) end + drugs_held = 0 + if !isempty(drugs) + # An observation that PINS (a measured drug, or a cofactor pin) wins; see + # above. One the confidence gate would discard does not: skipping the + # hold for it let a confidence-0 row silently unpin the drug (review). + # (The cofactor block above still uses haskey; changing it would move + # the default model, so it is left to its own change.) + dpins = Dict{String, Tuple{Float64, Float64}}( + u => (network.nodes[u].baseline * 100.0, 1.0) + for u in drugs + if !(haskey(observations, u) && observations[u][2] > OBS_CONFIDENCE_TOL)) + drugs_held = length(dpins) + observations = merge(observations, dpins) + end + conserved_held = 0 + if !isempty(conserved) + # An observation that PINS wins (a perturbed GAP such as NF1 stays + # perturbed); one the confidence gate would discard does not. + cpins = Dict{String, Tuple{Float64, Float64}}( + u => (network.nodes[u].baseline * 100.0, 1.0) + for u in conserved + if haskey(network.nodes, u) && + !(haskey(observations, u) && observations[u][2] > OBS_CONFIDENCE_TOL)) + conserved_held = length(cpins) + observations = merge(observations, cpins) + end # Initial guess: use observations where available, baseline elsewhere x0 = Dict{String, Float64}() @@ -175,6 +212,10 @@ function solve_steady_state( # Say which model ran: "on", "off", or inert because nothing told the solver # what contains what (a POSTed network or an older bundle). result.diagnostics["self_inhibitor_rule"] = self_rule + result.diagnostics["drug_rule"] = drug_rule + result.diagnostics["drugs_held"] = drugs_held + result.diagnostics["conserved_rule"] = cmode + result.diagnostics["conserved_held"] = conserved_held return result end diff --git a/test/asserting_suites.txt b/test/asserting_suites.txt index 8e20115..6b93149 100644 --- a/test/asserting_suites.txt +++ b/test/asserting_suites.txt @@ -21,3 +21,6 @@ test_solver_determinism test_loop_pool test_scc_break_roles test_self_inhibition +test_drugs +test_set_pool +test_conserved diff --git a/test/test_conserved.jl b/test/test_conserved.jl new file mode 100644 index 0000000..b4e2df9 --- /dev/null +++ b/test/test_conserved.jl @@ -0,0 +1,99 @@ +# specs/035: loop participants released unchanged by their loop (a conserved +# pool, like a cofactor) held at baseline under DS_CONSERVED_MODE=inert. +# +# Fixture = the RAF shape. Loop: F (e.g. F-actin) and the signal S enter R1 +# (AND) -> complex C (contains F) -> R2 (a release: consumes C) -> F and S2 -> +# back. F is released unchanged: CONSERVED. S2 is produced by R2 but C does not +# contain it: transformed, NOT conserved. G is a loop input with an outside +# producer (Z -> G): NOT conserved. +using Test +using DeltaSignal +const DS = DeltaSignal + +function with_env(f, pairs...) + old = Dict(k => get(ENV, k, nothing) for (k, _) in pairs) + for (k, v) in pairs + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + try + f() + finally + for (k, v) in old + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + end +end + +const BL = 0.01 +node(u) = u => DS.NetworkNode(u, "R-" * u, "protein", nothing, u, BL) +edge(s, t, et; and = true) = DS.LogicNetworkEdge(s, t, and, true, 1.0, et) + +function fixture() + nodes = Dict(node.(["X", "F", "S2", "G", "Z", "R1", "C", "R2", "T"])) + edges = [edge("X", "R1", "input"), edge("F", "R1", "input"), edge("S2", "R1", "input"), + edge("G", "R1", "input"), + edge("R1", "C", "output"; and = false), + edge("C", "R2", "input"), + edge("R2", "F", "output"; and = false), edge("R2", "S2", "output"; and = false), + edge("R2", "G", "output"; and = false), edge("Z", "G", "output"; and = false), + edge("C", "T", "input")] + cont = Dict("R-C" => Set(["R-F", "R-G"])) # C holds F and G, not S2 + DS.ReactionNetwork(nodes, edges, Dict{String, DS.SetExpansionMapping}(), Set{String}(), cont) +end + +solve(mode; obs = Dict("X" => (2.0, 1.0), "Z" => (1.0, 1.0))) = + with_env("DS_CONSERVED_MODE" => mode, "DS_SELF_INHIBITOR_WEIGHT" => "off") do + DS.solve_steady_state(fixture(), obs, DS.SteadyStateParams(1.0, 0.1, 500, 1e-6, "penalty")) + end + +@testset "conserved participants (specs/035)" begin + +@testset "detection: released unchanged, all supply in the loop" begin + @test DS.conserved_uuids(fixture()) == Set(["F"]) + # without the containment table nothing is known to be released unchanged + f = fixture() + bare = DS.ReactionNetwork(f.nodes, f.edges, f.set_mappings) + @test isempty(DS.conserved_uuids(bare)) +end + +@testset "detection is label-independent" begin + f = fixture() + rl(u) = "q_" * u + nodes = Dict(rl(k) => DS.NetworkNode(rl(k), v.reactome_id, v.entity_type, nothing, rl(k), v.baseline) + for (k, v) in f.nodes) + edges = [DS.LogicNetworkEdge(rl(e.parent_uuid), rl(e.child_uuid), e.is_and, e.is_positive, + e.stoichiometry, e.edge_type) for e in f.edges] + net = DS.ReactionNetwork(nodes, edges, f.set_mappings, Set{String}(), f.containment) + @test DS.conserved_uuids(net) == Set(["q_F"]) +end + +@testset "default off, byte-identical and reported" begin + with_env("DS_CONSERVED_MODE" => nothing) do + @test DS.conserved_mode() == "off" + end + r = solve(nothing) + @test r.diagnostics["conserved_rule"] == "off" + @test r.diagnostics["conserved_held"] == 0 +end + +@testset "inert holds F at baseline; an observation of F wins" begin + r = solve("inert") + @test r.diagnostics["conserved_rule"] == "inert" + @test r.diagnostics["conserved_held"] == 1 + @test r.node_activities["F"] / BL ≈ 1.0 rtol = 1e-9 + r = solve("inert"; obs = Dict("X" => (2.0, 1.0), "Z" => (1.0, 1.0), "F" => (0.0, 1.0))) + @test r.node_activities["F"] / BL ≈ 0.0 atol = 1e-12 + @test r.diagnostics["conserved_held"] == 0 + r = solve("inert"; obs = Dict("X" => (2.0, 1.0), "Z" => (1.0, 1.0), "F" => (0.0, 0.0))) + @test r.node_activities["F"] / BL ≈ 1.0 rtol = 1e-9 # a gated-out row does not unpin +end + +@testset "a typo is an error" begin + for bad in ("Inert", "on", "") + with_env("DS_CONSERVED_MODE" => bad) do + @test_throws ArgumentError DS.conserved_mode() + end + end +end + +end diff --git a/test/test_drugs.jl b/test/test_drugs.jl new file mode 100644 index 0000000..4178f59 --- /dev/null +++ b/test/test_drugs.jl @@ -0,0 +1,179 @@ +# specs/032: drug-derived entities held at baseline (DS_DRUG_MODE=inert). +# +# Fixture = the RAF shape traced in specs/032: a kinase X and a drug D (a root) +# assemble into the drug-bound complex C, and C inhibits the reaction by which +# X makes its target T (T = X AND W). Left to propagate, raising X raises C, +# which cancels the rise at T: 2x reads 1x. With drugs inert, C stays at +# baseline and T follows X. +using Test +using DeltaSignal +const DS = DeltaSignal + +function with_env(f, pairs...) + old = Dict(k => get(ENV, k, nothing) for (k, _) in pairs) + for (k, v) in pairs + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + try + f() + finally + for (k, v) in old + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + end +end + +const BL = 0.01 +node(u) = u => DS.NetworkNode(u, "R-" * u, "protein", nothing, u, BL) +edge(s, t, pos, et; and = true) = DS.LogicNetworkEdge(s, t, and, pos, 1.0, et) + +function fixture(drug_stids) + nodes = Dict(node.(["X", "D", "C", "W", "T"])) + edges = [edge("X", "C", true, "assembly"), edge("D", "C", true, "assembly"), + edge("X", "T", true, "input"), edge("W", "T", true, "input"), + edge("C", "T", false, "regulator"; and = false)] + DS.ReactionNetwork(nodes, edges, Dict{String, DS.SetExpansionMapping}(), Set{String}(), + Dict{String, Set{String}}(), drug_stids) +end + +function solve(net; mode, x = 2.0, extra = Dict{String, Tuple{Float64, Float64}}()) + with_env("DS_DRUG_MODE" => mode, "DS_SELF_INHIBITOR_WEIGHT" => "off", + "DS_ASSEMBLY_LIMITING" => nothing) do + obs = merge(Dict("X" => (x, 1.0), "W" => (1.0, 1.0)), extra) + DS.solve_steady_state(net, obs, DS.SteadyStateParams(1.0, 0.1, 500, 1e-6, "penalty")) + end +end +fold(r, u) = r.node_activities[u] / BL + +@testset "inert drugs (specs/032)" begin + +@testset "default is propagate, byte-identical to a network with no drug list" begin + with_env("DS_DRUG_MODE" => nothing) do + @test DS.drug_mode() == "propagate" + end + a = solve(fixture(Set(["R-C", "R-D"])); mode = nothing) + b = solve(fixture(nothing); mode = nothing) + @test a.node_activities == b.node_activities + @test a.diagnostics["drug_rule"] == "propagate" + @test a.diagnostics["drugs_held"] == 0 + # The mechanism the rule exists for: the drug-bound complex rises with X + # and its inhibition cancels X's rise at T. + @test fold(a, "C") ≈ 2.0 rtol = 1e-6 + @test fold(a, "T") ≈ 1.0 rtol = 1e-6 +end + +@testset "inert holds drug-derived nodes at baseline and T follows X" begin + for x in (2.0, 0.5, 0.0, 50.0) + r = solve(fixture(Set(["R-C", "R-D"])); mode = "inert", x = x) + @test fold(r, "C") ≈ 1.0 rtol = 1e-9 + @test fold(r, "D") ≈ 1.0 rtol = 1e-9 + @test fold(r, "T") ≈ min(x, 100.0) rtol = 1e-6 + @test r.diagnostics["drug_rule"] == "inert" + @test r.diagnostics["drugs_held"] == 2 + end +end + +@testset "a missing drug table is reported, not pretended" begin + r = solve(fixture(nothing); mode = "inert") + @test r.diagnostics["drug_rule"] == "inert: no drug table" + @test r.diagnostics["drugs_held"] == 0 + @test fold(r, "T") ≈ 1.0 rtol = 1e-6 # nothing held: same as propagate + r = solve(fixture(Set{String}()); mode = "inert") + @test r.diagnostics["drug_rule"] == "inert" # a table listing none is a real answer + @test r.diagnostics["drugs_held"] == 0 +end + +@testset "an explicit observation of a drug wins" begin + r = solve(fixture(Set(["R-C", "R-D"])); mode = "inert", x = 1.0, + extra = Dict("C" => (4.0, 1.0))) + @test fold(r, "C") ≈ 4.0 rtol = 1e-9 + @test fold(r, "T") ≈ 0.25 rtol = 1e-6 + @test r.diagnostics["drugs_held"] == 1 # D only +end + +@testset "a confidence-0 observation does not unpin a drug" begin + # The confidence gate discards it, so it must not stop the hold either. + r = solve(fixture(Set(["R-C", "R-D"])); mode = "inert", x = 2.0, + extra = Dict("C" => (4.0, 0.0))) + @test fold(r, "C") ≈ 1.0 rtol = 1e-9 + @test fold(r, "T") ≈ 2.0 rtol = 1e-6 + @test r.diagnostics["drugs_held"] == 2 +end + +@testset "the silo-bridge rebuild keeps the drug list" begin + # S is one entity split into a receiving copy S1 (X -> S1, no out-edge) and a + # feeding copy S2 (S2 -> T, no in-edge), the shape silo bridges join. Drugs + # are resolved before that rebuild, so this checks the BEHAVIOUR (still held + # with bridges on); the rebuild's own drug_stids pass-through is defensive + # and not read today (a mutation dropping it survives this test). + nodes = Dict(node.(["X", "D", "C", "W", "T"])) + nodes["S1"] = DS.NetworkNode("S1", "R-S", "protein", nothing, "S1", BL) + nodes["S2"] = DS.NetworkNode("S2", "R-S", "protein", nothing, "S2", BL) + edges = [edge("X", "C", true, "assembly"), edge("D", "C", true, "assembly"), + edge("X", "T", true, "input"), edge("W", "T", true, "input"), + edge("C", "T", false, "regulator"; and = false), + edge("X", "S1", true, "input"), edge("S2", "T", true, "input")] + net = DS.ReactionNetwork(nodes, edges, Dict{String, DS.SetExpansionMapping}(), Set{String}(), + Dict{String, Set{String}}(), Set(["R-C", "R-D"])) + with_env("DS_SILO_BRIDGE_MAX_REACH" => "50") do + @test !isempty(DS.silo_bridge_edges(net)) # the rebuild path runs + r = solve(net; mode = "inert") + @test r.diagnostics["drug_rule"] == "inert" + @test r.diagnostics["drugs_held"] == 2 + @test fold(r, "C") ≈ 1.0 rtol = 1e-9 + end +end + +@testset "parse_complete_network reads drugs.csv on both return paths" begin + mktempdir() do dir + write(joinpath(dir, "logic_network.csv"), + "source_id,target_id,pos_neg,and_or,edge_type,stoichiometry\n" * + "u-d,u-c,pos,and,assembly,1\nu-x,u-c,pos,and,assembly,1\n") + write(joinpath(dir, "stid_to_uuid_mapping.csv"), + "uuid,stable_id\nu-d,R-ALL-9\nu-x,R-HSA-1\nu-c,R-HSA-2\n") + ln, um = joinpath(dir, "logic_network.csv"), joinpath(dir, "stid_to_uuid_mapping.csv") + write(joinpath(dir, "drugs.csv"), "stable_id,schema_class,name,reactome_release\nR-ALL-9,ChemicalDrug,d,97\nR-HSA-2,Complex,c,97\n") + # no cofactors.csv: the first return path + net = DS.parse_complete_network(ln, um) + @test net.drug_stids == Set(["R-ALL-9", "R-HSA-2"]) + @test DS.drug_uuids(net) == Set(["u-d", "u-c"]) + # a cofactors.csv declaring one present cofactor: the second return path + write(joinpath(dir, "cofactors.csv"), + "stable_id,molecule,chebi_id,name,in_network,reactome_release\nR-HSA-1,ATP,1,x,1,97\n") + net2 = DS.parse_complete_network(ln, um) + @test net2.cofactor_stids == Set(["R-HSA-1"]) + @test net2.drug_stids == Set(["R-ALL-9", "R-HSA-2"]) + rm(joinpath(dir, "drugs.csv")) + @test DS.parse_complete_network(ln, um).drug_stids === nothing + end +end + +@testset "a typo in the mode is an error" begin + for bad in ("Inert", "inert ", "on", "") + with_env("DS_DRUG_MODE" => bad) do + @test_throws ArgumentError DS.drug_mode() + end + end +end + +@testset "drugs.csv parsing and JSON round trip" begin + mktempdir() do dir + ln = joinpath(dir, "logic_network.csv") + write(ln, "") + @test DS.parse_drug_list(ln) === nothing + write(joinpath(dir, "drugs.csv"), "stable_id,schema_class,name,reactome_release\n") + @test DS.parse_drug_list(ln) == Set{String}() + write(joinpath(dir, "drugs.csv"), + "stable_id,schema_class,name,reactome_release\nR-ALL-1,ChemicalDrug,x,97\nR-HSA-2,Complex,y,97\n") + @test DS.parse_drug_list(ln) == Set(["R-ALL-1", "R-HSA-2"]) + write(joinpath(dir, "drugs.csv"), "id\nR-ALL-1\n") + @test_throws ArgumentError DS.parse_drug_list(ln) + end + @test DS.drug_stids_json(fixture(nothing)) === nothing + @test DS.drug_stids_json(fixture(Set(["R-D", "R-C"]))) == ["R-C", "R-D"] + @test DS.drug_stids_from_json(nothing) === nothing + @test DS.drug_stids_from_json(["R-D"]) == Set(["R-D"]) + @test_throws ArgumentError DS.drug_stids_from_json("R-D") +end + +end diff --git a/test/test_set_pool.jl b/test/test_set_pool.jl new file mode 100644 index 0000000..2e272f3 --- /dev/null +++ b/test/test_set_pool.jl @@ -0,0 +1,107 @@ +# specs/033: a set-valued catalyst or regulator as ONE pool node. +# +# Fixture = the RAF shape: members M1, M2, M3 of a catalyst set feed the pool P +# by `set_member` edges; P catalyses T together with the input X (T = X AND P). +# The pool's mode decides how one member's change reaches T. +using Test +using DeltaSignal +const DS = DeltaSignal + +function with_env(f, pairs...) + old = Dict(k => get(ENV, k, nothing) for (k, _) in pairs) + for (k, v) in pairs + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + try + f() + finally + for (k, v) in old + v === nothing ? delete!(ENV, k) : (ENV[k] = v) + end + end +end + +const BL = 0.01 +node(u) = u => DS.NetworkNode(u, "R-" * u, "protein", nothing, u, BL) +edge(s, t, pos, et; and = true) = DS.LogicNetworkEdge(s, t, and, pos, 1.0, et) + +function fixture() + nodes = Dict(node.(["M1", "M2", "M3", "P", "X", "T", "A", "B", "O"])) + edges = [edge("M1", "P", true, "set_member"; and = false), + edge("M2", "P", true, "set_member"; and = false), + edge("M3", "P", true, "set_member"; and = false), + edge("X", "T", true, "input"), edge("P", "T", true, "catalyst"), + # an ordinary OR node, which is NOT a pool: must not depend on the mode + edge("A", "O", true, "output"; and = false), edge("B", "O", true, "output"; and = false)] + DS.ReactionNetwork(nodes, edges, Dict{String, DS.SetExpansionMapping}()) +end + +function solve(mode; m1 = 1.0, m2 = 1.0, m3 = 1.0, a = 1.0) + with_env("DS_SET_POOL_MODE" => mode, "DS_SELF_INHIBITOR_WEIGHT" => "off") do + obs = Dict("M1" => (m1, 1.0), "M2" => (m2, 1.0), "M3" => (m3, 1.0), + "X" => (1.0, 1.0), "A" => (a, 1.0), "B" => (1.0, 1.0)) + DS.solve_steady_state(fixture(), obs, DS.SteadyStateParams(1.0, 0.1, 500, 1e-6, "penalty")) + end +end +fold(r, u) = r.node_activities[u] / BL + +@testset "set pool nodes (specs/033)" begin + +@testset "set_pool_value, by hand" begin + b = BL + pv(folds, mode) = DS.set_pool_value([f * b for f in folds], b, mode) / b + @test pv([2.0, 3.0, 1.0], "product") ≈ 6.0 + @test pv([0.0, 3.0, 1.0], "product") == 0.0 + @test pv([80.0, 80.0, 1.0], "product") ≈ 100.0 # capped at 100x + @test pv([2.0, 0.25, 1.0], "extreme") ≈ 0.25 # |log .25| > |log 2| + @test pv([2.0, 0.5, 1.0], "extreme") ≈ 0.5 # a tie goes to the lower fold + @test pv([80.0, 1.0, 1.0], "extreme") ≈ 80.0 + @test pv([0.0, 50.0], "extreme") == 0.0 # a knockout wins outright + @test pv([1.0, 1.0, 1.0], "extreme") ≈ 1.0 + @test pv([8.0, 1.0, 1.0], "geomean") ≈ 2.0 + @test pv([0.0, 8.0], "geomean") == 0.0 + @test pv([0.0, 3.0, 1.0], "max") ≈ 3.0 # specs/038: the strongest member + @test pv([0.0, 1.0, 1.0], "max") ≈ 1.0 # a knockout is covered by its peers + @test pv([0.0, 0.0], "max") == 0.0 + @test_throws ArgumentError DS.set_pool_value([b], b, "mean") +end + +@testset "default is product" begin + with_env("DS_SET_POOL_MODE" => nothing) do + @test DS.resolve_reaction_eval_config().set_pool_mode == "product" + end +end + +@testset "one member up 2x, and one knocked out" begin + for (mode, up, ko) in (("product", 2.0, 0.0), ("extreme", 2.0, 0.0), + ("geomean", 2.0^(1/3), 0.0), ("mean", 4/3, 2/3), ("max", 2.0, 1.0)) + r = solve(mode; m1 = 2.0) + @test fold(r, "P") ≈ up rtol = 1e-6 + @test fold(r, "T") ≈ up rtol = 1e-6 + r = solve(mode; m1 = 0.0) + @test fold(r, "P") ≈ ko atol = 1e-9 + @test fold(r, "T") ≈ ko atol = 1e-9 + end +end + +@testset "several members up: product multiplies, extreme and geomean do not" begin + @test fold(solve("product"; m1 = 4.0, m2 = 4.0), "P") ≈ 16.0 rtol = 1e-6 + @test fold(solve("extreme"; m1 = 4.0, m2 = 4.0), "P") ≈ 4.0 rtol = 1e-6 + @test fold(solve("geomean"; m1 = 4.0, m2 = 4.0), "P") ≈ 16.0^(1/3) rtol = 1e-6 +end + +@testset "an ordinary OR node does not depend on the mode" begin + vals = [fold(solve(m; a = 3.0), "O") for m in ("product", "extreme", "geomean", "mean", "max")] + @test all(v -> v ≈ vals[1], vals) + @test vals[1] ≈ 2.0 rtol = 1e-6 # mean(3, 1): the OR rule +end + +@testset "a typo in the mode is an error" begin + for bad in ("Product", "maximum", "") + with_env("DS_SET_POOL_MODE" => bad) do + @test_throws ArgumentError DS.resolve_reaction_eval_config() + end + end +end + +end