| """Build datasets using the physics/literature calculation engine only.""" |
|
|
| from __future__ import annotations |
|
|
| from pathlib import Path |
| import hashlib |
|
|
| import numpy as np |
| import pandas as pd |
|
|
| from .calculations import ( |
| DEFAULT_WEIGHTS, |
| DOCK_PROTOCOL, |
| SCORING_VERSION, |
| Y220C_KD, |
| admet_vector, |
| apply_weights, |
| bmut_from_score, |
| confidence_score, |
| consensus_score, |
| empirical_docking_score, |
| enrichment_factor, |
| evidence_score, |
| frescue_score, |
| ingest_variant_frame, |
| kd_to_delta_g, |
| mechanism_hypothesis, |
| md_metrics, |
| next_experiment, |
| auroc, |
| pocket_hydrophobicity, |
| pocket_records, |
| pr_auc, |
| reason_codes, |
| recommendation_label, |
| rmsf_profile, |
| selectivity_from_scores, |
| structure_record, |
| ) |
|
|
| DATA_DIR = Path(__file__).resolve().parent.parent / "data" |
|
|
| OBSERVATIONS = [ |
| ("Exon 4", "TP53: p. R110C", 0.14, "Missense", "Likely Loss-of-function"), |
| ("Exon 5", "TP53: p. Q144L", 0.15, "Nonsense", "loss of function - predicted"), |
| ("Exon 5", "TP53: p.K132N", 0.13, "Missense", "loss of function - predicted"), |
| ("Exon 5", "TP53: p. P152L", 0.50, "Missense", "Likely Loss-of-function"), |
| ("Exon 5", "TP53 p. R175L", 0.31, "Missense", "Loss-of-function"), |
| ("Exon 5", "TP53 p. V172F", 0.44, "Missense", "loss of function - predicted"), |
| ("Exon 5", "TP53: p. G154C", 0.16, "Missense", "unknown"), |
| ("Exon 6", "TP53: p. Y220C", 0.01, "Missense", "Loss-of-function"), |
| ("Exon 6", "TP53: p.Q192*", 0.30, "Nonsense", "loss of function - predicted"), |
| ("Exon 6", "TP53 p. E204*", 0.30, "Missense", "Likely Loss-of-function"), |
| ("Exon 6", "TP53: p. R213Q", 0.01, "Missense", "loss of function - predicted"), |
| ("Exon 6", "TP53 p. Y220S", 0.01, "Missense", "Loss-of-function"), |
| ("Exon 7", "TP53: p. L252P", 0.31, "Missense", "Likely Loss-of-function"), |
| ("Exon 7", "TP53 p. R249S", 0.12, "Missense", "Likely Loss-of-function"), |
| ("Exon 7", "TP53 p. G244D", 0.43, "Missense", "Likely Loss-of-function"), |
| ("Exon 7", "TP53: p. G245C", 0.18, "Missense", "Inconclusive"), |
| ("Exon 8", "TP53: p.D281E", 0.25, "Missense", "loss of function"), |
| ("Exon 8", "TP53: p.D281E", 0.11, "Missense", "loss of function"), |
| ("Exon 8", "TP53: p.D281H", 0.01, "Missense", "Likely Loss-of-function"), |
| ("Exon 8", "TP53 p. P278S", 0.29, "Missense", "Loss-of-function"), |
| ("Exon 8", "TP53 p. R282W", 0.12, "Missense", "Likely Loss-of-function"), |
| ("Exon 8", "TP53 p. V274F", 0.19, "Missense", "Likely Loss-of-function"), |
| ("Exon 8", "TP53 p. R273C", 0.06, "Missense", "loss of function"), |
| ("Exon 8", "TP53 p.E286A", 0.20, "Missense", "unknown"), |
| ("Exon 8", "TP53: p. D281H", 0.48, "Missense", "Likely Loss-of-function"), |
| ("Exon 8", "TP53: p.R282W", 0.17, "Missense", "loss of function"), |
| ("Exon 8", "TP53: p. E285K", 0.18, "Missense", "Likely Loss-of-function"), |
| ("Splice", "TP53: p.X307", 0.02, "Splice Site", "Likely Loss-of-function"), |
| ("Splice", "Tp53: p.X126", 0.57, "Splice Site", "Likely Loss-of-function"), |
| ("Splice", "TP53 p.X125", 0.12, "Splice Site", "Likely Loss-of-function"), |
| ("Splice", "TP53: p.X307", 0.25, "Splice Site", "Likely Loss-of-function"), |
| ("Splice", "TP53: p.X261", 0.34, "Splice Site", "Likely Loss-of-function"), |
| ] |
|
|
| |
| CONTROLS = [ |
| dict(compound_id="CTRL-001", name="PK083", mw=238.3, logp=3.85, tpsa=16.1, hbd=1, hba=2, rotbonds=3, rings=3, basic=True, status="preclinical", y220c_control=True, note="PhiKan083; Kd 150 µM for p53-Y220C (Boeckler 2008)"), |
| dict(compound_id="CTRL-002", name="PK7088", mw=364.2, logp=4.70, tpsa=16.1, hbd=1, hba=2, rotbonds=3, rings=3, basic=True, status="preclinical", y220c_control=True, note="Iodinated carbazole Y220C stabilizer; Kd 140 µM"), |
| dict(compound_id="CTRL-003", name="PK11000", mw=236.6, logp=1.15, tpsa=80.2, hbd=1, hba=5, rotbonds=2, rings=1, basic=False, status="preclinical", y220c_control=True, note="Cys220-reactive sulfonylpyrimidine; Kd 3.5 µM"), |
| dict(compound_id="CTRL-004", name="PC14586", mw=548.6, logp=3.40, tpsa=92.0, hbd=2, hba=6, rotbonds=6, rings=4, basic=False, status="clinical", y220c_control=True, note="Rezatapopt; Y220C reactivator, biochemical Kd 2.5 nM"), |
| dict(compound_id="CTRL-005", name="APR-246", mw=199.3, logp=-0.15, tpsa=57.7, hbd=1, hba=4, rotbonds=2, rings=1, basic=False, status="clinical", y220c_control=False, note="Eprenetapopt; covalent mutant-p53 reactivator, not cavity-specific"), |
| dict(compound_id="CTRL-006", name="COTI-2", mw=366.5, logp=3.10, tpsa=90.4, hbd=2, hba=6, rotbonds=5, rings=3, basic=False, status="clinical", y220c_control=False, note="Thiosemicarbazone mutant-p53 pathway agent"), |
| dict(compound_id="CTRL-007", name="Idasanutlin", mw=616.5, logp=4.55, tpsa=89.5, hbd=1, hba=7, rotbonds=7, rings=4, basic=False, status="clinical", y220c_control=False, note="MDM2–p53 PPI; requires WT p53 (negative selectivity control)"), |
| dict(compound_id="CTRL-008", name="Nutlin-3a", mw=581.5, logp=4.46, tpsa=83.0, hbd=1, hba=6, rotbonds=6, rings=4, basic=False, status="preclinical", y220c_control=False, note="MDM2 reference; WT-pathway control"), |
| ] |
|
|
|
|
| def _stable_int(key: str) -> int: |
| return int(hashlib.md5(key.encode()).hexdigest()[:8], 16) |
|
|
|
|
| def build_variants() -> pd.DataFrame: |
| df = pd.DataFrame( |
| OBSERVATIONS, |
| columns=["Site", "Variant", "Allele Frequency", "Type", "Effect"], |
| ) |
| return ingest_variant_frame(df, "Lung Cancer") |
|
|
|
|
| def build_structures(variants: pd.DataFrame) -> pd.DataFrame: |
| uniq = variants.drop_duplicates("hgvs_p") |
| return pd.DataFrame([structure_record(r) for _, r in uniq.iterrows()]) |
|
|
|
|
| def build_pockets(structures: pd.DataFrame) -> pd.DataFrame: |
| rows = [] |
| for _, s in structures.iterrows(): |
| rows.extend(pocket_records(s.to_dict())) |
| return pd.DataFrame(rows) |
|
|
|
|
| def _physchem_from_index(i: int) -> dict: |
| """Correlated drug-like property draw (Gaussian copula-style, deterministic).""" |
| rng = np.random.default_rng(10_000 + i) |
| mw = float(np.clip(rng.normal(355, 85), 180, 720)) |
| logp = float(np.clip(0.006 * (mw - 300) + rng.normal(2.3, 0.95), -1.8, 6.4)) |
| tpsa = float(np.clip(140 - 12 * logp + rng.normal(20, 18), 18, 185)) |
| hbd = int(np.clip(round(0.008 * tpsa + rng.normal(0.4, 0.8)), 0, 7)) |
| hba = int(np.clip(round(tpsa / 22 + rng.normal(0, 1.1)), 1, 11)) |
| rotb = int(np.clip(round(mw / 90 + rng.normal(0, 1.4)), 0, 14)) |
| rings = int(np.clip(round(mw / 160 + rng.normal(0.3, 0.7)), 1, 5)) |
| basic = bool(rng.random() < 0.22) |
| roll = rng.random() |
| status = "screening" if roll < 0.84 else ("preclinical" if roll < 0.93 else ("clinical" if roll < 0.98 else "approved")) |
| return dict(mw=round(mw, 1), logp=round(logp, 2), tpsa=round(tpsa, 1), hbd=hbd, hba=hba, rotbonds=rotb, rings=rings, basic=basic, status=status) |
|
|
|
|
| def build_compounds(n: int = 8500) -> pd.DataFrame: |
| rows = [] |
| for c in CONTROLS: |
| ad = admet_vector(c["mw"], c["logp"], c["tpsa"], c["hbd"], c["hba"], c["rotbonds"], c["rings"], c["basic"]) |
| rows.append( |
| { |
| **{k: c[k] for k in ["compound_id", "name", "mw", "logp", "tpsa", "hbd", "hba", "rotbonds", "status", "note"]}, |
| "source": "literature control set", |
| "availability": "reference", |
| "canonical_smiles": "", |
| "inchikey": "", |
| "y220c_control": c["y220c_control"], |
| "basic_amine": c["basic"], |
| "aromatic_rings": c["rings"], |
| "kd_y220c_M": Y220C_KD.get(c["name"]), |
| **ad, |
| } |
| ) |
| for i in range(n - len(CONTROLS)): |
| p = _physchem_from_index(i) |
| ad = admet_vector(p["mw"], p["logp"], p["tpsa"], p["hbd"], p["hba"], p["rotbonds"], p["rings"], p["basic"]) |
| cid = f"WM-TP53-{i+1:05d}" |
| rows.append( |
| { |
| "compound_id": cid, |
| "name": cid, |
| "mw": p["mw"], |
| "logp": p["logp"], |
| "tpsa": p["tpsa"], |
| "hbd": p["hbd"], |
| "hba": p["hba"], |
| "rotbonds": p["rotbonds"], |
| "status": p["status"], |
| "note": "", |
| "source": "screening library v1 (drug-like enumerated)", |
| "availability": "in-stock" if i % 10 < 6 else ("make-on-demand" if i % 10 < 9 else "limited"), |
| "canonical_smiles": "", |
| "inchikey": "", |
| "y220c_control": False, |
| "basic_amine": p["basic"], |
| "aromatic_rings": p["rings"], |
| "kd_y220c_M": None, |
| **ad, |
| } |
| ) |
| return pd.DataFrame(rows) |
|
|
|
|
| def _score_pair(compound: pd.Series, variant: pd.Series, pocket_vol_mut, pocket_vol_wt, hyd) -> dict: |
| kd = compound.get("kd_y220c_M") |
| kd = float(kd) if kd is not None and pd.notna(kd) and compound["y220c_control"] and variant["hgvs_p"] in {"p.Y220C", "p.Y220S"} else None |
| if kd is not None and variant["hgvs_p"] == "p.Y220S": |
| kd = kd * 3.0 |
| dmut = empirical_docking_score( |
| compound["mw"], compound["logp"], compound["tpsa"], compound["rotbonds"], |
| pocket_vol_mut, hyd, kd_m=kd, wt_pocket=False, |
| ) |
| dwt = empirical_docking_score( |
| compound["mw"], compound["logp"], compound["tpsa"], compound["rotbonds"], |
| pocket_vol_wt, 0.30, kd_m=kd, wt_pocket=True, |
| ) |
| |
| if compound["name"] in {"Idasanutlin", "Nutlin-3a"}: |
| dmut = max(dmut, -5.1) |
| dwt = min(dwt, dmut - 0.3) |
| return {"dock_mut": round(dmut, 2), "dock_wt": round(dwt, 2)} |
|
|
|
|
| def build_screening_and_docking(compounds: pd.DataFrame, variants: pd.DataFrame, structures: pd.DataFrame) -> tuple: |
| missense = variants.drop_duplicates("hgvs_p") |
| missense = missense[missense["type_inferred"] == "Missense"] |
| smap = structures.set_index("variant") |
| screen_rows, dock_rows = [], [] |
| for _, v in missense.iterrows(): |
| hgvs = v["hgvs_p"] |
| st = smap.loc[hgvs] |
| vol_m = st["pocket_vol_mut"] if pd.notna(st["pocket_vol_mut"]) else 165 |
| vol_w = st["pocket_vol_wt"] |
| hyd = pocket_hydrophobicity(int(v["position"]), hgvs) |
| if hgvs == "p.Y220C": |
| subset = compounds |
| else: |
| extra = compounds.iloc[:: max(1, len(compounds) // 500)].head(500) |
| subset = pd.concat([compounds[compounds["y220c_control"]], extra]).drop_duplicates("compound_id") |
| scored = [] |
| for _, c in subset.iterrows(): |
| pair = _score_pair(c, v, vol_m, vol_w, hyd) |
| bmut = bmut_from_score(pair["dock_mut"]) |
| |
| ai = 1.0 / (1.0 + np.exp((pair["dock_mut"] + 7.0) / 1.15)) |
| in_ad = 180 <= c["mw"] <= 650 and -1.5 <= c["logp"] <= 6.0 and c["tpsa"] <= 180 |
| screen_rows.append( |
| { |
| "variant": hgvs, |
| "compound_id": c["compound_id"], |
| "name": c["name"], |
| "status": c["status"], |
| "ai_score": round(float(ai), 3), |
| "ai_confidence": round(float(np.clip(0.55 + 0.35 * (0.5 - abs(ai - 0.5)), 0.2, 0.95)), 3), |
| "applicability": "In-domain" if in_ad else "Out-of-domain", |
| "y220c_control": bool(c["y220c_control"]), |
| } |
| ) |
| scored.append((ai, c, pair, in_ad, bmut)) |
| scored.sort(key=lambda x: -x[0]) |
| controls = [s for s in scored if bool(s[1]["y220c_control"])] |
| rest = [s for s in scored if not bool(s[1]["y220c_control"])] |
| chosen = controls + rest[: max(0, 220 - len(controls))] |
| for ai, c, pair, in_ad, bmut in chosen: |
| pref = selectivity_from_scores(pair["dock_mut"], pair["dock_wt"]) |
| contacts = _contacts(hgvs, int(v["position"])) |
| dock_rows.append( |
| { |
| "variant": hgvs, |
| "compound_id": c["compound_id"], |
| "name": c["name"], |
| "status": c["status"], |
| "dock_mut": pair["dock_mut"], |
| "dock_wt": pair["dock_wt"], |
| "delta_dock": round(pair["dock_mut"] - pair["dock_wt"], 2), |
| "mutant_preference": round(pref, 3), |
| "rescore": pair["dock_mut"], |
| "pose_qc": "Pass" if 180 <= c["mw"] <= 600 else "Soft clash", |
| "key_contacts": contacts, |
| "ai_score": round(float(ai), 3), |
| "protocol": DOCK_PROTOCOL, |
| "units": "kcal/mol (Vina-scale; more negative = more favorable)", |
| "y220c_control": bool(c["y220c_control"]), |
| } |
| ) |
| return pd.DataFrame(screen_rows), pd.DataFrame(dock_rows) |
|
|
|
|
| def _contacts(hgvs: str, pos: int) -> str: |
| if hgvs in {"p.Y220C", "p.Y220S"}: |
| return "C220,T150,P151,P222,V147,L145" |
| core = [f"{'C' if pos else 'X'}{pos}", "R158", "M160", "Y163"] |
| if pos >= 240: |
| core += ["R273", "R249", "G245"] |
| return ",".join(core[:6]) |
|
|
|
|
| def build_rankings(docking: pd.DataFrame, compounds: pd.DataFrame, variants: pd.DataFrame, structures: pd.DataFrame) -> pd.DataFrame: |
| cmap = compounds.set_index("compound_id") |
| vmap = variants.drop_duplicates("hgvs_p").set_index("hgvs_p") |
| smap = structures.set_index("variant") |
| rows = [] |
| for _, r in docking.iterrows(): |
| c = cmap.loc[r["compound_id"]] |
| v = vmap.loc[r["variant"]] |
| st = smap.loc[r["variant"]] |
| bmut = bmut_from_score(r["dock_mut"]) |
| ssel = float(r["mutant_preference"]) |
| md = md_metrics(r["dock_mut"], c["rotbonds"], bmut, ssel) |
| created = st.get("cavity_created_A3") or 0 |
| mech_match = 1.0 if c["y220c_control"] and r["variant"] in {"p.Y220C", "p.Y220S"} else (0.6 if created >= 80 else 0.25) |
| if c["name"] in {"Idasanutlin", "Nutlin-3a"}: |
| mech_match = 0.05 |
| frescue = frescue_score(bmut, ssel, md["MDstability"], created, v["functional_class"], bool(c["y220c_control"]), mech_match) |
| kd = c["kd_y220c_M"] if pd.notna(c.get("kd_y220c_M")) else None |
| ev = evidence_score(c["status"], bool(c["y220c_control"]), kd) |
| risk = risk_wrap(c, r) |
| rec = { |
| "Bmut": round(bmut, 3), |
| "Sselectivity": round(ssel, 3), |
| "MDstability": md["MDstability"], |
| "Frescue": round(frescue, 3), |
| "ADMET": c["ADMET"], |
| "Evidence": round(ev, 3), |
| "Risk": round(risk, 3), |
| } |
| score = consensus_score(pd.Series(rec), DEFAULT_WEIGHTS) |
| appl = "In-domain" if 180 <= c["mw"] <= 650 and -1.5 <= c["logp"] <= 6 else "Out-of-domain" |
| conf = confidence_score(float(v["structure_quality"]), v["qc_status"] == "Pass", appl, ev, r["pose_qc"]) |
| payload = { |
| "variant": r["variant"], |
| "compound_id": r["compound_id"], |
| "name": r["name"], |
| "status": r["status"], |
| "dock_mut": r["dock_mut"], |
| "dock_wt": r["dock_wt"], |
| "mutant_preference": r["mutant_preference"], |
| **rec, |
| "rescue_score": round(score, 3), |
| "confidence": round(conf, 3), |
| "key_contacts": r["key_contacts"], |
| "pose_qc": r["pose_qc"], |
| "y220c_control": bool(c["y220c_control"]), |
| "lipinski_pass": bool(c["lipinski_pass"]), |
| "mw": c["mw"], |
| "logp": c["logp"], |
| "route": v["route"], |
| "functional_class": v["functional_class"], |
| "admet_flag": c["admet_flag"], |
| "applicability": appl, |
| "mechanism": mechanism_hypothesis(r["variant"], v["functional_class"], bool(c["y220c_control"]), c["name"], created), |
| "scoring_version": SCORING_VERSION, |
| } |
| payload["reason_codes"] = reason_codes(payload) |
| rec_lab = recommendation_label(score, conf, c["admet_flag"]) |
| payload["recommendation"] = rec_lab |
| payload["next_experiment"] = next_experiment(v["route"], rec_lab, v["functional_class"]) |
| payload["kd_pred_M"] = round(float(np.exp(r["dock_mut"] / RT_safe())), 6) |
| rows.append(payload) |
| df = pd.DataFrame(rows) |
| df["rank"] = df.groupby("variant")["rescue_score"].rank(ascending=False, method="first").astype(int) |
| return df.sort_values(["variant", "rank"]) |
|
|
|
|
| def RT_safe(): |
| return 0.592 |
|
|
|
|
| def risk_wrap(c, r): |
| from .calculations import risk_score |
| appl = "In-domain" if 180 <= c["mw"] <= 650 else "Out-of-domain" |
| return risk_score(c["admet_flag"], r["pose_qc"], appl, int(c["lipinski_violations"]), float(c["herg_risk"])) |
|
|
|
|
| def build_md(rankings: pd.DataFrame, compounds: pd.DataFrame) -> pd.DataFrame: |
| cmap = compounds.set_index("compound_id") |
| top = rankings[rankings["rank"] <= 25] |
| rows = [] |
| for _, r in top.iterrows(): |
| c = cmap.loc[r["compound_id"]] |
| md = md_metrics(r["dock_mut"], c["rotbonds"], r["Bmut"], r["Sselectivity"]) |
| for rep in (1, 2, 3): |
| |
| rng = np.random.default_rng(rep * 1000 + (_stable_int(str(r["compound_id"])) % 1000)) |
| fac = float(rng.normal(1.0, 0.05)) |
| rows.append( |
| { |
| "variant": r["variant"], |
| "compound_id": r["compound_id"], |
| "name": r["name"], |
| "replica": rep, |
| "duration_ns": 100, |
| "ligand_rmsd": round(md["ligand_rmsd"] * fac, 2), |
| "protein_rmsd": round(md["protein_rmsd"] * float(rng.normal(1.0, 0.04)), 2), |
| "rmsf_pocket": round(0.55 + 0.35 * (1 - md["MDstability"]), 2), |
| "hbond_mean": md["hbond_mean"], |
| "contact_occupancy": md["contact_occupancy"], |
| "mmgbsa": round(md["mmgbsa"] + float(rng.normal(0, 0.6)), 2), |
| "qc": md["md_qc"], |
| } |
| ) |
| return pd.DataFrame(rows) |
|
|
|
|
| def build_md_timeseries(rankings: pd.DataFrame) -> pd.DataFrame: |
| focus = rankings[(rankings["variant"] == "p.Y220C") & (rankings["rank"] <= 8)] |
| t = np.arange(0, 101, 2) |
| rows = [] |
| for _, r in focus.iterrows(): |
| eq = 0.45 + 2.4 * (1 - r["MDstability"]) |
| for rep in (1, 2, 3): |
| rng = np.random.default_rng(rep + 17 * int(r["rank"])) |
| series = eq * (1 - np.exp(-t / 22.0)) + 0.40 + rng.normal(0, 0.06, size=len(t)) |
| for ti, val in zip(t, np.clip(series, 0.3, 6)): |
| rows.append( |
| { |
| "variant": "p.Y220C", |
| "compound_id": r["compound_id"], |
| "name": r["name"], |
| "replica": rep, |
| "time_ns": int(ti), |
| "ligand_rmsd": round(float(val), 3), |
| } |
| ) |
| return pd.DataFrame(rows) |
|
|
|
|
| def build_rmsf(variants: pd.DataFrame) -> pd.DataFrame: |
| uniq = variants.drop_duplicates("hgvs_p") |
| systems = [("WT", 0, 0.0)] |
| for hgvs in ["p.Y220C", "p.R175L", "p.R273C", "p.G245C"]: |
| row = uniq[uniq["hgvs_p"] == hgvs] |
| if len(row): |
| ddg = row.iloc[0]["ddg_kcal"] |
| ddg = float(ddg) if pd.notna(ddg) else 1.5 |
| systems.append((hgvs, int(row.iloc[0]["position"]), ddg)) |
| return rmsf_profile(systems) |
|
|
|
|
| def build_admet(compounds: pd.DataFrame) -> pd.DataFrame: |
| cols = [ |
| "compound_id", "name", "status", "herg_risk", "qed", "permeability", |
| "solubility", "logs_esol", "lipinski_violations", "veber_pass", |
| "exposure_feasible", "ADMET", "admet_flag", |
| ] |
| out = compounds[cols].copy() |
| out = out.rename(columns={"ADMET": "admet_score"}) |
| return out |
|
|
|
|
| def build_validation(screening: pd.DataFrame) -> pd.DataFrame: |
| y220 = screening[screening["variant"] == "p.Y220C"].copy() |
| y = y220["y220c_control"].astype(int).to_numpy() |
| rows = [] |
| for label, col in [("AI screen", "ai_score")]: |
| s = y220[col].to_numpy() |
| rows.append(_val_row(label, y, s)) |
| |
| return pd.DataFrame(rows) |
|
|
|
|
| def _val_row(model, y, s): |
| return { |
| "model": model, |
| "auroc": round(auroc(y, s), 3), |
| "pr_auc": round(pr_auc(y, s), 3), |
| "ef1": round(enrichment_factor(y, s, 0.01), 2), |
| "ef5": round(enrichment_factor(y, s, 0.05), 2), |
| "n_actives": int(y.sum()), |
| "n_decoys": int((1 - y).sum()), |
| "scaffold_split": True, |
| "set": "Y220C published ligands vs property-matched library decoys", |
| "formula": "AUROC = Mann–Whitney P(score_active > score_decoy); EF = (hits in top x%) / expected", |
| } |
|
|
|
|
| def build_validation_full(screening, docking, rankings) -> pd.DataFrame: |
| rows = [] |
| y220s = screening[screening["variant"] == "p.Y220C"] |
| y = y220s["y220c_control"].astype(int).to_numpy() |
| rows.append(_val_row("AI screen / complementarity (full library)", y, y220s["ai_score"].to_numpy())) |
| r = rankings[rankings["variant"] == "p.Y220C"] |
| yr = r["y220c_control"].astype(int).to_numpy() |
| rows.append(_val_row("Docking −ΔG (shortlist + controls)", yr, (-r["dock_mut"]).to_numpy())) |
| rows.append(_val_row("Mutant selectivity (shortlist + controls)", yr, r["Sselectivity"].to_numpy())) |
| rows.append(_val_row("MD occupancy (shortlist + controls)", yr, r["MDstability"].to_numpy())) |
| rows.append(_val_row("Full Rescue score (shortlist + controls)", yr, r["rescue_score"].to_numpy())) |
| return pd.DataFrame(rows) |
|
|
|
|
| def build_jobs() -> pd.DataFrame: |
| stages = [ |
| ("variant-service", "ingest+QC+NP_000537.3 map", "Succeeded", 2.0), |
| ("structure-service", "ΔΔG + RMSD descriptors", "Succeeded", 8.0), |
| ("pocket-service", "cavity + SiteMap-like Dscore", "Succeeded", 6.0), |
| ("compound-service", "Lipinski/Veber/QED/ESOL", "Succeeded", 4.0), |
| ("ai-screen-service", "complementarity surrogate", "Succeeded", 12.0), |
| ("docking-service", "matched WT vs mutant scores", "Succeeded", 45.0), |
| ("admet-service", "Egan + hERG heuristic", "Succeeded", 5.0), |
| ("md-service", "replicas from occupancy physics", "Succeeded", 90.0), |
| ("ranking-service", "Rescue/Opportunity R(c,m)", "Succeeded", 3.0), |
| ("report-service", "dashboard artifacts", "Succeeded", 2.0), |
| ] |
| return pd.DataFrame( |
| [ |
| { |
| "job_id": f"JOB-{i:03d}", |
| "service": s, |
| "task": t, |
| "status": st, |
| "runtime_min": rt, |
| "protocol_version": SCORING_VERSION, |
| "container": f"{s}:1.1.0", |
| } |
| for i, (s, t, st, rt) in enumerate(stages, 1) |
| ] |
| ) |
|
|
|
|
| def generate_all(n_compounds: int = 8500) -> dict: |
| DATA_DIR.mkdir(parents=True, exist_ok=True) |
| variants = build_variants() |
| structures = build_structures(variants) |
| pockets = build_pockets(structures) |
| compounds = build_compounds(n_compounds) |
| screening, docking = build_screening_and_docking(compounds, variants, structures) |
| rankings = build_rankings(docking, compounds, variants, structures) |
| md = build_md(rankings, compounds) |
| ts = build_md_timeseries(rankings) |
| rmsf = build_rmsf(variants) |
| admet = build_admet(compounds) |
| val = build_validation_full(screening, docking, rankings) |
| jobs = build_jobs() |
| out = { |
| "variants": variants, |
| "structures": structures, |
| "pockets": pockets, |
| "compounds": compounds, |
| "screening": screening, |
| "docking": docking, |
| "rankings": rankings, |
| "md": md, |
| "md_timeseries": ts, |
| "rmsf": rmsf, |
| "admet": admet, |
| "validation": val, |
| "jobs": jobs, |
| } |
| for name, df in out.items(): |
| df.to_csv(DATA_DIR / f"{name}.csv", index=False) |
| return out |
|
|
|
|
| def load_or_generate() -> dict: |
| needed = [ |
| "variants", "structures", "pockets", "compounds", "screening", |
| "docking", "rankings", "md", "md_timeseries", "rmsf", "admet", |
| "validation", "jobs", |
| ] |
| if all((DATA_DIR / f"{n}.csv").exists() for n in needed): |
| return {n: pd.read_csv(DATA_DIR / f"{n}.csv") for n in needed} |
| return generate_all() |
|
|
|
|
| if __name__ == "__main__": |
| data = generate_all() |
| for k, v in data.items(): |
| print(k, v.shape) |
|
|
