Use this Skill for cheminformatics with RDKit: SMILES/InChI parsing, Morgan fingerprints, Tanimoto similarity, Murcko scaffold decomposition, substructure search, and chemical space visualization.
Use this Skill for cheminformatics with RDKit: SMILES/InChI parsing, Morgan fingerprints, Tanimoto similarity, Murcko scaffold decomposition, substructure search, and chemical space visualization.
TL;DR — Full cheminformatics pipeline with RDKit: parse SMILES/InChI/SDF,
compute molecular properties (Lipinski), generate Morgan fingerprints, calculate
Tanimoto similarity, decompose Murcko scaffolds, run SMARTS substructure searches,
and visualize chemical space with PCA.
When to Use
Use this Skill whenever you need to:
Parse molecular structures from SMILES strings, InChI identifiers, or SDF files
RDKit represents molecules as Mol objects containing atoms, bonds, ring systems,
and stereochemistry. The canonical entry points are:
SMILES (Simplified Molecular Input Line Entry System): human-readable ASCII string,
e.g. CCO for ethanol.
InChI (IUPAC International Chemical Identifier): standardized, canonical; preferred
for database exchange.
SDF/MOL (Structure Data File): 2D/3D coordinates + properties; standard for compound
libraries.
Morgan Fingerprints (ECFP)
Morgan/circular fingerprints encode the chemical environment around each atom up to a
given radius. radius=2 with nBits=2048 corresponds to ECFP4, the industry standard
for virtual screening and ML models.
Tanimoto Similarity
Tanimoto(A, B) = |A ∩ B| / |A ∪ B|
Values range from 0 (no common bits) to 1 (identical fingerprints). A threshold of
0.4 is commonly used for scaffold hopping; 0.85+ indicates very close analogs.
Murcko Scaffolds
The Murcko scaffold (Bemis-Murcko decomposition) strips side chains to retain the ring
systems and their connecting linkers. It is the standard framework for scaffold frequency
analysis in medicinal chemistry.
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors
from rdkit.Chem.Scaffolds import MurckoScaffold
from rdkit.Chem import Draw, DataStructs
import pandas as pd
import numpy as np
defmol_from_smiles(smiles: str):
"""Parse a SMILES string; returns None if invalid."""
mol = Chem.MolFromSmiles(smiles)
if mol isNone:
print(f"Invalid SMILES: {smiles}")
return mol
defmol_from_inchi(inchi: str):
"""Parse an InChI string to a RDKit Mol object."""
mol = Chem.MolFromInchi(inchi)
if mol isNone:
print(f"Invalid InChI: {inchi}")
return mol
defmols_from_sdf(sdf_path: str) -> list:
"""
Read all molecules from an SDF file.
Args:
sdf_path: Path to .sdf file.
Returns:
List of (mol, name) tuples; entries where mol is None are skipped.
"""
supplier = Chem.SDMolSupplier(sdf_path, removeHs=True, sanitize=True)
molecules = []
for mol in supplier:
if mol isnotNone:
name = mol.GetProp("_Name") if mol.HasProp("_Name") else"unnamed"
molecules.append((mol, name))
print(f"Loaded {len(molecules)} valid molecules from {sdf_path}")
return molecules
# Quick demo
smiles_list = [
"CC(=O)Oc1ccccc1C(=O)O", # Aspirin"CC12CCC3C(C1CCC2O)CCC4=CC(=O)CCC34C", # Testosterone"c1ccc2c(c1)cc1ccc3cccc4ccc2c1c34", # Pyrene"CN1C=NC2=C1C(=O)N(C(=O)N2C)C", # Caffeine"CC(C)Cc1ccc(cc1)C(C)C(=O)O", # Ibuprofen
]
mols = [mol_from_smiles(s) for s in smiles_list]
mols = [m for m in mols if m isnotNone]
print(f"Parsed {len(mols)} molecules")
Step 2 — Calculate Molecular Properties and Lipinski Filtering
defcalculate_properties(mol) -> dict:
"""
Calculate key physicochemical descriptors for drug-likeness assessment.
Returns:
Dictionary with MW, LogP, TPSA, HBD, HBA, RotBonds, RingCount,
and Lipinski Rule-of-Five compliance flag.
"""
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
tpsa = rdMolDescriptors.CalcTPSA(mol)
hbd = rdMolDescriptors.CalcNumHBD(mol) # H-bond donors
hba = rdMolDescriptors.CalcNumHBA(mol) # H-bond acceptors
rot = rdMolDescriptors.CalcNumRotatableBonds(mol)
rings = rdMolDescriptors.CalcNumRings(mol)
arom = rdMolDescriptors.CalcNumAromaticRings(mol)
# Lipinski Rule-of-Five (oral bioavailability proxy)
ro5 = (mw <= 500) and (logp <= 5) and (hbd <= 5) and (hba <= 10)
# Veber rules (add-on for permeability)
veber = (rot <= 10) and (tpsa <= 140)
return {
"MW": round(mw, 2),
"LogP": round(logp, 2),
"TPSA": round(tpsa, 2),
"HBD": hbd,
"HBA": hba,
"RotBonds": rot,
"Rings": rings,
"AromaticRings": arom,
"Lipinski_RO5": ro5,
"Veber": veber,
}
defbatch_properties(smiles_df: pd.DataFrame, smiles_col: str = "smiles") -> pd.DataFrame:
"""
Calculate properties for a DataFrame of SMILES strings.
Args:
smiles_df: DataFrame containing SMILES strings.
smiles_col: Column name with SMILES data.
Returns:
DataFrame with original columns plus all property columns.
"""
records = []
for _, row in smiles_df.iterrows():
smi = row[smiles_col]
mol = Chem.MolFromSmiles(str(smi))
if mol isNone:
records.append({k: Nonefor k in [
"MW", "LogP", "TPSA", "HBD", "HBA",
"RotBonds", "Rings", "AromaticRings", "Lipinski_RO5", "Veber"
]})
else:
records.append(calculate_properties(mol))
props_df = pd.DataFrame(records)
return pd.concat([smiles_df.reset_index(drop=True), props_df], axis=1)
# Demonstration
names = ["Aspirin", "Testosterone", "Pyrene", "Caffeine", "Ibuprofen"]
demo_df = pd.DataFrame({"name": names, "smiles": smiles_list})
result_df = batch_properties(demo_df)
print(result_df[["name", "MW", "LogP", "TPSA", "HBD", "HBA", "Lipinski_RO5"]].to_string(index=False))
Step 3 — Morgan Fingerprints and Tanimoto Similarity
defget_morgan_fp(mol, radius: int = 2, n_bits: int = 2048):
"""
Generate ECFP4-equivalent Morgan fingerprint as a bit vector.
Args:
mol: RDKit Mol object.
radius: Circular neighborhood radius (2 = ECFP4, 3 = ECFP6).
n_bits: Fingerprint bit-vector length.
Returns:
RDKit ExplicitBitVect object.
"""return AllChem.GetMorganFingerprintAsBitVect(mol, radius=radius, nBits=n_bits)
deftanimoto_matrix(mols: list) -> np.ndarray:
"""
Compute symmetric Tanimoto similarity matrix for a list of molecules.
Args:
mols: List of RDKit Mol objects.
Returns:
numpy array of shape (n, n) with pairwise Tanimoto values.
"""
fps = [get_morgan_fp(m) for m in mols]
n = len(fps)
matrix = np.ones((n, n))
for i inrange(n):
for j inrange(i + 1, n):
sim = DataStructs.TanimotoSimilarity(fps[i], fps[j])
matrix[i, j] = sim
matrix[j, i] = sim
return matrix
defbulk_similarity_search(
query_mol,
library_mols: list,
library_names: list,
threshold: float = 0.4,
radius: int = 2,
n_bits: int = 2048,
) -> pd.DataFrame:
"""
Find all molecules in a library with Tanimoto >= threshold to the query.
Args:
query_mol: Query RDKit Mol object.
library_mols: List of library RDKit Mol objects.
library_names: Names corresponding to library_mols.
threshold: Minimum Tanimoto similarity to include.
radius: Morgan radius.
n_bits: Fingerprint bits.
Returns:
DataFrame sorted by similarity descending, columns: name, smiles, tanimoto.
"""
query_fp = get_morgan_fp(query_mol, radius, n_bits)
lib_fps = [get_morgan_fp(m, radius, n_bits) for m in library_mols]
sims = DataStructs.BulkTanimotoSimilarity(query_fp, lib_fps)
hits = []
for name, mol, sim inzip(library_names, library_mols, sims):
if sim >= threshold:
hits.append({
"name": name,
"smiles": Chem.MolToSmiles(mol),
"tanimoto": round(sim, 4),
})
return pd.DataFrame(hits).sort_values("tanimoto", ascending=False).reset_index(drop=True)
# Similarity matrix example
sim_mat = tanimoto_matrix(mols)
print("Tanimoto similarity matrix:")
print(pd.DataFrame(sim_mat, index=names, columns=names).round(3).to_string())
Advanced Usage
Murcko Scaffold Decomposition
defget_murcko_scaffold(mol, generic: bool = False) -> str:
"""
Extract the Murcko scaffold SMILES from a molecule.
Args:
mol: RDKit Mol object.
generic: If True, return the generic scaffold (all atoms -> C, all bonds -> single).
Returns:
SMILES string of the scaffold, or '' if extraction fails.
"""try:
scaffold_mol = MurckoScaffold.GetScaffoldForMol(mol)
if generic:
scaffold_mol = MurckoScaffold.MakeScaffoldGeneric(scaffold_mol)
return Chem.MolToSmiles(scaffold_mol)
except Exception as e:
print(f"Scaffold error: {e}")
return""defscaffold_frequency_analysis(
smiles_list: list,
names: list = None,
top_n: int = 10,
) -> pd.DataFrame:
"""
Count scaffold frequencies across a compound library.
Args:
smiles_list: List of SMILES strings.
names: Optional compound names (same length as smiles_list).
top_n: Number of top scaffolds to return.
Returns:
DataFrame: scaffold_smiles, count, fraction, example_compound.
"""from collections import defaultdict
scaffold_to_compounds = defaultdict(list)
names = names or [f"cpd_{i}"for i inrange(len(smiles_list))]
for smi, name inzip(smiles_list, names):
mol = Chem.MolFromSmiles(smi)
if mol isNone:
continue
scaffold = get_murcko_scaffold(mol)
scaffold_to_compounds[scaffold].append(name)
records = []
total = len(smiles_list)
for scaffold_smi, compounds in scaffold_to_compounds.items():
records.append({
"scaffold_smiles": scaffold_smi,
"count": len(compounds),
"fraction": round(len(compounds) / total, 4),
"example_compound": compounds[0],
})
df = pd.DataFrame(records).sort_values("count", ascending=False).head(top_n)
return df.reset_index(drop=True)
# Example: scaffold analysis on a small library
scaffolds = scaffold_frequency_analysis(smiles_list, names)
print(scaffolds[["scaffold_smiles", "count", "fraction"]].to_string(index=False))
SMARTS Substructure Search
defsubstructure_search(
library_mols: list,
library_names: list,
smarts_pattern: str,
) -> pd.DataFrame:
"""
Filter a compound library by a SMARTS substructure pattern.
Args:
library_mols: List of RDKit Mol objects.
library_names: Names for each molecule.
smarts_pattern: SMARTS string to match against.
Returns:
DataFrame of matching compounds with match atom indices.
"""
query = Chem.MolFromSmarts(smarts_pattern)
if query isNone:
raise ValueError(f"Invalid SMARTS pattern: {smarts_pattern}")
hits = []
for name, mol inzip(library_names, library_mols):
if mol.HasSubstructMatch(query):
matches = mol.GetSubstructMatches(query)
hits.append({
"name": name,
"smiles": Chem.MolToSmiles(mol),
"n_matches": len(matches),
"match_atoms": str(matches[0]),
})
return pd.DataFrame(hits)
# Search for carboxylic acids
aromatic_ring_smarts = "c1ccccc1"# benzene ring
carboxylic_acid_smarts = "C(=O)[OH]"# carboxylic acid
sulfonamide_smarts = "S(=O)(=O)N"# sulfonamideprint("Benzene ring hits:")
print(substructure_search(mols, names, aromatic_ring_smarts))
print("\nCarboxylic acid hits:")
print(substructure_search(mols, names, carboxylic_acid_smarts))
Chemical Space Visualization with PCA
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
defchemical_space_pca(
mols: list,
labels: list = None,
color_by: list = None,
output_path: str = None,
) -> pd.DataFrame:
"""
Project molecules into 2D chemical space using PCA of Morgan fingerprints.
Args:
mols: List of RDKit Mol objects.
labels: Compound names for hover/annotation.
color_by: Numeric values for color-coding points (e.g., pIC50).
output_path: Path to save the PNG plot.
Returns:
DataFrame with columns: name, PC1, PC2, (color_value).
"""# Build fingerprint matrix
fps = [get_morgan_fp(m) for m in mols]
fp_array = np.array([list(fp.ToBitString()) for fp in fps], dtype=float)
# PCA
pca = PCA(n_components=2, random_state=42)
coords = pca.fit_transform(fp_array)
labels = labels or [f"cpd_{i}"for i inrange(len(mols))]
df_pca = pd.DataFrame({
"name": labels,
"PC1": coords[:, 0],
"PC2": coords[:, 1],
})
explained = pca.explained_variance_ratio_ * 100
fig, ax = plt.subplots(figsize=(8, 6))
scatter_kwargs = dict(s=60, alpha=0.8, edgecolors="k", linewidths=0.4)
if color_by isnotNone:
sc = ax.scatter(df_pca["PC1"], df_pca["PC2"], c=color_by,
cmap="RdYlGn", **scatter_kwargs)
plt.colorbar(sc, ax=ax, label="Property value")
df_pca["color_value"] = color_by
else:
ax.scatter(df_pca["PC1"], df_pca["PC2"], color="#4C72B0", **scatter_kwargs)
for _, row in df_pca.iterrows():
ax.annotate(row["name"], (row["PC1"], row["PC2"]),
fontsize=7, ha="center", va="bottom",
xytext=(0, 4), textcoords="offset points")
ax.set_xlabel(f"PC1 ({explained[0]:.1f}%)")
ax.set_ylabel(f"PC2 ({explained[1]:.1f}%)")
ax.set_title("Chemical Space — PCA of Morgan Fingerprints (ECFP4)")
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=150)
print(f"Saved PCA plot to {output_path}")
plt.show()
return df_pca
# Visualize chemical space with LogP as color
logp_values = [Descriptors.MolLogP(m) for m in mols]
pca_df = chemical_space_pca(mols, labels=names, color_by=logp_values,
output_path="chemical_space.png")
print(pca_df)
Examples
Example 1 — Batch Property Calculation from SMILES CSV
