From 7f7a326b46bf974ee706d54e1e216309f85454a7 Mon Sep 17 00:00:00 2001 From: Halleluyah Oludele Date: Tue, 25 Aug 2026 20:56:10 +0100 Subject: [PATCH] feat(analysis): block-wise ablation, manuscript preview, repaired binding sites Separate the docking arm into blocks that carry different information rather than ablating it as one lump. Eleven columns describe the docked *ligand's* 3D shape and need no receptor to compute; 132 are receptor-dependent. The combined ablation could not distinguish them, so it could not answer whether receptor structure helps. Measured over three seeds x three folds under stratified, scaffold-grouped and receptor-grouped CV, the receptor-dependent block contributes -0.002 / +0.000 / +0.004 macro-F1 -- each inside the seed spread of its own baseline. Also: - build_preview_html.py renders a reading preview with figures placed at first mention; the DOCX pipeline puts them under the legends, which makes a draft hard to read. - binding_sites.json carries the repaired boxes for the seven receptors whose written chain was a G-protein subunit, a nanobody or a different receptor, plus a repaired_receptors record. main still shipped the stale boxes. Co-Authored-By: Claude Opus 5 --- .gitignore | 1 + data/processed/binding_sites.json | 213 ++++++++++++++++-------------- scripts/ablation_by_block.py | 119 +++++++++++++++++ scripts/build_preview_html.py | 191 +++++++++++++++++++++++++++ 4 files changed, 422 insertions(+), 102 deletions(-) create mode 100644 scripts/ablation_by_block.py create mode 100644 scripts/build_preview_html.py diff --git a/.gitignore b/.gitignore index 0359025..8ac0516 100644 --- a/.gitignore +++ b/.gitignore @@ -169,3 +169,4 @@ scripts/build_natcomms_docx.js scripts/build_imrad_md.py scripts/apa_out.json scripts/figures/ +scripts/check_refs.py diff --git a/data/processed/binding_sites.json b/data/processed/binding_sites.json index 9561ec2..2457397 100644 --- a/data/processed/binding_sites.json +++ b/data/processed/binding_sites.json @@ -22,25 +22,25 @@ "uniprot": "O14842" }, { + "uniprot": "O15552", "biasdb_name": "FFA2 receptor", - "center_x": 121.3982, - "center_y": 114.1356, - "center_z": 113.9273, - "cocrystal_ligand_n_atoms": null, - "cocrystal_ligand_resname": null, - "confidence": "ok", + "pdb_id": "8T3S", + "structure_source": "pdb", "method": "consensus", - "notes": "p2rank top + fpocket top agreement = 2.19 \u00c5", - "pdb_id": "9K1D", - "pocket_prediction_agrees": true, - "pocket_prediction_distance_angstroms": 2.188294288573317, - "pocket_score_fpocket": 0.385, - "pocket_score_p2rank": 1.89, + "center_x": 96.4309, + "center_y": 96.1629, + "center_z": 138.8659, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "O15552" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 4.706148663426533, + "pocket_prediction_agrees": true, + "pocket_score_p2rank": 16.99, + "pocket_score_fpocket": 0.01, + "confidence": "ok", + "notes": "p2rank top + fpocket top agreement = 4.71 \u00c5" }, { "biasdb_name": "\u03b22-adrenoceptor", @@ -253,25 +253,25 @@ "uniprot": "P20309" }, { + "uniprot": "P21452", "biasdb_name": "NK2 receptor", - "center_x": 135.1467, - "center_y": 154.615, - "center_z": 106.6592, - "cocrystal_ligand_n_atoms": null, - "cocrystal_ligand_resname": null, - "confidence": "ok", - "method": "consensus", - "notes": "p2rank top + fpocket top agreement = 1.54 \u00c5", "pdb_id": "7XWO", - "pocket_prediction_agrees": true, - "pocket_prediction_distance_angstroms": 1.5371655525782293, - "pocket_score_fpocket": 0.697, - "pocket_score_p2rank": 34.42, + "structure_source": "pdb", + "method": "primary_only", + "center_x": 114.1459, + "center_y": 103.232, + "center_z": 143.3739, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P21452" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 9.265273875342334, + "pocket_prediction_agrees": false, + "pocket_score_p2rank": 35.31, + "pocket_score_fpocket": 0.024, + "confidence": "low_confidence", + "notes": "p2rank top + fpocket top agreement = 9.27 \u00c5" }, { "biasdb_name": "S1P1 receptor", @@ -295,25 +295,25 @@ "uniprot": "P21453" }, { + "uniprot": "P21462", "biasdb_name": "FPR1", - "center_x": 129.6976, - "center_y": 85.7293, - "center_z": 119.9248, - "cocrystal_ligand_n_atoms": null, - "cocrystal_ligand_resname": null, - "confidence": "ok", - "method": "consensus", - "notes": "p2rank top + fpocket top agreement = 0.91 \u00c5", - "pdb_id": "7VFX", - "pocket_prediction_agrees": true, - "pocket_prediction_distance_angstroms": 0.9125193969244265, - "pocket_score_fpocket": 0.507, - "pocket_score_p2rank": 6.55, + "pdb_id": "7EUO", + "structure_source": "pdb", + "method": "primary_only", + "center_x": 108.4012, + "center_y": 115.0815, + "center_z": 76.1135, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P21462" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 10.812407348402004, + "pocket_prediction_agrees": false, + "pocket_score_p2rank": 24.05, + "pocket_score_fpocket": 0.002, + "confidence": "low_confidence", + "notes": "p2rank top + fpocket top agreement = 10.81 \u00c5" }, { "biasdb_name": "CB1 receptor", @@ -547,25 +547,25 @@ "uniprot": "P29275" }, { + "uniprot": "P30518", "biasdb_name": "V2 receptor", - "center_x": 98.2838, - "center_y": 112.9913, - "center_z": 127.3186, - "cocrystal_ligand_n_atoms": null, - "cocrystal_ligand_resname": null, - "confidence": "low_confidence", + "pdb_id": "7DW9", + "structure_source": "pdb", "method": "primary_only", - "notes": "p2rank top + fpocket top agreement = 13.55 \u00c5", - "pdb_id": "7KH0", - "pocket_prediction_agrees": false, - "pocket_prediction_distance_angstroms": 13.545152622333719, - "pocket_score_fpocket": 0.084, - "pocket_score_p2rank": 1.53, + "center_x": 49.3278, + "center_y": 70.0202, + "center_z": 72.9365, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P30518" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 7.609059392513593, + "pocket_prediction_agrees": false, + "pocket_score_p2rank": 16.23, + "pocket_score_fpocket": 0.144, + "confidence": "low_confidence", + "notes": "p2rank top + fpocket top agreement = 7.61 \u00c5" }, { "biasdb_name": "A1 receptor", @@ -820,25 +820,25 @@ "uniprot": "P41143" }, { + "uniprot": "P41145", "biasdb_name": "\u03ba receptor", - "center_x": 154.11657097226097, - "center_y": 141.57423800513857, - "center_z": 122.21078582037063, - "cocrystal_ligand_n_atoms": 42, - "cocrystal_ligand_resname": "Q6Q", - "confidence": "ok", + "pdb_id": "8DZP", + "structure_source": "pdb", "method": "cocrystal_ligand", - "notes": "box centered on bounded Q6Q (42 heavy atoms)", - "pdb_id": "7UL2", - "pocket_prediction_agrees": null, - "pocket_prediction_distance_angstroms": null, - "pocket_score_fpocket": null, - "pocket_score_p2rank": null, + "center_x": 129.84212838449787, + "center_y": 120.11083910542149, + "center_z": 93.83938672465663, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P41145" + "cocrystal_ligand_resname": "U99", + "cocrystal_ligand_n_atoms": 31, + "pocket_prediction_distance_angstroms": null, + "pocket_prediction_agrees": null, + "pocket_score_p2rank": null, + "pocket_score_fpocket": null, + "confidence": "ok", + "notes": "box centered on bounded U99 (31 heavy atoms)" }, { "biasdb_name": "NOP receptor", @@ -904,25 +904,25 @@ "uniprot": "P41595" }, { + "uniprot": "P41597", "biasdb_name": "CCR2", - "center_x": 3.083374999463558, - "center_y": -6.391812473535538, - "center_z": -23.6166250705719, - "cocrystal_ligand_n_atoms": 16, - "cocrystal_ligand_resname": "TYS", - "confidence": "ok", - "method": "cocrystal_ligand", - "notes": "box centered on bounded TYS (16 heavy atoms)", - "pdb_id": "7P8X", - "pocket_prediction_agrees": null, - "pocket_prediction_distance_angstroms": null, - "pocket_score_fpocket": null, - "pocket_score_p2rank": null, + "pdb_id": "7XA3", + "structure_source": "pdb", + "method": "consensus", + "center_x": 113.6946, + "center_y": 138.318, + "center_z": 130.1826, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P41597" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 4.824759347532267, + "pocket_prediction_agrees": true, + "pocket_score_p2rank": 23.91, + "pocket_score_fpocket": 0.427, + "confidence": "ok", + "notes": "p2rank top + fpocket top agreement = 4.82 \u00c5" }, { "biasdb_name": "MC3 receptor", @@ -988,25 +988,25 @@ "uniprot": "P49682" }, { + "uniprot": "P51681", "biasdb_name": "CCR5", - "center_x": -3.751437494531274, - "center_y": -11.872749924659729, - "center_z": -18.457124888896942, - "cocrystal_ligand_n_atoms": 16, - "cocrystal_ligand_resname": "TYS", - "confidence": "ok", - "method": "cocrystal_ligand", - "notes": "box centered on bounded TYS (16 heavy atoms)", - "pdb_id": "5YY4", - "pocket_prediction_agrees": null, - "pocket_prediction_distance_angstroms": null, - "pocket_score_fpocket": null, - "pocket_score_p2rank": null, + "pdb_id": "7O7F", + "structure_source": "pdb", + "method": "primary_only", + "center_x": 192.7681, + "center_y": 179.8396, + "center_z": 230.1375, "size_x": 22.0, "size_y": 22.0, "size_z": 22.0, - "structure_source": "pdb", - "uniprot": "P51681" + "cocrystal_ligand_resname": null, + "cocrystal_ligand_n_atoms": null, + "pocket_prediction_distance_angstroms": 13.56752954970882, + "pocket_prediction_agrees": false, + "pocket_score_p2rank": 26.63, + "pocket_score_fpocket": 0.002, + "confidence": "low_confidence", + "notes": "p2rank top + fpocket top agreement = 13.57 \u00c5" }, { "biasdb_name": "PAR2", @@ -1298,5 +1298,14 @@ "primary_only": 2 }, "receptors": 61 - } + }, + "repaired_receptors": [ + "O15552", + "P21452", + "P21462", + "P30518", + "P41145", + "P41597", + "P51681" + ] } \ No newline at end of file diff --git a/scripts/ablation_by_block.py b/scripts/ablation_by_block.py new file mode 100644 index 0000000..06cde2c --- /dev/null +++ b/scripts/ablation_by_block.py @@ -0,0 +1,119 @@ +"""Block-wise ablation: separate ligand 3D shape from receptor-derived features. + +`reviewer_extras.ablation_drop_structural` drops one combined "structural" +block, but that block mixes three different information sources: descriptors of +the docked *ligand's* 3D shape (which need no receptor), pair-level docking +energetics, and receptor-contact fingerprints. A ligand is a structure too, so +the combined Delta cannot answer "does the receptor help?" — only "does the +whole docking arm help?". This runs each block separately. + +Run: + PYTHONPATH=src python scripts/ablation_by_block.py +""" +from __future__ import annotations + +import json as _json +import warnings +from pathlib import Path + +import joblib +import numpy as np +import pandas as pd +from sklearn.metrics import f1_score +from sklearn.model_selection import GroupKFold, StratifiedKFold + +from cancerag.ml.preprocessing import get_X_y_groups, build_full_pipeline +from cancerag.ml.model_training import MODEL_FACTORIES, _combined_weight + +warnings.filterwarnings("ignore") + +POSE_3D = ("Asphericity", "Eccentricity", "InertialShapeFactor", "NPR1", "NPR2", + "PMI1", "PMI2", "PMI3", "RadiusOfGyration", "SpherocityIndex", + "pose_3d_missing") +# needs the receptor to compute at all +RECEPTOR_DEPENDENT = ("vina_", "ifp_", "n_residues_contacted", "n_total_contacts", + "n_poses", "ifp_missing", "ifp_no_contacts") + +SEEDS = (42, 7, 13) +FOLDS = 3 + + +def cols_matching(X, prefixes): + return [c for c in X.columns + if any(c.startswith(p) or c == p for p in prefixes)] + + +def main() -> None: + df = pd.read_parquet("data/processed/ml_ready_dataset.parquet") + le = joblib.load("data/processed/ml_preprocessed/label_encoder.joblib") + X, y, sw, le, _ = get_X_y_groups(df, label_encoder=le) + n_classes = len(le.classes_) + winner = _json.loads( + Path("data/processed/ml_models/selection_decision.json").read_text())["chosen"] + factory = MODEL_FACTORIES[winner] + + pose = cols_matching(X, POSE_3D) + recep = cols_matching(X, RECEPTOR_DEPENDENT) + print(f"model={winner} total={X.shape[1]} " + f"ligand-3D-shape={len(pose)} receptor-dependent={len(recep)}") + + # The published ablation ran StratifiedKFold while the manuscript described + # it as scaffold-grouped. Stratified is the optimistic regime this paper + # argues against, so the headline Delta must be reported on the grouped + # splits as well -- especially receptor-grouped, which is the only split + # that asks whether receptor features help on an unseen receptor. + variants = { + "full": X, + "drop_receptor_dependent": X.drop(columns=recep), # keeps ligand 3D shape + "drop_ligand_3d_shape": X.drop(columns=pose), # keeps receptor features + "drop_both": X.drop(columns=recep + pose), + } + + scaf = pd.factorize(df["murcko_scaffold"])[0] if "murcko_scaffold" in df.columns \ + else pd.factorize(df["scaffold"])[0] + recep = pd.factorize(df["receptor_uniprot"])[0] + SPLITS = { + "stratified": lambda Xv, seed: StratifiedKFold( + n_splits=FOLDS, shuffle=True, random_state=seed).split(Xv, y), + "scaffold_grouped": lambda Xv, seed: GroupKFold(n_splits=FOLDS).split(Xv, y, scaf), + "receptor_grouped": lambda Xv, seed: GroupKFold(n_splits=FOLDS).split(Xv, y, recep), + } + + rows = [] + for split_name, splitter in SPLITS.items(): + for name, Xv in variants.items(): + for seed in SEEDS: + f1s = [] + for tr, te in splitter(Xv, seed): + model = factory(n_classes=n_classes, random_state=seed) + pipe = build_full_pipeline(model) + last = pipe.steps[-1][0] + pipe.fit(Xv.iloc[tr], y[tr], + **{f"{last}__sample_weight": _combined_weight(y[tr], sw[tr])}) + f1s.append(float(f1_score(y[te], pipe.predict(Xv.iloc[te]), + average="macro", zero_division=0))) + rows.append({"split": split_name, "variant": name, "seed": seed, + "n_features": int(Xv.shape[1]), + "macro_f1": float(np.mean(f1s))}) + print(f" [{split_name}] done") + + out = pd.DataFrame(rows) + outdir = Path("data/processed/ml_models/extras"); outdir.mkdir(parents=True, exist_ok=True) + out.to_csv(outdir / "ablation_by_block.csv", index=False) + + print("\n=== mean macro-F1 over 3 seeds x 3 folds ===") + for split_name in SPLITS: + sub = out[out.split == split_name] + base = sub[sub.variant == "full"].macro_f1 + print(f"\n--- {split_name} ---") + for name in variants: + m = sub[sub.variant == name].macro_f1 + d = m.mean() - base.mean() + # seed spread of the baseline, for judging whether Delta is signal + print(f" {name:<26} {m.mean():.4f} +/- {m.std():.4f}" + f" (Delta vs full: {d:+.4f})") + print(f" [baseline seed spread: +/-{base.std():.4f}]") + + +if __name__ == "__main__": + main() diff --git a/scripts/build_preview_html.py b/scripts/build_preview_html.py new file mode 100644 index 0000000..a773bb2 --- /dev/null +++ b/scripts/build_preview_html.py @@ -0,0 +1,191 @@ +"""Build a self-contained HTML reading preview with figures placed inline. + +The DOCX builder embeds figures only under the figure-legend sections at the +end, which is the Nature submission convention but makes the manuscript hard to +read: a Results paragraph cites "Fig. 3" and the reader has to scroll past the +Methods to see it. The markdown carries no image tags at all, so a plain +markdown preview shows no figures whatever. + +This renderer places each figure at its first in-text mention in Results, keeps +the legend text with it, and inlines every image as a data URI so the file works +with no network and no sibling asset directory. + +Run: + PYTHONPATH=src python scripts/build_preview_html.py +""" +from __future__ import annotations + +import base64 +import re +import sys +from pathlib import Path + +import mistune + +REPO = Path(__file__).resolve().parent.parent +FIG_DIR = REPO / "manuscript" / "figures" + +def figure_files(md: str) -> dict[int, list[str]]: + """Figure number -> image files, read from each legend's tag. + + The mapping travels with the legend so a renumbered document stays correct + without a second table to keep in step. + """ + out: dict[int, list[str]] = {} + for m in re.finditer(r"\*\*Figure (\d+) \|(?:(?!\*\*Figure )[\s\S])*?", md): + out[int(m.group(1))] = [f.strip() for f in m.group(2).split(",") if f.strip()] + return out + + +def data_uri(name: str) -> str | None: + p = FIG_DIR / name + if not p.exists(): + return None + return "data:image/png;base64," + base64.b64encode(p.read_bytes()).decode() + + +def legends(md: str) -> dict[int, tuple[str, str]]: + """Figure number -> (title, remaining legend text).""" + out: dict[int, tuple[str, str]] = {} + for m in re.finditer(r"^\*\*Figure (\d+) \|(.+?)\*\*(.*?)$", md, re.M | re.S): + out[int(m.group(1))] = (m.group(2).strip(), m.group(3).strip()) + return out + + +def figure_block(num: int, leg: dict, files: dict) -> str: + imgs = [u for u in (data_uri(f) for f in files.get(num, [])) if u] + if not imgs: + return "" + title, rest = leg.get(num, ("", "")) + cap = f"Figure {num}" + if title: + cap += f" | {mistune.html(title).strip()[3:-4]}" + rest = re.sub(r"", "", rest).strip() + body = f"
{mistune.html(rest)}
" if rest else "" + tags = "".join(f"Figure {num}" for u in imgs) + return (f"
{tags}" + f"
{cap}{body}
") + + +def main() -> None: + src = Path(sys.argv[1]) + dst = Path(sys.argv[2]) + md = src.read_text() + leg = legends(md) + files = figure_files(md) + + # Render body and legend sections separately: the legend sections keep their + # figures too, so the end-of-paper layout a reviewer expects still works. + html = mistune.html(re.sub(r"", "", md)) + + # Place each figure at its first in-text citation. Anchor on the closing tag + # of the paragraph containing the mention so a figure never lands mid-sentence. + placed: set[int] = set() + + def place(match: re.Match) -> str: + para = match.group(0) + add = [] + for mm in re.finditer(r"Fig(?:ure)?\.?\s*(\d+)", para): + num = int(mm.group(1)) + if num in placed: + continue + blk = figure_block(num, leg, files) + if blk: + placed.add(num) + add.append(blk) + return para + "".join(add) + + # Paragraphs and headings, and only before the legend sections. Headings + # must be scanned too: several Results sections cite their figure only in + # the heading ("... (Fig. 5)"), so a paragraph-only pass drops the two + # main figures the paper is built around. + cut = html.find("

Figure legends") + head, tail = (html[:cut], html[cut:]) if cut != -1 else (html, "") + head = re.sub(r"<(p|h2|h3)>(?:(?!).)*?", place, head, flags=re.S) + + # legend sections: show the image under each legend as well + def legend_img(match: re.Match) -> str: + para = match.group(0) + mm = re.search(r"Figure (\d+) \|", para) or re.search(r"Figure (\d+)", para) + if not mm: + return para + num = int(mm.group(1)) + tags = "".join(f"" for u in + (data_uri(f) for f in files.get(num, [])) if u) + return para + (f"
{tags}
" if tags else "") + + tail = re.sub(r"

(?:(?!

).)*?

", legend_img, tail, flags=re.S) + + title = re.search(r"^# (.+)$", md, re.M) + title = title.group(1) if title else src.stem + + dst.write_text(TEMPLATE.replace("{{TITLE}}", title) + .replace("{{BODY}}", head + tail) + .replace("{{NFIG}}", str(len(placed)))) + print(f"{dst} ({dst.stat().st_size/1_000_000:.1f} MB, {len(placed)} figures placed inline)") + + +TEMPLATE = """ + + +{{TITLE}} +
+ +{{BODY}} +
+ + +""" + +if __name__ == "__main__": + main()