""" Explain what a Bemis-Murcko scaffold actually is, using real molecules from the datasets in the post -- drawn by RDKit, not illustrated. Three figures: scaffold_chain.svg molecule -> Murcko framework -> generic framework scaffold_sink.svg several real, unlike molecules that all collapse to benzene scaffold_hole.svg real acyclic molecules, which have no scaffold at all RDKit draws in black; we rewrite that to `currentColor` so the figures follow the page's light/dark theme like every other figure in the post. """ import gzip import re import pandas as pd from rdkit import Chem, DataStructs, RDLogger from rdkit.Chem import rdFingerprintGenerator from rdkit.Chem.Draw import rdMolDraw2D from rdkit.Chem.Scaffolds import MurckoScaffold RDLogger.DisableLog("rdApp.*") MFP = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048) CELL_W, CELL_H = 210, 165 def draw(smiles, w=CELL_W, h=CELL_H): """Render one molecule to a theme-aware SVG fragment.""" mol = Chem.MolFromSmiles(smiles) if mol is None or mol.GetNumAtoms() == 0: return None d = rdMolDraw2D.MolDraw2DSVG(w, h) o = d.drawOptions() o.clearBackground = False o.useBWAtomPalette() o.bondLineWidth = 2 rdMolDraw2D.PrepareAndDrawMolecule(d, mol) d.FinishDrawing() svg = d.GetDrawingText() svg = svg[svg.index("", "", svg) # strip the outer ")] inner = inner.replace("#000000", "currentColor").replace("#000", "currentColor") inner = re.sub(r"]*style='[^']*fill:#FFFFFF[^']*'[^>]*/>", "", inner) return inner def cell(x, y, smiles, title, sub, w=CELL_W, h=CELL_H): g = [f''] g.append(f'') g.append(f'{title}') art = draw(smiles, w, h - 44) if art: g.append(f'{art}') else: g.append(f'' "(empty string)") g.append(f'{sub}') g.append("") return "\n".join(g) def arrow(x, y, label): return (f'' f'' f'' f'{label}' "") STYLE = """""" HEAD = ('") return "\n".join(o) # ---------------------------------------------------------------- fig 2 def pick_benzene_molecules(path, smi_col, k=4): """Real molecules from a real dataset whose Murcko scaffold is bare benzene.""" opener = gzip.open if path.endswith(".gz") else open with opener(path, "rt") as fh: df = pd.read_csv(fh) hits, seen = [], set() for s in df[smi_col].astype(str): m = Chem.MolFromSmiles(s) if m is None: continue cs = Chem.MolToSmiles(m) if cs in seen: continue seen.add(cs) try: if MurckoScaffold.MurckoScaffoldSmiles(mol=m) != "c1ccccc1": continue except Exception: continue hits.append((cs, m)) total = len(hits) # show only genuinely decorated members -- bare benzene itself is in the # group but makes a trivial-looking panel hits = [(cs, m) for cs, m in hits if cs != "c1ccccc1"] # choose maximally dissimilar members so the point is honest, not cherry-picked fps = [MFP.GetFingerprint(m) for _, m in hits] chosen = [0] while len(chosen) < k and len(chosen) < len(hits): best, bestscore = None, 2.0 for i in range(len(hits)): if i in chosen: continue worst = max(DataStructs.TanimotoSimilarity(fps[i], fps[j]) for j in chosen) if worst < bestscore: best, bestscore = i, worst chosen.append(best) return [hits[i][0] for i in chosen], total def fig_sink(): smis, total = pick_benzene_molecules("data/tox21.csv.gz", "smiles", k=4) W = 40 + 4 * (CELL_W + 12) + 260 H = CELL_H + 82 o = [HEAD.format(w=W, h=H, alt="Four chemically unlike Tox21 molecules that all " "reduce to the same benzene scaffold"), STYLE] o.append('These four Tox21 molecules are ' 'chemically unalike, and Murcko says they are the same group.') for i, s in enumerate(smis): o.append(cell(40 + i * (CELL_W + 12), 30, s, f"molecule {i+1}", s[:30])) ax = 40 + 4 * (CELL_W + 12) o.append(arrow(ax - 4, 30 + CELL_H / 2, "")) o.append(cell(ax + 44, 30, "c1ccccc1", "one scaffold", "c1ccccc1", w=140)) o.append(f'{total:,} Tox21 molecules land here') o.append("") return "\n".join(o) # ---------------------------------------------------------------- fig 3 def fig_hole(): opener = gzip.open("data/tox21.csv.gz", "rt") df = pd.read_csv(opener) acyc = [] for s in df["smiles"].astype(str): m = Chem.MolFromSmiles(s) if m is None or "." in s: continue try: if MurckoScaffold.MurckoScaffoldSmiles(mol=m) == "": acyc.append(Chem.MolToSmiles(m)) except Exception: continue if len(acyc) > 400: break picks = [s for s in acyc if 6 <= len(s) <= 26][:3] or acyc[:3] W = 40 + 3 * (CELL_W + 12) + 210 H = CELL_H + 76 o = [HEAD.format(w=W, h=H, alt="Real acyclic Tox21 molecules, which have no " "Murcko scaffold at all"), STYLE] o.append('No ring means no framework. ' 'RDKit returns the empty string, and every implementation then ' 'invents its own answer.') for i, s in enumerate(picks): o.append(cell(40 + i * (CELL_W + 12), 30, s, f"acyclic molecule", s[:30])) ax = 40 + 3 * (CELL_W + 12) o.append(arrow(ax - 4, 30 + CELL_H / 2, "")) o.append(cell(ax + 44, 30, "", "its scaffold", '""', w=160)) o.append("") return "\n".join(o) if __name__ == "__main__": for name, fn in [("scaffold_chain", fig_chain), ("scaffold_sink", fig_sink), ("scaffold_hole", fig_hole)]: open(f"{name}.svg", "w").write(fn()) print("wrote", name + ".svg")