Building a Scaffold-Split Random Forest QSAR Co-Scientist for EGFR Inhibitor Discovery Using ChEMBL, RDKit, SHAP, and BRICS
In this tutorial, we construct an end-to-end autonomous AI co-scientist workflow for next-generation EGFR inhibitor discovery, specializing in the C797S osimertinib-resistance mutation in non-small cell lung most cancers. We begin by resolving the organic goal via ChEMBL and UniProt, then mine curated EGFR IC50 bioactivity data and convert them into a clear pIC50 modeling dataset. We use RDKit to standardize molecules, take away salts, combination replicate measurements, compute Morgan fingerprints, extract physicochemical descriptors, and analyze scaffold variety in order that our mannequin learns from chemically significant representations reasonably than uncooked strings. From there, we prepare a scaffold-split Random Forest QSAR mannequin, consider its potential to generalize to unseen chemotypes, interpret potency-driving options with SHAP or mannequin importances, and visualize influential molecular substructures. Finally, we transfer past prediction into generative design by recombining BRICS fragments from potent drug-like actives, scoring the ensuing digital analogs throughout efficiency, drug-likeness, synthesizability, novelty, and developability gates, and cross-checking the shortlisted candidates towards 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], test=False)
for mod, pkg in [("rdkit", "rdkit"), ("shap", "shap"), ("requests", "requests")]:
attempt:
importlib.import_module(mod)
besides 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.*")
attempt:
from rdkit.Chem.MolStandardize import rdMolStandardize
_HAS_STD = True
besides Exception:
_HAS_STD = False
attempt:
from rdkit.Chem import RDConfig
sys.path.append(os.path.be a part of(RDConfig.RDContribDir, "SA_Score"))
import sascorer
_HAS_SA = True
besides Exception:
_HAS_SA = False
TARGET_CHEMBL_ID = "CHEMBL203"
TARGET_QUERY = "Epidermal development issue 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/knowledge"
HDRS = {"Accept": "utility/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 okay in vary(tries):
attempt:
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
besides Exception:
cross
time.sleep(1.5 * (okay + 1))
return None
banner("[1/9] TARGET INTELLIGENCE (ChEMBL + UniProt)")
print("Question: What goal are we drugging, and why is it laborious?n")
def ic50_count(tid):
js = http_json(f"{BASE}/exercise", {"target_chembl_id": tid, "standard_type": "IC50",
"pchembl_value__isnull": "false", "restrict": 1, "format": "json"})
attempt:
return int(js["page_meta"]["total_count"])
besides 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}/goal/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("rating", 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:
greatest = max(scored, key=lambda x: x[1])[0]
target_id, target_name = greatest["target_chembl_id"], greatest.get("pref_name", TARGET_QUERY)
print(f" Auto-resolved '{TARGET_QUERY}' by knowledge quantity -> {target_id}")
else:
target_id = FALLBACK_CHEMBL_ID
print(f" Auto-resolve discovered no knowledge; falling again to {FALLBACK_CHEMBL_ID}")
det = http_json(f"{BASE}/goal/{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 goal : {target_name}")
print(f" ChEMBL ID : {target_id}")
print(f" UniProt : {uniprot_acc}")
if uniprot_acc:
uni = http_json(f"https://relaxation.uniprot.org/uniprotkb/{uniprot_acc}.json")
if uni:
attempt:
fn = subsequent(c["texts"][0]["value"] for c in uni.get("feedback", [])
if c.get("commentType") == "FUNCTION")
print("n Function (UniProt):")
print(" ", (fn[:340] + " ...") if len(fn) > 340 else fn)
besides Exception:
cross
print("""
Resistance context: 1st/2nd/Third-gen EGFR TKIs lose efficiency as soon as tumours purchase the
C797S mutation, which abolishes the covalent cysteine anchor exploited by osimertinib.
Goal of this run: be taught the chemistry of recognized EGFR inhibitors and suggest NOVEL,
drug-like analogs as beginning factors for a C797S-active 4th-generation collection.""")
We start by getting ready the complete scientific computing atmosphere and putting in any lacking chemistry, modeling, plotting, and API dependencies required by the workflow. We configure the EGFR goal settings, outline modeling constants, initialize reproducible random seeds, and create helper features for banners and strong JSON API calls. We then resolve the ChEMBL goal, retrieve UniProt context when accessible, and body the organic motivation round EGFR C797S resistance.
Mining ChEMBL Bioactivity Data
banner("[2/9] BIOACTIVITY MINING (ChEMBL actions -> pIC50)")
def pull_activities(tid, cap):
url, rows = f"{BASE}/exercise", []
params = {"target_chembl_id": tid, "standard_type": "IC50",
"pchembl_value__isnull": "false", "restrict": 1000, "format": "json"}
js = http_json(url, params)
pages = 0
whereas js and pages < 60:
rows.prolong(js.get("actions", []))
pages += 1
if len(rows) >= cap:
break
nxt = js.get("page_meta", {}).get("subsequent")
if not nxt:
break
nurl = nxt if nxt.startswith("http") else "https://www.ebi.ac.uk" + nxt
js = http_json(nurl)
return rows[:cap]
uncooked = pull_activities(target_id, MAX_ACTIVITIES)
print(f" Pulled {len(uncooked)} uncooked IC50 data with a curated pChEMBL worth.")
recs = []
for a in uncooked:
smi, pv = a.get("canonical_smiles"), a.get("pchembl_value")
if not smi or pv in (None, ""):
proceed
if a.get("standard_relation") != "=":
proceed
if a.get("standard_units") not in ("nM", None):
proceed
attempt:
pv = float(pv)
besides Exception:
proceed
recs.append({"chembl_id": a.get("molecule_chembl_id"), "smiles": smi, "pIC50": pv})
raw_df = pd.DataBody(recs, columns=["chembl_id", "smiles", "pIC50"])
print(f" After high quality filters: {len(raw_df)} measurements.")
if len(raw_df) == 0:
print("n STOP: no usable IC50 knowledge was retrieved for this goal.n"
" Fix: set TARGET_CHEMBL_ID to a goal that has inhibitor datan"
" (e.g. CHEMBL203=EGFR, CHEMBL5251=BTK, CHEMBL2971=JAK2),n"
" or set TARGET_CHEMBL_ID="" to auto-resolve TARGET_QUERY by knowledge quantity.")
increase SystemExit("No bioactivity knowledge for the chosen goal.")
banner("[3/9] MOLECULAR CURATION (standardize, de-salt, combination)")
_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
attempt:
if _HAS_STD:
m = _lfc.select(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)
besides 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:
proceed
canon.append({"smiles": cs, "pIC50": r["pIC50"], "chembl_id": r["chembl_id"]})
keep_mol[cs] = m
cdf = pd.DataBody(canon, columns=["smiles", "pIC50", "chembl_id"])
knowledge = (cdf.groupby("smiles")
.agg(pIC50=("pIC50", "median"), n=("pIC50", "measurement"),
chembl_id=("chembl_id", "first")).reset_index())
if len(knowledge) > MAX_UNIQUE:
knowledge = knowledge.pattern(MAX_UNIQUE, random_state=RANDOM_STATE).reset_index(drop=True)
knowledge["mol"] = knowledge["smiles"].map(keep_mol)
n_active = int((knowledge["pIC50"] >= ACTIVE_PIC50).sum())
print(f" Unique curated molecules : {len(knowledge)}")
print(f" Potent actives (IC50<=100nM): {n_active} ({100*n_active/len(knowledge):.1f}%)")
print(f" pIC50 vary: {knowledge.pIC50.min():.2f} - {knowledge.pIC50.max():.2f} "
f"(median {knowledge.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):
attempt:
Xd[i, j] = fn(m)
besides Exception:
Xd[i, j] = 0.0
return np.nan_to_num(np.hstack([Xb, Xd]).astype(np.float32))
X = featurize(listing(knowledge["mol"]))
y = knowledge["pIC50"].values
print(f" Feature matrix: {X.form[0]} molecules x {X.form[1]} options "
f"({NBITS} ECFP bits + {len(DESC)} descriptors)")
We mine curated IC50 bioactivity measurements from ChEMBL and convert the uncooked exercise data into a usable pIC50 dataset. We filter out incomplete, non-exact, or inconsistent measurements in order that the downstream QSAR mannequin trains on cleaner efficiency values. We then standardize the molecules with RDKit, take away salts or smaller fragments, combination duplicate molecules by median pIC50, and convert every molecule into Morgan fingerprint bits plus interpretable physicochemical descriptors.
Scaffold Analysis and QSAR
banner("[4/9] CHEMICAL SPACE & SCAFFOLD ANALYSIS")
def murcko(m):
attempt:
s = MurckoScaffold.GetScaffoldForMol(m)
cs = Chem.MolToSmiles(s)
return cs if cs else "(acyclic)"
besides Exception:
return "(error)"
knowledge["scaffold"] = knowledge["mol"].map(murcko)
top_scaf = knowledge["scaffold"].value_counts().head(10)
print(" Top recurring Murcko scaffolds (chemotype households):")
for i, (sc, c) in enumerate(top_scaf.gadgets(), 1):
print(f" {i:2nd}. 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 house (ECFP-PCA), colored by efficiency",
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, coloration="#3b7dd8", edgecolor="white")
ax[1].axvline(ACTIVE_PIC50, coloration="crimson", ls="--", label=f"energetic minimize (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], coloration="#6c5ce7")
ax[2].set(title="Top 10 scaffolds", xlabel="rely")
plt.tight_layout(); plt.savefig("fig1_chemical_space.png", dpi=120); plt.present()
banner("[5/9] INTERPRETABLE QSAR MODEL (scaffold-split, leakage-free)")
def scaffold_split(scaffolds, test_frac=0.2, seed=RANDOM_STATE):
teams = {}
for i, s in enumerate(scaffolds):
teams.setdefault(s, []).append(i)
order = listing(teams.values())
random.Random(seed).shuffle(order)
n_test = int(len(scaffolds) * test_frac)
take a look at, prepare = [], []
for g so as:
(take a look at if len(take a look at) < n_test else prepare).prolong(g)
return np.array(sorted(prepare)), np.array(sorted(take a look at))
tr, te = scaffold_split(knowledge["scaffold"].values, 0.2)
mannequin = RandomForestRegressor(n_estimators=400, max_features="sqrt",
n_jobs=-1, random_state=RANDOM_STATE)
mannequin.match(X[tr], y[tr])
pred = mannequin.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.distinctive(ycls)) == 2 else float("nan")
print(f" Held-out (new-scaffold) efficiency on {len(te)} molecules:")
print(f" R^2 = {r2:.3f}")
print(f" RMSE (pIC50) = {rmse:.3f} (~{rmse:.2f} log items)")
print(f" Spearman rho = {rho:.3f}")
print(f" ROC-AUC energetic = {auc:.3f} (rating potent vs weak)")
model_full = RandomForestRegressor(n_estimators=400, max_features="sqrt",
n_jobs=-1, random_state=RANDOM_STATE).match(X, y)
We analyze the curated chemical house by extracting Murcko scaffolds and figuring out the most typical chemotype households within the dataset. We challenge Morgan fingerprints with PCA to visualise how EGFR inhibitors distribute throughout chemical house and how efficiency varies throughout that panorama. We then prepare a Random Forest QSAR mannequin utilizing a scaffold break up, which helps us consider whether or not the mannequin generalizes to unseen molecular scaffolds reasonably than memorizing shut analogs.
Interpretability and BRICS Generation
banner("[6/9] MODEL INTERPRETABILITY (which substructures drive efficiency?)")
top_feat_idx, shap_ok = None, False
attempt:
import shap
samp = np.random.RandomState(RANDOM_STATE).alternative(len(te), min(300, len(te)), substitute=False)
expl = shap.TreeExplainer(mannequin)
sv = expl.shap_values(X[te][samp])
if isinstance(sv, listing):
sv = sv[0]
imp = np.abs(sv).imply(0)
shap_ok = True
print(" Using SHAP TreeExplainer (imply |SHAP| over held-out molecules).")
besides Exception as e:
imp = model_full.feature_importances_
print(f" SHAP unavailable ({sort(e).__name__}); utilizing 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 so as if i >= NBITS][:8]
top_bits = [i for i in order if i < NBITS][:6]
print("n Most influential physicochemical descriptors:")
for identify, v in top_desc:
print(f" {identify:12s} significance={v:.4f}")
train_smis = listing(knowledge["smiles"].iloc[tr])
def bit_exemplar(bit):
for smi in train_smis:
m = Chem.MolFromSmiles(smi)
if m is None:
proceed
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 = listing(env)
for bidx in env:
b = m.GetBondWithIdx(bidx)
atoms.replace((b.GetBeginAtomIdx(), b.GetEndAtomIdx()))
attempt:
return Draw.MolToImage(m, measurement=(300, 240),
highlightAtoms=listing(atoms), highlightBonds=bonds)
besides Exception:
return None
return None
attempt:
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 mannequin associates with efficiency", y=1.02)
plt.tight_layout(); plt.savefig("fig2_potency_substructures.png", dpi=120, bbox_inches="tight")
plt.present()
besides Exception as e:
print(f" (substructure drawing skipped: {sort(e).__name__})")
banner("[7/9] GENERATIVE DESIGN (BRICS fragment recombination -> novel analogs)")
seed = knowledge[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 father or mother molecules.")
frags = set()
for m in seed["mol"]:
attempt:
frags.replace(BRICS.BRICSDecompose(m))
besides Exception:
cross
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.")
recognized = set(knowledge["smiles"])
gen = {}
attempt:
for i, prod in enumerate(BRICS.BRICSBuild(frag_mols, scrambleReagents=True, maxDepth=2)):
if i >= BRICS_MAX_TRIES:
break
attempt:
prod.UpdatePropertyCache(strict=False)
Chem.SanitizeMol(prod)
cs = Chem.MolToSmiles(prod)
besides Exception:
proceed
if cs in recognized or cs in gen:
proceed
mw = Descriptors.MolWt(prod)
if 250 <= mw <= 600 and 8 <= prod.GetNumHeavyAtoms() <= 45:
gen[cs] = prod
besides Exception as e:
print(f" (BRICS construct ended early: {sort(e).__name__})")
print(f" Generated {len(gen)} distinctive, novel, size-reasonable digital molecules.")
We interpret the skilled QSAR mannequin by estimating which descriptors and fingerprint bits contribute most strongly to predicted efficiency. We use SHAP when accessible; in any other case, we fall again to Random Forest function importances to maintain the workflow strong. We additionally visualize consultant molecular substructures related to influential ECFP bits, then start generative design by decomposing potent drug-like father or mother molecules into BRICS fragments and recombining them to generate novel digital analogs.
Multi-Parameter Candidate Scoring
banner("[8/9] MULTI-PARAMETER PRIORITISATION")
gsmiles = listing(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, hello, hard_lo=None, hard_hi=None):
hl = hard_lo if hard_lo just isn't None else lo
hh = hard_hi if hard_hi just isn't None else hello
if x < lo:
return float(np.clip((x - hl) / (lo - hl + 1e-9), 0, 1))
if x > hello:
return float(np.clip((hh - x) / (hh - hello + 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)
attempt:
qed = QED.qed(m)
besides 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)
rating = (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, rating=rating))
_CANDCOLS = ["smiles", "pred_pIC50", "MolWt", "MolLogP", "TPSA", "HBD", "HBA",
"QED", "SA", "novelty", "lipinski", "veber_ok", "score"]
cand = pd.DataBody(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("rating", ascending=False).reset_index(drop=True)
print(f" {len(gate)} of {len(cand)} generated molecules handed 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 recognized EGFR inhibitors)n")
shortlist = gate.head(N_SHORTLIST).copy()
We rating the generated molecules with the complete QSAR mannequin and consider them throughout efficiency, molecular weight, lipophilicity, hydrogen bonding, polar floor space, rotatable bonds, QED, artificial accessibility, and novelty. We calculate a multi-parameter desirability rating that balances predicted efficiency with drug-likeness, synthesizability, and structural distance from recognized EGFR inhibitors. We then apply laborious developability gates and maintain solely the strongest candidates for the ultimate shortlist.
PubChem Novelty Cross-Check
banner("[9/9] NOVELTY CROSS-CHECK (PubChem) & FINAL SHORTLIST")
def pubchem_cid(smi):
m = Chem.MolFromSmiles(smi)
attempt:
ik = Chem.MolToInchiKey(m)
besides Exception:
return "inchikey_error"
if not ik:
return "inchikey_error"
js = http_json(f"https://pubchem.ncbi.nlm.nih.gov/relaxation/pug/compound/inchikey/{ik}/cids/JSON")
time.sleep(0.25)
attempt:
return f"CID {js['IdentifierList']['CID'][0]}"
besides 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"]
fairly = shortlist[show_cols].copy()
fairly.insert(0, "rank", vary(1, len(fairly) + 1))
pd.set_option("show.width", 200, "show.max_colwidth", 40)
print("nTOP CANDIDATE 4th-GENERATION EGFR-INHIBITOR STARTING POINTS:n")
print(fairly.spherical(3).to_string(index=False))
attempt:
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")
attempt:
from IPython.show import show
show(grid)
besides Exception:
plt.determine(figsize=(14, 9)); plt.imshow(grid); plt.axis("off"); plt.present()
besides Exception as e:
print(f" (candidate drawing skipped: {sort(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(knowledge)} curated, de-duplicated EGFR inhibitors from ChEMBL
Learned mannequin : RandomForest QSAR, scaffold-split R^2={r2:.2f}, ROC-AUC={auc:.2f}
-> generalises to unseen chemotypes, not memorising analogs
Key drivers : {", ".be a part of(n for n, _ in top_desc[:4])} + particular ECFP substructures
Invented : {len(gen)} novel digital analogs by way of BRICS fragment recombination
Prioritised : {len(gate)} handed developability gates; prime {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 house + efficiency + scaffolds)
- fig2_potency_substructures.png (SHAP-implicated substructures)
- fig3_top_candidates.png (constructions of the shortlist)
NEXT STEPS a wet-lab workforce would take: dock the shortlist into the EGFR(L858R/T790M/C797S)
triple-mutant construction, prioritise allosteric binders, test artificial routes, and
assay the highest ~5 for C797S efficiency and selectivity vs wild-type EGFR.
Reminder: that is an academic in-silico speculation generator, not a validated drug
pipeline. Predictions require experimental affirmation.""")
print("nDONE.")
We cross-check the shortlisted molecules towards PubChem utilizing an InChIKey lookup to find out whether or not every candidate is understood or putatively novel. We current the ultimate ranked desk of efficiency, drug-likeness, novelty, artificial accessibility, and PubChem standing, and then draw the chosen molecular constructions in a grid. We save the shortlist and figures to disk and shut the workflow with an autonomous analysis abstract that clearly separates computational hypotheses from experimentally validated drug candidates.
Conclusion
In conclusion, we accomplished the tutorial with a full in silico discovery loop that begins with public EGFR bioactivity knowledge and ends with a prioritized shortlist of novel candidate molecules for additional experimental investigation. We don’t merely prepare a efficiency mannequin; we curated the underlying chemistry, protected analysis with a scaffold break up, inspected mannequin drivers, generated new analogs, and ranked them utilizing a multi-parameter goal that displays life like medicinal chemistry trade-offs. The workflow additionally produces helpful artifacts, together with chemical-space plots, substructure-importance visualizations, candidate-structure grids, and a CSV shortlist, enabling us to overview each the computational proof and the proposed molecular designs. By conserving the pipeline CPU-friendly and API-key-free, we made superior drug-discovery automation accessible inside a customary Colab atmosphere whereas nonetheless preserving scientific warning: the generated EGFR inhibitor hypotheses should not validated medicine and require docking, synthesis planning, selectivity profiling, and wet-lab assays earlier than any actual therapeutic declare could be made. Also, we demonstrated how an autonomous AI co-scientist can mix goal intelligence, QSAR modeling, interpretability, fragment-based era, and developability scoring into a coherent analysis workflow for resistance-aware kinase inhibitor discovery.
Check out the Full Codes here. Also, be at liberty to observe us on Twitter and don’t overlook to hitch 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 accomplice with us for selling your GitHub Repo OR Hugging Face Page OR Product Release OR Webinar and so on.? Connect with us
The publish Building a Scaffold-Split Random Forest QSAR Co-Scientist for EGFR Inhibitor Discovery Using ChEMBL, RDKit, SHAP, and BRICS appeared first on MarkTechPost.
