fetch_benchmark_ligands.py

"""
fetch_benchmark_ligands.py — build a REAL, REPRODUCIBLE 150-ligand screen set.

Pulls genuine BRAF (CHEMBL5145) inhibitors from ChEMBL, keeps the most potent
unique molecules, maps each to its PubChem CID via InChIKey, and writes a
`compound_list.txt` ready for the AutoDock pipeline.

Why this instead of a hand-typed list of IDs:
  * Reproducible — the same query returns the same molecules for anyone.
  * Real — every CID is a measured BRAF ligand, not a made-up number.
  * Meaningful — the set is dominated by known BRAF inhibitors, so a correct
    pipeline should rank the potent/approved ones (vemurafenib, dabrafenib,
    encorafenib) near the top. Rediscovering known drugs from a blind screen
    is exactly the "it works on real data" signal you want.

Requires (already installed in your envs): chembl_webresource_client, pubchempy.
Run it in an env that has both (e.g. `base`). Needs internet. Takes a few
minutes (it queries PubChem once per molecule).

Usage:
    cd "C:\\DRUG DESIGN PROJECTS\\AutoDock_Pipeline"
    python fetch_benchmark_ligands.py --n 150
    # -> writes data/ligands/compound_list.txt
"""

from __future__ import annotations

import argparse
import time
from pathlib import Path

from chembl_webresource_client.new_client import new_client
import pubchempy as pcp


TARGET_CHEMBL_ID = "CHEMBL5145"   # canonical human BRAF (matches the QSAR target)
ACTIVITY_TYPE = "IC50"


def most_potent_unique(target: str, activity_type: str) -> list[tuple[str, float]]:
    """Return [(molecule_chembl_id, best_pchembl)] sorted most-potent-first."""
    activity = new_client.activity
    rows = activity.filter(
        target_chembl_id=target,
        standard_type=activity_type,
        pchembl_value__isnull=False,
    )
    best: dict[str, float] = {}
    for r in rows:
        cid = r.get("molecule_chembl_id")
        pv = r.get("pchembl_value")
        if not cid or pv is None:
            continue
        pv = float(pv)
        if cid not in best or pv > best[cid]:
            best[cid] = pv
    return sorted(best.items(), key=lambda kv: kv[1], reverse=True)


def chembl_to_pubchem_cid(chembl_id: str) -> int | None:
    """Map a ChEMBL molecule to a PubChem CID via its standard InChIKey."""
    molecule = new_client.molecule
    record = molecule.get(chembl_id)
    structures = (record or {}).get("molecule_structures") or {}
    inchikey = structures.get("standard_inchi_key")
    if not inchikey:
        return None
    hits = pcp.get_compounds(inchikey, "inchikey")
    if hits and hits[0].cid:
        return int(hits[0].cid)
    return None


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--n", type=int, default=150, help="How many ligands.")
    parser.add_argument(
        "--out",
        type=Path,
        default=Path("data/ligands/compound_list.txt"),
        help="Output compound list (PubChem CIDs, one per line).",
    )
    args = parser.parse_args()

    print(f"Querying ChEMBL for {ACTIVITY_TYPE} actives of {TARGET_CHEMBL_ID} (BRAF)...")
    ranked = most_potent_unique(TARGET_CHEMBL_ID, ACTIVITY_TYPE)
    print(f"  {len(ranked)} unique molecules with a pChEMBL value.")

    top = ranked[: args.n]
    print(f"Mapping the top {len(top)} to PubChem CIDs (one lookup each)...")

    resolved: list[tuple[int, str, float]] = []
    for i, (chembl_id, pv) in enumerate(top, 1):
        try:
            cid = chembl_to_pubchem_cid(chembl_id)
        except Exception as exc:  # network / lookup hiccup: skip, keep going
            cid = None
            print(f"  [{i}/{len(top)}] {chembl_id}: lookup error ({exc})")
        if cid:
            resolved.append((cid, chembl_id, pv))
            print(f"  [{i}/{len(top)}] {chembl_id}  pChEMBL={pv:.2f}  -> CID {cid}")
        else:
            print(f"  [{i}/{len(top)}] {chembl_id}: no PubChem CID (skipped)")
        time.sleep(0.2)  # be polite to PubChem

    args.out.parent.mkdir(parents=True, exist_ok=True)
    with args.out.open("w", encoding="utf-8") as fh:
        fh.write("# BRAF (CHEMBL5145) benchmark ligands - real, reproducible.\n")
        fh.write("# Built by fetch_benchmark_ligands.py from ChEMBL IC50 actives.\n")
        fh.write(f"# {len(resolved)} PubChem CIDs, most-potent-first.\n")
        for cid, chembl_id, pv in resolved:
            fh.write(f"{cid}\n")

    print(f"\nWrote {len(resolved)} PubChem CIDs to {args.out}")
    print("Now run the pipeline:  python run_pipeline.py")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())