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")]:

importlib.import_module(mod)

except Exception:

print(f"Installing {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.*")

from rdkit.Chem.MolStandardize import rdMolStandardize

_HAS_STD = True

except Exception:

_HAS_STD = False

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):

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:

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"})

return int(js["page_meta"]["total_count"])

except Exception:

target_id, target_name, uniprot_acc = None, TARGET_QUERY, None

if TARGET_CHEMBL_ID:

target_id = TARGET_CHEMBL_ID

srch = http_json(f"{BASE}/target/search", {"q": TARGET_QUERY, "format": "json"})

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]

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}")

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"]

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")

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:

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)

while js and pages < 60:

rows.extend(js.get("activities", []))

if len(rows) >= cap:

nxt = js.get("page_meta", {}).get("next")

if not nxt:

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.")

for a in raw:

smi, pv = a.get("canonical_smiles"), a.get("pchembl_value")

if not smi or pv in (None, ""):

if a.get("standard_relation") != "=":

if a.get("standard_units") not in ("nM", None):

pv = float(pv)

except Exception:

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

if _HAS_STD:

m = _lfc.choose(m); m = _unc.uncharge(m)

frags = Chem.GetMolFrags(m, asMols=True, sanitizeFrags=True)

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:

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)

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):

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)")

We mine curated IC50 bioactivity measurements from ChEMBL and convert the raw activity records into a usable pIC50 dataset. We filter out incomplete, non-exact, or inconsistent measurements so that the downstream QSAR model trains on cleaner potency values. We then standardize the molecules with RDKit, remove salts or smaller fragments, aggregate duplicate molecules by median pIC50, and convert each molecule into Morgan fingerprint bits plus interpretable physicochemical descriptors.

Scaffold Analysis and QSAR

Copy CodeCopiedUse a different Browser

banner("[4/9] CHEMICAL SPACE & SCAFFOLD ANALYSIS")

def murcko(m):

s = MurckoScaffold.GetScaffoldForMol(m)

cs = Chem.MolToSmiles(s)

return cs if cs else "(acyclic)"

except Exception:

return "(error)"

data["scaffold"] = data["mol"].map(murcko)

top_scaf = data["scaffold"].value_counts().head(10)

print(" Top recurring Murcko scaffolds (chemotype families):")

for i, (sc, c) in enumerate(top_scaf.items(), 1):

print(f" {i:2d}. n={c:4d} {sc[:70]}")

pca = PCA(n_components=2, random_state=RANDOM_STATE)

emb = pca.fit_transform(X[:, :NBITS])

fig, ax = plt.subplots(1, 3, figsize=(18, 4.8))

sc0 = ax[0].scatter(emb[:, 0], emb[:, 1], c=y, cmap="viridis", s=10, alpha=0.6)

ax[0].set(title="Chemical space (ECFP-PCA), coloured by potency",

xlabel=f"PC1 ({pca.explained_variance_ratio_[0]*100:.0f}%)",

ylabel=f"PC2 ({pca.explained_variance_ratio_[1]*100:.0f}%)")

plt.colorbar(sc0, ax=ax[0], label="pIC50")

ax[1].hist(y, bins=40, color="#3b7dd8", edgecolor="white")

ax[1].axvline(ACTIVE_PIC50, color="crimson", ls="--", label=f"active cut (pIC50={ACTIVE_PIC50})")

ax[1].set(title="Potency distribution", xlabel="pIC50", ylabel="molecules"); ax[1].legend()

ax[2].barh([s[:22] + "..." for s in top_scaf.index[::-1]], top_scaf.values[::-1], color="#6c5ce7")

ax[2].set(title="Top 10 scaffolds", xlabel="count")

plt.tight_layout(); plt.savefig("fig1_chemical_space.png", dpi=120); plt.show()

banner("[5/9] INTERPRETABLE QSAR MODEL (scaffold-split, leakage-free)")

def scaffold_split(scaffolds, test_frac=0.2, seed=RANDOM_STATE):

groups = {}

for i, s in enumerate(scaffolds):

groups.setdefault(s, []).append(i)

order = list(groups.values())

random.Random(seed).shuffle(order)

n_test = int(len(scaffolds) * test_frac)

test, train = [], []

for g in order:

