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

In this tutorial, we build an end-to-end autonomous AI co-scientist workflow for next-generation EGFR inhibitor discovery, focusing on the C797S osimertinib-resistance mutation in non-small cell lung cancer. We start by resolving the biological target through ChEMBL and UniProt, then mine curated EGFR IC50 bioactivity records and convert them into a clean pIC50 modeling dataset. We use RDKit to standardize molecules, remove salts, aggregate replicate measurements, compute Morgan fingerprints, extract physicochemical descriptors, and analyze scaffold diversity so that our model learns from chemically meaningful representations rather than raw strings. From there, we train a scaffold-split Random Forest QSAR model, evaluate its ability to generalize to unseen chemotypes, interpret potency-driving features with SHAP or model importances, and visualize influential molecular substructures. Finally, we move beyond prediction into generative design by recombining BRICS fragments from potent drug-like actives, scoring the resulting virtual analogs across potency, drug-likeness, synthesizability, novelty, and developability gates, and cross-checking the shortlisted candidates against PubChem.

EGFR Target Setup

Copy CodeCopiedUse a different Browser

import sys, subprocess, importlib, warnings, time, os, random, json warnings.filterwarnings("ignore") def _pip(*pkgs): subprocess.run([sys.executable, "-m", "pip", "install", "-q", *pkgs], check=False) for mod, pkg in [("rdkit", "rdkit"), ("shap", "shap"), ("requests", "requests")]: try: importlib.import_module(mod) except Exception: print(f"Installing {pkg} ...") _pip(pkg) import numpy as np import pandas as pd import requests import matplotlib.pyplot as plt from scipy.stats import spearmanr from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import r2_score, mean_squared_error, roc_auc_score from sklearn.decomposition import PCA from rdkit import Chem, DataStructs, RDLogger from rdkit.Chem import Descriptors, Draw, QED, rdMolDescriptors, BRICS, rdFingerprintGenerator from rdkit.Chem.Scaffolds import MurckoScaffold RDLogger.DisableLog("rdApp.*") try: from rdkit.Chem.MolStandardize import rdMolStandardize _HAS_STD = True except Exception: _HAS_STD = False try: from rdkit.Chem import RDConfig sys.path.append(os.path.join(RDConfig.RDContribDir, "SA_Score")) import sascorer _HAS_SA = True except Exception: _HAS_SA = False TARGET_CHEMBL_ID = "CHEMBL203" TARGET_QUERY = "Epidermal growth factor receptor" FALLBACK_CHEMBL_ID = "CHEMBL203" NBITS, RADIUS = 2048, 2 RANDOM_STATE = 42 MAX_ACTIVITIES = 9000 MAX_UNIQUE = 4000 ACTIVE_PIC50 = 7.0 BRICS_MAX_TRIES = 4000 N_FRAG_PARENTS = 60 N_SHORTLIST = 12 np.random.seed(RANDOM_STATE); random.seed(RANDOM_STATE) BASE = "https://www.ebi.ac.uk/chembl/api/data" HDRS = {"Accept": "application/json", "User-Agent": "ai-coscientist-tutorial/1.0"} def banner(title): print("\n" + "=" * 86 + f"\n {title}\n" + "=" * 86) def http_json(url, params=None, tries=3, timeout=45): for k in range(tries): try: r = requests.get(url, params=params, headers=HDRS, timeout=timeout) if r.status_code == 200: return r.json() if r.status_code == 404: return None except Exception: pass time.sleep(1.5 * (k + 1)) return None banner("[1/9] TARGET INTELLIGENCE (ChEMBL + UniProt)") print("Question: What target are we drugging, and why is it hard?\n") def ic50_count(tid): js = http_json(f"{BASE}/activity", {"target_chembl_id": tid, "standard_type": "IC50", "pchembl_value__isnull": "false", "limit": 1, "format": "json"}) try: return int(js["page_meta"]["total_count"]) except Exception: return 0 target_id, target_name, uniprot_acc = None, TARGET_QUERY, None if TARGET_CHEMBL_ID: target_id = TARGET_CHEMBL_ID else: srch = http_json(f"{BASE}/target/search", {"q": TARGET_QUERY, "format": "json"}) cands = [] if srch: for t in srch.get("targets", []): if t.get("organism") == "Homo sapiens" and t.get("target_type") == "SINGLE PROTEIN": cands.append(t) cands = sorted(cands, key=lambda t: float(t.get("score", 0)), reverse=True)[:8] scored = [(t, ic50_count(t["target_chembl_id"])) for t in cands] scored = [(t, n) for t, n in scored if n > 0] if scored: best = max(scored, key=lambda x: x[1])[0] target_id, target_name = best["target_chembl_id"], best.get("pref_name", TARGET_QUERY) print(f" Auto-resolved '{TARGET_QUERY}' by data volume -> {target_id}") else: target_id = FALLBACK_CHEMBL_ID print(f" Auto-resolve found no data; falling back to {FALLBACK_CHEMBL_ID}") det = http_json(f"{BASE}/target/{target_id}", {"format": "json"}) if det and det.get("pref_name"): target_name = det["pref_name"] if det: for comp in det.get("target_components", []): if comp.get("accession"): uniprot_acc = comp["accession"]; break print(f" Resolved target : {target_name}") print(f" ChEMBL ID : {target_id}") print(f" UniProt : {uniprot_acc}") if uniprot_acc: uni = http_json(f"https://rest.uniprot.org/uniprotkb/{uniprot_acc}.json") if uni: try: fn = next(c["texts"][0]["value"] for c in uni.get("comments", []) if c.get("commentType") == "FUNCTION") print("\n Function (UniProt):") print(" ", (fn[:340] + " ...") if len(fn) > 340 else fn) except Exception: pass print(""" Resistance context: 1st/2nd/3rd-gen EGFR TKIs lose potency once tumours acquire the C797S mutation, which abolishes the covalent cysteine anchor exploited by osimertinib. Goal of this run: learn the chemistry of known EGFR inhibitors and propose NOVEL, drug-like analogs as starting points for a C797S-active 4th-generation series.""")

