Building a Scaffold-Split Random Forest QSAR Co-Scientist for EGFR Inhibitor Discovery Using ChEMBL, RDKit, SHAP, and BRICS

banner("[2/9] BIOACTIVITY MINING  (ChEMBL activities -> pIC50)")
def pull_activities(tid, cap):
   url, rows = f"{BASE}/activity", []
   params = {"target_chembl_id": tid, "standard_type": "IC50",
             "pchembl_value__isnull": "false", "limit": 1000, "format": "json"}
   js = http_json(url, params)
   pages = 0
   while js and pages < 60:
       rows.extend(js.get("activities", []))
       pages += 1
       if len(rows) >= cap:
           break
       nxt = js.get("page_meta", {}).get("next")
       if not nxt:
           break
       nurl = nxt if nxt.startswith("http") else "https://www.ebi.ac.uk" + nxt
       js = http_json(nurl)
   return rows[:cap]
raw = pull_activities(target_id, MAX_ACTIVITIES)
print(f"  Pulled {len(raw)} raw IC50 records with a curated pChEMBL value.")
recs = []
for a in raw:
   smi, pv = a.get("canonical_smiles"), a.get("pchembl_value")
   if not smi or pv in (None, ""):
       continue
   if a.get("standard_relation") != "=":
       continue
   if a.get("standard_units") not in ("nM", None):
       continue
   try:
       pv = float(pv)
   except Exception:
       continue
   recs.append({"chembl_id": a.get("molecule_chembl_id"), "smiles": smi, "pIC50": pv})
raw_df = pd.DataFrame(recs, columns=["chembl_id", "smiles", "pIC50"])
print(f"  After quality filters: {len(raw_df)} measurements.")
if len(raw_df) == 0:
   print("\n  STOP: no usable IC50 data was retrieved for this target.\n"
         "  Fix: set TARGET_CHEMBL_ID to a target that has inhibitor data\n"
         "       (e.g. CHEMBL203=EGFR, CHEMBL5251=BTK, CHEMBL2971=JAK2),\n"
         "       or set TARGET_CHEMBL_ID=\"\" to auto-resolve TARGET_QUERY by data volume.")
   raise SystemExit("No bioactivity data for the selected target.")
banner("[3/9] MOLECULAR CURATION  (standardize, de-salt, aggregate)")
_lfc = rdMolStandardize.LargestFragmentChooser() if _HAS_STD else None
_unc = rdMolStandardize.Uncharger() if _HAS_STD else None
def standardize(smi):
   m = Chem.MolFromSmiles(smi)
   if m is None:
       return None, None
   try:
       if _HAS_STD:
           m = _lfc.choose(m); m = _unc.uncharge(m)
       else:
           frags = Chem.GetMolFrags(m, asMols=True, sanitizeFrags=True)
           if frags:
               m = max(frags, key=lambda x: x.GetNumHeavyAtoms())
       return m, Chem.MolToSmiles(m)
   except Exception:
       return None, None
canon, keep_mol = [], {}
for _, r in raw_df.iterrows():
   m, cs = standardize(r["smiles"])
   if cs is None or m.GetNumHeavyAtoms() < 6:
       continue
   canon.append({"smiles": cs, "pIC50": r["pIC50"], "chembl_id": r["chembl_id"]})
   keep_mol[cs] = m
cdf = pd.DataFrame(canon, columns=["smiles", "pIC50", "chembl_id"])
data = (cdf.groupby("smiles")
       .agg(pIC50=("pIC50", "median"), n=("pIC50", "size"),
            chembl_id=("chembl_id", "first")).reset_index())
if len(data) > MAX_UNIQUE:
   data = data.sample(MAX_UNIQUE, random_state=RANDOM_STATE).reset_index(drop=True)
data["mol"] = data["smiles"].map(keep_mol)
n_active = int((data["pIC50"] >= ACTIVE_PIC50).sum())
print(f"  Unique curated molecules : {len(data)}")
print(f"  Potent actives (IC50<=100nM): {n_active}  ({100*n_active/len(data):.1f}%)")
print(f"  pIC50 range: {data.pIC50.min():.2f} - {data.pIC50.max():.2f}  "
     f"(median {data.pIC50.median():.2f})")
mfpgen = rdFingerprintGenerator.GetMorganGenerator(radius=RADIUS, fpSize=NBITS)
DESC = [("MolWt", Descriptors.MolWt), ("MolLogP", Descriptors.MolLogP),
       ("TPSA", Descriptors.TPSA), ("HBD", Descriptors.NumHDonors),
       ("HBA", Descriptors.NumHAcceptors), ("RotB", Descriptors.NumRotatableBonds),
       ("AromRings", Descriptors.NumAromaticRings), ("FracCSP3", Descriptors.FractionCSP3),
       ("HeavyAtoms", Descriptors.HeavyAtomCount),
       ("NumRings", lambda m: rdMolDescriptors.CalcNumRings(m))]
FEAT_NAMES = [f"bit_{i}" for i in range(NBITS)] + [n for n, _ in DESC]
def fp_array(m):
   a = np.zeros((NBITS,), dtype=np.int8)
   DataStructs.ConvertToNumpyArray(mfpgen.GetFingerprint(m), a)
   return a
def featurize(mols):
   Xb = np.zeros((len(mols), NBITS), dtype=np.int8)
   Xd = np.zeros((len(mols), len(DESC)), dtype=np.float32)
   for i, m in enumerate(mols):
       Xb[i] = fp_array(m)
       for j, (_, fn) in enumerate(DESC):
           try:
               Xd[i, j] = fn(m)
           except Exception:
               Xd[i, j] = 0.0
   return np.nan_to_num(np.hstack([Xb, Xd]).astype(np.float32))
X = featurize(list(data["mol"]))
y = data["pIC50"].values
print(f"  Feature matrix: {X.shape[0]} molecules x {X.shape[1]} features "
     f"({NBITS} ECFP bits + {len(DESC)} descriptors)")