(test if len(test) < n_test else train).extend(g)

return np.array(sorted(train)), np.array(sorted(test))

tr, te = scaffold_split(data["scaffold"].values, 0.2)

model = RandomForestRegressor(n_estimators=400, max_features="sqrt",

n_jobs=-1, random_state=RANDOM_STATE)

model.fit(X[tr], y[tr])

pred = model.predict(X[te])

r2 = r2_score(y[te], pred)

rmse = mean_squared_error(y[te], pred) ** 0.5

rho = spearmanr(y[te], pred).statistic

ycls = (y[te] >= 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

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):

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 < NBITS][:6]

print("\n Most influential physicochemical descriptors:")

for name, v in top_desc:

print(f" {name:12s} importance={v:.4f}")

train_smis = list(data["smiles"].iloc[tr])

def bit_exemplar(bit):

for smi in train_smis:

m = Chem.MolFromSmiles(smi)

if m is None:

ao = rdFingerprintGenerator.AdditionalOutput(); ao.AllocateBitInfoMap()

_ = mfpgen.GetFingerprint(m, additionalOutput=ao)

bi = ao.GetBitInfoMap()

if bit in bi and len(bi[bit]):

atom, rad = bi[bit][0]

atoms, bonds = {atom}, []

if rad > 0:

env = Chem.FindAtomEnvironmentOfRadiusN(m, rad, atom)

bonds = list(env)

for bidx in env:

b = m.GetBondWithIdx(bidx)

atoms.update((b.GetBeginAtomIdx(), b.GetEndAtomIdx()))

return Draw.MolToImage(m, size=(300, 240),

highlightAtoms=list(atoms), highlightBonds=bonds)

except Exception:

return None

return None

imgs = [(b, bit_exemplar(b)) for b in top_bits]

imgs = [(b, im) for b, im in imgs if im is not None]

fig, ax = plt.subplots(1, len(imgs), figsize=(3.1 * len(imgs), 3.3))

if len(imgs) == 1:

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")

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 <= 500)].sort_values("pIC50", ascending=False).head(N_FRAG_PARENTS)

print(f" Seeding generative design with {len(seed)} potent, drug-like parent molecules.")

frags = set()

for m in seed["mol"]:

frags.update(BRICS.BRICSDecompose(m))

except Exception:

frag_mols = [f for f in (Chem.MolFromSmiles(s) for s in frags) if f is not None]

print(f" Fragment pool: {len(frag_mols)} BRICS fragments.")

known = set(data["smiles"])

for i, prod in enumerate(BRICS.BRICSBuild(frag_mols, scrambleReagents=True, maxDepth=2)):

if i >= BRICS_MAX_TRIES:

prod.UpdatePropertyCache(strict=False)

Chem.SanitizeMol(prod)

cs = Chem.MolToSmiles(prod)

except Exception:

if cs in known or cs in gen:

mw = Descriptors.MolWt(prod)

if 250 <= mw <= 600 and 8 <= prod.GetNumHeavyAtoms() <= 45:

gen[cs] = prod

except Exception as e:

print(f" (BRICS build ended early: {type(e).__name__})")

print(f" Generated {len(gen)} unique, novel, size-reasonable virtual molecules.")

We interpret the trained QSAR model by estimating which descriptors and fingerprint bits contribute most strongly to predicted potency. We use SHAP when available; otherwise, we fall back to Random Forest feature importances to keep the workflow robust. We also visualize representative molecular substructures associated with influential ECFP bits, then begin generative design by decomposing potent drug-like parent molecules into BRICS fragments and recombining them to generate novel virtual analogs.

Multi-Parameter Candidate Scoring

Copy CodeCopiedUse a different Browser

banner("[8/9] MULTI-PARAMETER PRIORITISATION")

gsmiles = list(gen.keys())

gmols = [gen[s] for s in gsmiles]

gX = featurize(gmols)

gpred = model_full.predict(gX)

train_fps = [mfpgen.GetFingerprint(m) for m in data["mol"]]

def novelty(m):

sims = DataStructs.BulkTanimotoSimilarity(mfpgen.GetFingerprint(m), train_fps)

return 1.0 - (max(sims) if sims else 0.0)

def desirability(x, lo, hi, hard_lo=None, hard_hi=None):