We begin by preparing the full scientific computing environment and installing any missing chemistry, modeling, plotting, and API dependencies required by the workflow. We configure the EGFR target settings, define modeling constants, initialize reproducible random seeds, and create helper functions for banners and robust JSON API calls. We then resolve the ChEMBL target, retrieve UniProt context when available, and frame the biological motivation around EGFR C797S resistance.

Mining ChEMBL Bioactivity Data

Copy CodeCopiedUse a different Browser

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 = 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() 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= ACTIVE_PIC50).astype(int) auc = roc_auc_score(ycls, pred) if len(np.unique(ycls)) == 2 else float("nan") print(f" Held-out (new-scaffold) performance on {len(te)} molecules:") print(f" R^2 = {r2:.3f}") print(f" RMSE (pIC50) = {rmse:.3f} (~{rmse:.2f} log units)") print(f" Spearman rho = {rho:.3f}") print(f" ROC-AUC active = {auc:.3f} (ranking potent vs weak)") model_full = RandomForestRegressor(n_estimators=400, max_features="sqrt", n_jobs=-1, random_state=RANDOM_STATE).fit(X, y)

We analyze the curated chemical space by extracting Murcko scaffolds and identifying the most common chemotype families in the dataset. We project Morgan fingerprints with PCA to visualize how EGFR inhibitors distribute across chemical space and how potency varies across that landscape. We then train a Random Forest QSAR model using a scaffold split, which helps us evaluate whether the model generalizes to unseen molecular scaffolds rather than memorizing close analogs.

Interpretability and BRICS Generation

Copy CodeCopiedUse a different Browser

banner("[6/9] MODEL INTERPRETABILITY (which substructures drive potency?)") top_feat_idx, shap_ok = None, False try: import shap samp = np.random.RandomState(RANDOM_STATE).choice(len(te), min(300, len(te)), replace=False) expl = shap.TreeExplainer(model) sv = expl.shap_values(X[te][samp]) if isinstance(sv, list): sv = sv[0] imp = np.abs(sv).mean(0) shap_ok = True print(" Using SHAP TreeExplainer (mean |SHAP| over held-out molecules).") except Exception as e: imp = model_full.feature_importances_ print(f" SHAP unavailable ({type(e).name}); using RandomForest importances.") order = np.argsort(imp)[::-1] top_feat_idx = [i for i in order[:25]] top_desc = [(FEAT_NAMES[i], imp[i]) for i in order if i >= NBITS][:8] top_bits = [i for i in order if i 0: env = Chem.FindAtomEnvironmentOfRadiusN(m, rad, atom) bonds = list(env) for bidx in env: b = m.GetBondWithIdx(bidx) atoms.update((b.GetBeginAtomIdx(), b.GetEndAtomIdx())) try: return Draw.MolToImage(m, size=(300, 240), highlightAtoms=list(atoms), highlightBonds=bonds) except Exception: return None return None try: imgs = [(b, bit_exemplar(b)) for b in top_bits] imgs = [(b, im) for b, im in imgs if im is not None] if imgs: fig, ax = plt.subplots(1, len(imgs), figsize=(3.1 * len(imgs), 3.3)) if len(imgs) == 1: ax = [ax] for a, (b, im) in zip(ax, imgs): a.imshow(im); a.axis("off"); a.set_title(f"ECFP bit {b}\n(rank imp.)", fontsize=9) plt.suptitle("Substructures the model associates with potency", y=1.02) plt.tight_layout(); plt.savefig("fig2_potency_substructures.png", dpi=120, bbox_inches="tight") plt.show() except Exception as e: print(f" (substructure drawing skipped: {type(e).name})") banner("[7/9] GENERATIVE DESIGN (BRICS fragment recombination -> novel analogs)") seed = data[data["pIC50"] >= ACTIVE_PIC50].copy() seed["mw"] = seed["mol"].map(Descriptors.MolWt) seed = seed[(seed.mw >= 250) & (seed.mw = BRICS_MAX_TRIES: break try: prod.UpdatePropertyCache(strict=False) Chem.SanitizeMol(prod) cs = Chem.MolToSmiles(prod) except Exception: continue if cs in known or cs in gen: continue mw = Descriptors.MolWt(prod) if 250 hi: return float(np.clip((hh - x) / (hh - hi + 1e-9), 0, 1)) return 1.0 rows = [] for smi, m, pp in zip(gsmiles, gmols, gpred): mw, lp = Descriptors.MolWt(m), Descriptors.MolLogP(m) hbd, hba = Descriptors.NumHDonors(m), Descriptors.NumHAcceptors(m) tpsa, rotb = Descriptors.TPSA(m), Descriptors.NumRotatableBonds(m) try: qed = QED.qed(m) except Ex

[truncated for AI cost control]