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
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
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)")
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
banner("[4/9] CHEMICAL SPACE & SCAFFOLD ANALYSIS")
def murcko(m):
try:
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
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 < 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:
continue
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()))
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 <= 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"]:
try:
frags.update(BRICS.BRICSDecompose(m))
except Exception:
pass
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"])
gen = {}
try:
for i, prod in enumerate(BRICS.BRICSBuild(frag_mols, scrambleReagents=True, maxDepth=2)):
if i >= 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 <= 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
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
if x < lo:
return float(np.clip((x - hl) / (lo - hl + 1e-9), 0, 1))
if x > 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 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
banner("[9/9] NOVELTY CROSS-CHECK (PubChem) & FINAL SHORTLIST")
def pubchem_cid(smi):
m = Chem.MolFromSmiles(smi)
try:
ik = Chem.MolToInchiKey(m)
except Exception:
return "inchikey_error"
if not ik:
return "inchikey_error"
js = http_json(f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ik}/cids/JSON")
time.sleep(0.25)
try:
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))
try:
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")
try:
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.
Conclusion
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.
from MarkTechPost https://ift.tt/bDNtwSz
via IFTTT
Comments
Post a Comment