hl = hard_lo if hard_lo is not None else lo

hh = hard_hi if hard_hi is not None else hi

return float(np.clip((x - hl) / (lo - hl + 1e-9), 0, 1))

return float(np.clip((hh - x) / (hh - hi + 1e-9), 0, 1))

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)

qed = QED.qed(m)

except Exception:

qed = np.nan

sa = sascorer.calculateScore(m) if _HAS_SA else np.nan

lip = int(mw <= 500) + int(lp <= 5) + int(hbd <= 5) + int(hba <= 10)

veber = (rotb <= 10) and (tpsa <= 140)

nov = novelty(m)

d_pot = desirability(pp, 7.5, 12, hard_lo=5.5)

d_mw = desirability(mw, 250, 500, hard_lo=150, hard_hi=650)

d_lp = desirability(lp, 1, 4, hard_lo=-1, hard_hi=6)

d_sa = desirability(-(sa if not np.isnan(sa) else 3), -3.5, -1, hard_lo=-6)

score = (0.40 * d_pot + 0.20 * (qed if not np.isnan(qed) else 0.5) +

0.10 * d_mw + 0.10 * d_lp + 0.10 * d_sa + 0.10 * nov)

rows.append(dict(smiles=smi, pred_pIC50=pp, MolWt=mw, MolLogP=lp, TPSA=tpsa,

HBD=hbd, HBA=hba, QED=qed, SA=sa, novelty=nov,

lipinski=lip, veber_ok=veber, score=score))

_CANDCOLS = ["smiles", "pred_pIC50", "MolWt", "MolLogP", "TPSA", "HBD", "HBA",

"QED", "SA", "novelty", "lipinski", "veber_ok", "score"]

cand = pd.DataFrame(rows, columns=_CANDCOLS)

gate = cand[(cand.pred_pIC50 >= 6.5) & (cand.MolWt.between(250, 600)) &

(cand.lipinski >= 3) & (cand.veber_ok) &

((cand.SA <= 6) | cand.SA.isna()) & (cand.novelty >= 0.35)]

gate = gate.sort_values("score", ascending=False).reset_index(drop=True)

print(f" {len(gate)} of {len(cand)} generated molecules passed the developability gate.")

print(f" (gate: predicted pIC50>=6.5, MW 250-600, <=1 Lipinski violation, Veber OK,")

print(f" SA<=6, novelty>=0.35 vs all known EGFR inhibitors)\n")

shortlist = gate.head(N_SHORTLIST).copy()

We score the generated molecules with the full QSAR model and evaluate them across potency, molecular weight, lipophilicity, hydrogen bonding, polar surface area, rotatable bonds, QED, synthetic accessibility, and novelty. We calculate a multi-parameter desirability score that balances predicted potency with drug-likeness, synthesizability, and structural distance from known EGFR inhibitors. We then apply hard developability gates and keep only the strongest candidates for the final shortlist.

PubChem Novelty Cross-Check

Copy CodeCopiedUse a different Browser

banner("[9/9] NOVELTY CROSS-CHECK (PubChem) & FINAL SHORTLIST")

def pubchem_cid(smi):

m = Chem.MolFromSmiles(smi)

ik = Chem.MolToInchiKey(m)

except Exception:

return "inchikey_error"

return "inchikey_error"

js = http_json(f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ik}/cids/JSON")

time.sleep(0.25)

return f"CID {js['IdentifierList']['CID'][0]}"

except Exception:

return "NOT in PubChem (putatively new)"

shortlist["pubchem"] = shortlist["smiles"].map(pubchem_cid)

show_cols = ["pred_pIC50", "MolWt", "MolLogP", "QED", "SA", "novelty",

"lipinski", "score", "pubchem"]

pretty = shortlist[show_cols].copy()

pretty.insert(0, "rank", range(1, len(pretty) + 1))

pd.set_option("display.width", 200, "display.max_colwidth", 40)

print("\nTOP CANDIDATE 4th-GENERATION EGFR-INHIBITOR STARTING POINTS:\n")

print(pretty.round(3).to_string(index=False))

mols = [Chem.MolFromSmiles(s) for s in shortlist["smiles"]]

legs = [f"#{i+1} pIC50~{r.pred_pIC50:.1f} | QED {r.QED:.2f} | nov {r.novelty:.2f}"

for i, r in shortlist.reset_index().iterrows()]

grid = Draw.MolsToGridImage(mols, molsPerRow=4, subImgSize=(300, 250), legends=legs)

grid.save("fig3_top_candidates.png")

from IPython.display import display

display(grid)

except Exception:

plt.figure(figsize=(14, 9)); plt.imshow(grid); plt.axis("off"); plt.show()

except Exception as e:

print(f" (candidate drawing skipped: {type(e).__name__})")

shortlist.to_csv("egfr_coscientist_candidates.csv", index=False)

banner("AUTONOMOUS RESEARCH SUMMARY")

n_new = int((shortlist["pubchem"].str.startswith("NOT")).sum())

print(f""" Target : {target_name} ({target_id}) -- overcoming C797S resistance

Evidence base : {len(data)} curated, de-duplicated EGFR inhibitors from ChEMBL

Learned model : RandomForest QSAR, scaffold-split R^2={r2:.2f}, ROC-AUC={auc:.2f}

-> generalises to unseen chemotypes, not memorising analogs

Key drivers : {", ".join(n for n, _ in top_desc[:4])} + specific ECFP substructures

Invented : {len(gen)} novel virtual analogs via BRICS fragment recombination

Prioritised : {len(gate)} passed developability gates; top {len(shortlist)} shortlisted

Novelty audit : {n_new}/{len(shortlist)} shortlisted molecules are absent from PubChem

Artifacts written to disk:

- egfr_coscientist_candidates.csv (full scored shortlist)

- fig1_chemical_space.png (chemical space + potency + scaffolds)

- fig2_potency_substructures.png (SHAP-implicated substructures)

- fig3_top_candidates.png (structures of the shortlist)

NEXT STEPS a wet-lab team would take: dock the shortlist into the EGFR(L858R/T790M/C797S)

triple-mutant structure, prioritise allosteric binders, check synthetic routes, and

assay the top ~5 for C797S potency and selectivity vs wild-type EGFR.

Reminder: this is an educational in-silico hypothesis generator, not a validated drug

pipeline. Predictions require experimental confirmation.""")

print("\nDONE.")

We cross-check the shortlisted molecules against PubChem using an InChIKey lookup to determine whether each candidate is known or putatively novel. We present the final ranked table of potency, drug-likeness, novelty, synthetic accessibility, and PubChem status, and then draw the selected molecular structures in a grid. We save the shortlist and figures to disk and close the workflow with an autonomous research summary that clearly separates computational hypotheses from experimentally validated drug candidates.

In conclusion, we completed the tutorial with a full in silico discovery loop that starts with public EGFR bioactivity data and ends with a prioritized shortlist of novel candidate molecules for further experimental investigation. We do not simply train a potency model; we curated the underlying chemistry, protected evaluation with a scaffold split, inspected model drivers, generated new analogs, and ranked them using a multi-parameter objective that reflects realistic medicinal chemistry trade-offs. The workflow also produces useful artifacts, including chemical-space plots, substructure-importance visualizations, candidate-structure grids, and a CSV shortlist, enabling us to review both the computational evidence and the proposed molecular designs. By keeping the pipeline CPU-friendly and API-key-free, we made advanced drug-discovery automation accessible within a standard Colab environment while still preserving scientific caution: the generated EGFR inhibitor hypotheses are not validated drugs and require docking, synthesis planning, selectivity profiling, and wet-lab assays before any real therapeutic claim can be made. Also, we demonstrated how an autonomous AI co-scientist can combine target intelligence, QSAR modeling, interpretability, fragment-based generation, and developability scoring into a coherent research workflow for resistance-aware kinase inhibitor discovery.

Check out the Full Codes here. Also, feel free to follow us on Twitter and don’t forget to join our 150k+ML SubReddit and Subscribe to our Newsletter. Wait! are you on telegram? now you can join us on telegram as well.

Need to partner with us for promoting your GitHub Repo OR Hugging Face Page OR Product Release OR Webinar etc.? Connect with us

The post Building a Scaffold-Split Random Forest QSAR Co-Scientist for EGFR Inhibitor Discovery Using ChEMBL, RDKit, SHAP, and BRICS appeared first on MarkTechPost.