RDKit - Cheminformatics and Drug Discovery
RDKit is the industry-standard open-source toolkit for cheminformatics. It provides comprehensive tools for molecular manipulation, descriptor calculation, fingerprinting, substructure searching, and 3D molecular modeling. RDKit is used extensively in pharmaceutical companies for drug discovery and virtual screening.
When to Use
- Reading and writing chemical file formats (SMILES, SDF, MOL2, PDB).
- Calculating molecular descriptors and drug-like properties (Lipinski's Rule of Five).
- Generating molecular fingerprints for similarity searching.
- Substructure searching and chemical pattern matching (SMARTS).
- 3D conformer generation and molecular alignment.
- Virtual screening of compound libraries.
- Pharmacophore modeling and shape similarity.
- QSAR (Quantitative Structure-Activity Relationship) modeling.
- Reaction enumeration and retrosynthesis.
- Visualizing chemical structures in 2D and 3D.
- Building machine learning models for molecular property prediction.
Reference Documentation
Official docs: https://www.rdkit.org/docs/ RDKit Book: https://www.rdkit.org/docs/RDKit_Book.html GitHub: https://github.com/rdkit/rdkit Search patterns: rdkit.Chem, rdkit.Chem.Descriptors, rdkit.Chem.AllChem, rdkit.DataStructs
Core Principles
Molecular Representation
RDKit represents molecules as graphs where atoms are nodes and bonds are edges. The core object is Mol, which can be created from SMILES, SDF files, or built programmatically.
SMILES (Simplified Molecular Input Line Entry System)
A text-based notation for chemical structures. Example: CCO is ethanol, c1ccccc1 is benzene. RDKit can parse and generate SMILES strings.
Fingerprints for Similarity
Molecular fingerprints are binary vectors encoding structural features. They enable fast similarity searching and clustering of large compound libraries.
Lazy Evaluation
Many RDKit operations are lazy - properties are computed only when needed. This makes operations on large libraries very efficient.
Quick Reference
Installation
# Via conda (recommended)
conda install -c conda-forge rdkit
# Via pip
pip install rdkit
# For visualization
pip install rdkit pillowStandard Imports
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Draw, Lipinski
from rdkit.Chem import rdFingerprintGenerator
from rdkit import DataStructs
import numpy as np
import pandas as pdBasic Pattern - SMILES to Molecule
from rdkit import Chem
# 1. Create molecule from SMILES
smiles = "CC(=O)OC1=CC=CC=C1C(=O)O" # Aspirin
mol = Chem.MolFromSmiles(smiles)
# 2. Check if molecule is valid
if mol is None:
print("Invalid SMILES")
else:
print(f"Molecular formula: {Chem.rdMolDescriptors.CalcMolFormula(mol)}")
print(f"Molecular weight: {Descriptors.MolWt(mol):.2f}")Basic Pattern - Calculate Properties
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")
# Calculate drug-like properties
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)
print(f"MW: {mw:.2f}, LogP: {logp:.2f}, HBD: {hbd}, HBA: {hba}")
# Check Lipinski's Rule of Five
lipinski_pass = (mw <= 500 and logp <= 5 and hbd <= 5 and hba <= 10)
print(f"Lipinski compliant: {lipinski_pass}")Critical Rules
✅ DO
- Always Validate Molecules - Check
mol is not Noneafter parsing SMILES/files to catch invalid structures. - Use Canonical SMILES - Use
Chem.MolToSmiles(mol)to get canonical (standardized) SMILES for comparison. - Sanitize Molecules - RDKit auto-sanitizes by default (valence checking, aromaticity). Keep it enabled unless you have a specific reason.
- Generate 3D Coordinates - Use
AllChem.EmbedMolecule()before 3D operations like alignment or docking. - Use Fingerprints for Large Libraries - For similarity searching in millions of compounds, fingerprints are 1000x faster than direct comparison.
- Specify Random Seeds - For reproducible conformer generation, always set
randomSeed. - Handle Stereochemistry - Use
Chem.AssignStereochemistry()to properly assign R/S and E/Z labels. - Batch Processing - Use generators or chunking for processing millions of molecules to avoid memory issues.
❌ DON'T
- Don't Ignore Invalid Molecules - Always handle the case when
MolFromSmiles()returnsNone. - Don't Compare SMILES Strings Directly - Two different SMILES can represent the same molecule. Use canonical SMILES or InChI.
- Don't Skip Kekulization - For aromatic systems, ensure proper Kekulé structure assignment.
- Don't Use Descriptors for Similarity - Use fingerprints (much faster and more appropriate).
- Don't Forget Hydrogens - Add explicit hydrogens with
Chem.AddHs()when needed for 3D operations. - Don't Overuse 3D Minimization - Energy minimization is slow; only use when necessary (docking, visualization).
Anti-Patterns (NEVER)
from rdkit import Chem
from rdkit.Chem import AllChem
# ❌ BAD: Not checking if molecule is valid
smiles = "INVALID_SMILES"
mol = Chem.MolFromSmiles(smiles)
mw = Descriptors.MolWt(mol) # Crashes!
# ✅ GOOD: Always validate
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
mw = Descriptors.MolWt(mol)
else:
print("Invalid SMILES")
# ❌ BAD: Comparing SMILES strings directly
smiles1 = "CC(C)C" # isobutane
smiles2 = "C(C)CC" # same molecule, different SMILES
if smiles1 == smiles2: # False, but same molecule!
print("Same")
# ✅ GOOD: Use canonical SMILES
mol1 = Chem.MolFromSmiles(smiles1)
mol2 = Chem.MolFromSmiles(smiles2)
can1 = Chem.MolToSmiles(mol1)
can2 = Chem.MolToSmiles(mol2)
if can1 == can2: # True
print("Same molecule")
# ❌ BAD: 3D operations without 3D coordinates
mol = Chem.MolFromSmiles("CCO")
AllChem.AlignMol(mol, ref_mol) # Fails! No 3D coords
# ✅ GOOD: Generate 3D coordinates first
mol = Chem.MolFromSmiles("CCO")
AllChem.EmbedMolecule(mol)
AllChem.AlignMol(mol, ref_mol)Molecular I/O and Conversion
SMILES Parsing
from rdkit import Chem
# Parse SMILES
mol = Chem.MolFromSmiles("CCO")
# Parse SMILES with sanitization control
mol = Chem.MolFromSmiles("CCO", sanitize=True) # Default
# Generate canonical SMILES
canonical = Chem.MolToSmiles(mol)
# Generate isomeric SMILES (includes stereochemistry)
iso_smiles = Chem.MolToSmiles(mol, isomericSmiles=True)
# Generate SMILES without stereochemistry
non_iso = Chem.MolToSmiles(mol, isomericSmiles=False)
# Handle invalid SMILES
smiles_list = ["CCO", "INVALID", "c1ccccc1"]
mols = []
for smi in smiles_list:
mol = Chem.MolFromSmiles(smi)
if mol is not None:
mols.append(mol)
else:
print(f"Failed to parse: {smi}")Reading SDF Files
from rdkit import Chem
# Read single molecule from file
mol = Chem.MolFromMolFile("molecule.mol")
# Read multiple molecules from SDF
suppl = Chem.SDMolSupplier("compounds.sdf")
# Iterate through molecules
for mol in suppl:
if mol is None:
continue
smiles = Chem.MolToSmiles(mol)
print(f"SMILES: {smiles}")
# Access SDF properties
if mol.HasProp("_Name"):
name = mol.GetProp("_Name")
print(f"Name: {name}")
# Read with removeHs=False to keep explicit hydrogens
suppl = Chem.SDMolSupplier("compounds.sdf", removeHs=False)Writing SDF Files
from rdkit import Chem
# Write single molecule
mol = Chem.MolFromSmiles("CCO")
writer = Chem.SDWriter("output.sdf")
writer.write(mol)
writer.close()
# Write multiple molecules
mols = [Chem.MolFromSmiles(s) for s in ["CCO", "c1ccccc1", "CC(=O)O"]]
writer = Chem.SDWriter("output.sdf")
for mol in mols:
if mol is not None:
writer.write(mol)
writer.close()
# Add properties to molecules
mol = Chem.MolFromSmiles("CCO")
mol.SetProp("_Name", "Ethanol")
mol.SetProp("Activity", "10.5")
writer = Chem.SDWriter("output.sdf")
writer.write(mol)
writer.close()InChI and InChIKey
from rdkit import Chem
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O") # Aspirin
# Generate InChI (unique chemical identifier)
inchi = Chem.MolToInchi(mol)
print(f"InChI: {inchi}")
# Generate InChIKey (hashed InChI, good for database lookups)
inchikey = Chem.MolToInchiKey(mol)
print(f"InChIKey: {inchikey}")
# Parse InChI
mol_from_inchi = Chem.MolFromInchi(inchi)Molecular Descriptors
Common Descriptors
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski, Crippen
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O") # Aspirin
# Basic properties
mw = Descriptors.MolWt(mol)
num_atoms = mol.GetNumAtoms()
num_heavy_atoms = Lipinski.HeavyAtomCount(mol)
# Lipinski's Rule of Five parameters
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)
rotatable_bonds = Lipinski.NumRotatableBonds(mol)
# Topological descriptors
tpsa = Descriptors.TPSA(mol) # Topological polar surface area
rings = Lipinski.RingCount(mol)
aromatic_rings = Lipinski.NumAromaticRings(mol)
# Complexity
bertz_ct = Descriptors.BertzCT(mol) # Molecular complexity
print(f"""
Molecular Weight: {mw:.2f}
LogP: {logp:.2f}
HBD: {hbd}
HBA: {hba}
TPSA: {tpsa:.2f}
Rotatable Bonds: {rotatable_bonds}
Aromatic Rings: {aromatic_rings}
""")Calculate All Descriptors
from rdkit import Chem
from rdkit.Chem import Descriptors
mol = Chem.MolFromSmiles("CCO")
# Get all available descriptors
descriptor_names = [desc[0] for desc in Descriptors.descList]
# Calculate all descriptors
descriptors = {}
for name in descriptor_names:
calc = getattr(Descriptors, name)
descriptors[name] = calc(mol)
print(f"Total descriptors: {len(descriptors)}")
print(f"First 5: {list(descriptors.items())[:5]}")Drug-Likeness Filters
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski
def check_lipinski(mol):
"""Check Lipinski's Rule of Five."""
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)
rules = {
'MW <= 500': mw <= 500,
'LogP <= 5': logp <= 5,
'HBD <= 5': hbd <= 5,
'HBA <= 10': hba <= 10
}
passed = all(rules.values())
return passed, rules
def check_veber(mol):
"""Check Veber's rules for oral bioavailability."""
rotatable = Lipinski.NumRotatableBonds(mol)
tpsa = Descriptors.TPSA(mol)
rules = {
'Rotatable bonds <= 10': rotatable <= 10,
'TPSA <= 140': tpsa <= 140
}
passed = all(rules.values())
return passed, rules
# Usage
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O") # Aspirin
lipinski_pass, lipinski_rules = check_lipinski(mol)
veber_pass, veber_rules = check_veber(mol)
print(f"Lipinski: {lipinski_pass}")
print(f"Veber: {veber_pass}")Molecular Fingerprints
Morgan Fingerprints (Circular)
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs
# Create molecule
mol = Chem.MolFromSmiles("CCO")
# Generate Morgan fingerprint (ECFP4)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
# Convert to numpy array
import numpy as np
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
# Generate count-based fingerprint (for feature importance)
fp_counts = AllChem.GetMorganFingerprint(mol, radius=2)
# Get feature info
info = {}
fp_info = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048, bitInfo=info)
print(f"Number of on-bits: {len(info)}")Fingerprint Similarity
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs
# Create molecules
mol1 = Chem.MolFromSmiles("CCO")
mol2 = Chem.MolFromSmiles("CCCO")
mol3 = Chem.MolFromSmiles("c1ccccc1")
# Generate fingerprints
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, radius=2)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, radius=2)
fp3 = AllChem.GetMorganFingerprintAsBitVect(mol3, radius=2)
# Calculate Tanimoto similarity
sim_12 = DataStructs.TanimotoSimilarity(fp1, fp2)
sim_13 = DataStructs.TanimotoSimilarity(fp1, fp3)
print(f"Ethanol vs Propanol: {sim_12:.3f}") # High similarity
print(f"Ethanol vs Benzene: {sim_13:.3f}") # Low similarity
# Calculate Dice similarity
dice_12 = DataStructs.DiceSimilarity(fp1, fp2)
# Bulk similarity (compare one to many)
fps = [fp1, fp2, fp3]
similarities = DataStructs.BulkTanimotoSimilarity(fp1, fps)
print(f"Bulk similarities: {similarities}")Other Fingerprint Types
from rdkit import Chem
from rdkit.Chem import AllChem, RDKFingerprint
mol = Chem.MolFromSmiles("CCO")
# RDKit fingerprint (topological)
fp_rdkit = Chem.RDKFingerprint(mol)
# Atom pair fingerprint
fp_atompair = AllChem.GetHashedAtomPairFingerprintAsBitVect(mol)
# Topological torsion fingerprint
fp_torsion = AllChem.GetHashedTopologicalTorsionFingerprintAsBitVect(mol)
# MACCS keys (166-bit structural keys)
from rdkit.Chem import MACCSkeys
fp_maccs = MACCSkeys.GenMACCSKeys(mol)Substructure Searching
SMARTS Pattern Matching
from rdkit import Chem
# Define molecule
mol = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O") # Aspirin
# SMARTS pattern for carboxylic acid
pattern = Chem.MolFromSmarts("C(=O)O")
# Check if substructure exists
has_match = mol.HasSubstructMatch(pattern)
print(f"Contains carboxylic acid: {has_match}")
# Get matching atoms
matches = mol.GetSubstructMatches(pattern)
print(f"Number of matches: {len(matches)}")
print(f"Matching atom indices: {matches}")
# Common SMARTS patterns
patterns = {
'Carboxylic acid': 'C(=O)O',
'Ester': 'C(=O)O[C,c]',
'Amide': 'C(=O)N',
'Alcohol': '[OH][C,c]',
'Primary amine': '[NH2][C,c]',
'Aromatic ring': 'c1ccccc1'
}
for name, smarts in patterns.items():
pattern = Chem.MolFromSmarts(smarts)
if mol.HasSubstructMatch(pattern):
print(f"Contains {name}")Substructure Filtering
from rdkit import Chem
# Library of molecules
smiles_list = [
"CC(=O)O", # Acetic acid
"CCO", # Ethanol
"c1ccccc1C(=O)O", # Benzoic acid
"CCCC", # Butane
]
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
# Filter molecules containing carboxylic acid
pattern = Chem.MolFromSmarts("C(=O)O")
filtered = [mol for mol in mols if mol.HasSubstructMatch(pattern)]
print(f"Molecules with carboxylic acid: {len(filtered)}/{len(mols)}")
# Filter by multiple patterns (AND logic)
pattern1 = Chem.MolFromSmarts("c1ccccc1") # Aromatic ring
pattern2 = Chem.MolFromSmarts("C(=O)O") # Carboxylic acid
aromatic_acids = [
mol for mol in mols
if mol.HasSubstructMatch(pattern1) and mol.HasSubstructMatch(pattern2)
]Replace Substructures
from rdkit import Chem
from rdkit.Chem import AllChem
# Replace carboxylic acid with ester
mol = Chem.MolFromSmiles("CC(=O)O") # Acetic acid
# Define replacement
rxn = AllChem.ReactionFromSmarts('[C:1](=O)O>>[C:1](=O)OC')
# Apply reaction
products = rxn.RunReactants((mol,))
if products:
product = products[0][0]
print(f"Product: {Chem.MolToSmiles(product)}") # Methyl acetate3D Conformer Generation
Generate 3D Coordinates
from rdkit import Chem
from rdkit.Chem import AllChem
# Create molecule
mol = Chem.MolFromSmiles("CCO")
# Add hydrogens (required for 3D)
mol = Chem.AddHs(mol)
# Generate 3D coordinates
result = AllChem.EmbedMolecule(mol, randomSeed=42)
if result == 0: # Success
print("3D coordinates generated")
else:
print("Failed to generate 3D coordinates")
# Optimize geometry with MMFF force field
AllChem.MMFFOptimizeMolecule(mol)
# Get atomic positions
conf = mol.GetConformer()
for i in range(mol.GetNumAtoms()):
pos = conf.GetAtomPosition(i)
print(f"Atom {i}: ({pos.x:.3f}, {pos.y:.3f}, {pos.z:.3f})")Multiple Conformers
from rdkit import Chem
from rdkit.Chem import AllChem
mol = Chem.MolFromSmiles("CCCC") # Butane
mol = Chem.AddHs(mol)
# Generate multiple conformers
conf_ids = AllChem.EmbedMultipleConfs(
mol,
numConfs=10,
randomSeed=42,
pruneRmsThresh=0.5 # Remove similar conformers
)
print(f"Generated {len(conf_ids)} conformers")
# Optimize each conformer
for conf_id in conf_ids:
AllChem.MMFFOptimizeMolecule(mol, confId=conf_id)
# Get energies
props = AllChem.MMFFGetMoleculeProperties(mol)
for conf_id in conf_ids:
ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=conf_id)
energy = ff.CalcEnergy()
print(f"Conformer {conf_id}: {energy:.2f} kcal/mol")Molecular Alignment
from rdkit import Chem
from rdkit.Chem import AllChem
# Reference molecule
ref_mol = Chem.MolFromSmiles("c1ccccc1C") # Toluene
ref_mol = Chem.AddHs(ref_mol)
AllChem.EmbedMolecule(ref_mol)
# Probe molecule
probe_mol = Chem.MolFromSmiles("c1ccccc1CC") # Ethylbenzene
probe_mol = Chem.AddHs(probe_mol)
AllChem.EmbedMolecule(probe_mol)
# Align probe to reference
rmsd = AllChem.AlignMol(probe_mol, ref_mol)
print(f"RMSD: {rmsd:.3f} Å")
# Get aligned coordinates
# Now probe_mol has coordinates aligned to ref_molMolecular Visualization
2D Drawings
from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
# Single molecule
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O") # Aspirin
img = Draw.MolToImage(mol, size=(300, 300))
plt.imshow(img)
plt.axis('off')
plt.show()
# Multiple molecules
mols = [Chem.MolFromSmiles(s) for s in ["CCO", "c1ccccc1", "CC(=O)O"]]
legends = ["Ethanol", "Benzene", "Acetic acid"]
img = Draw.MolsToGridImage(
mols,
molsPerRow=3,
subImgSize=(200, 200),
legends=legends
)
plt.imshow(img)
plt.axis('off')
plt.show()Highlight Substructures
from rdkit import Chem
from rdkit.Chem import Draw
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O") # Aspirin
# Highlight carboxylic acid group
pattern = Chem.MolFromSmarts("C(=O)O")
match = mol.GetSubstructMatch(pattern)
# Draw with highlighted atoms
img = Draw.MolToImage(mol, highlightAtoms=match, size=(300, 300))Save to File
from rdkit import Chem
from rdkit.Chem import Draw
mol = Chem.MolFromSmiles("CCO")
# Save as PNG
Draw.MolToFile(mol, "molecule.png", size=(300, 300))
# Save as SVG (vector graphics)
from rdkit.Chem.Draw import rdMolDraw2D
drawer = rdMolDraw2D.MolDraw2DSVG(300, 300)
drawer.DrawMolecule(mol)
drawer.FinishDrawing()
svg = drawer.GetDrawingText()
with open("molecule.svg", "w") as f:
f.write(svg)Practical Workflows
1. Virtual Screening Pipeline
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Lipinski
from rdkit import DataStructs
import pandas as pd
def screen_library(library_file, reference_smiles, similarity_threshold=0.7):
"""Screen compound library for similar, drug-like molecules."""
# Reference molecule and fingerprint
ref_mol = Chem.MolFromSmiles(reference_smiles)
ref_fp = AllChem.GetMorganFingerprintAsBitVect(ref_mol, radius=2)
# Results
hits = []
# Read library
suppl = Chem.SDMolSupplier(library_file)
for i, mol in enumerate(suppl):
if mol is None:
continue
# Step 1: Drug-likeness filter (Lipinski)
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)
if not (mw <= 500 and logp <= 5 and hbd <= 5 and hba <= 10):
continue
# Step 2: Similarity filter
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2)
similarity = DataStructs.TanimotoSimilarity(ref_fp, fp)
if similarity < similarity_threshold:
continue
# Step 3: PAINS filter (Pan-Assay Interference Compounds)
from rdkit.Chem import FilterCatalog
params = FilterCatalog.FilterCatalogParams()
params.AddCatalog(FilterCatalog.FilterCatalogParams.FilterCatalogs.PAINS)
catalog = FilterCatalog.FilterCatalog(params)
if catalog.HasMatch(mol):
continue
# Passed all filters
hits.append({
'id': i,
'smiles': Chem.MolToSmiles(mol),
'similarity': similarity,
'mw': mw,
'logp': logp
})
# Convert to DataFrame
df_hits = pd.DataFrame(hits)
df_hits = df_hits.sort_values('similarity', ascending=False)
return df_hits
# Usage
# hits = screen_library('compounds.sdf', 'CC(=O)OC1=CC=CC=C1C(=O)O')
# print(f"Found {len(hits)} hits")2. Diversity Selection
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs
import numpy as np
def select_diverse_set(smiles_list, n_select=100):
"""Select diverse subset using MaxMin algorithm."""
# Generate fingerprints
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
fps = [AllChem.GetMorganFingerprintAsBitVect(m, radius=2) for m in mols if m]
if len(fps) < n_select:
return list(range(len(fps)))
# MaxMin diversity picking
from rdkit.SimDivFilters import MaxMinPicker
def distance_function(i, j):
return 1 - DataStructs.TanimotoSimilarity(fps[i], fps[j])
picker = MaxMinPicker()
picks = picker.LazyPick(
distance_function,
len(fps),
n_select,
seed=42
)
return list(picks)
# Usage
smiles_list = ["CCO", "CCCO", "c1ccccc1", "CC(=O)O", "CCCCCCCC"]
diverse_indices = select_diverse_set(smiles_list, n_select=3)
diverse_smiles = [smiles_list[i] for i in diverse_indices]
print(f"Selected: {diverse_smiles}")3. QSAR Model Building
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import numpy as np
def build_qsar_model(smiles_list, activities):
"""Build QSAR model from SMILES and activities."""
# Generate fingerprints
fps = []
valid_activities = []
for smi, act in zip(smiles_list, activities):
mol = Chem.MolFromSmiles(smi)
if mol is not None:
# Morgan fingerprint as features
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps.append(arr)
valid_activities.append(act)
X = np.array(fps)
y = np.array(valid_activities)
# Split data
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# Train model
model = RandomForestRegressor(n_estimators=100, random_state=42)
model.fit(X_train, y_train)
# Evaluate
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)
print(f"Train R²: {r2_score(y_train, y_pred_train):.3f}")
print(f"Test R²: {r2_score(y_test, y_pred_test):.3f}")
print(f"Test RMSE: {mean_squared_error(y_test, y_pred_test, squared=False):.3f}")
return model
# Usage
# smiles = ["CCO", "CCCO", "c1ccccc1", "CC(=O)O"]
# activities = [5.2, 4.8, 6.1, 5.5] # pIC50 values
# model = build_qsar_model(smiles, activities)4. Scaffold Analysis
from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import Counter
def analyze_scaffolds(smiles_list):
"""Analyze Murcko scaffolds in compound set."""
scaffolds = []
for smi in smiles_list:
mol = Chem.MolFromSmiles(smi)
if mol is not None:
# Get Murcko scaffold
scaffold = MurckoScaffold.GetScaffoldForMol(mol)
scaffold_smi = Chem.MolToSmiles(scaffold)
scaffolds.append(scaffold_smi)
# Count scaffolds
scaffold_counts = Counter(scaffolds)
print(f"Unique scaffolds: {len(scaffold_counts)}")
print("\nTop 5 scaffolds:")
for scaffold, count in scaffold_counts.most_common(5):
print(f"{scaffold}: {count}")
return scaffold_counts
# Usage
# smiles_list = ["c1ccccc1CC", "c1ccccc1CCC", "c1ccc(O)cc1"]
# scaffolds = analyze_scaffolds(smiles_list)5. Reaction Enumeration
from rdkit import Chem
from rdkit.Chem import AllChem
def enumerate_amide_coupling(acids, amines):
"""Enumerate all possible amide products."""
# Define reaction SMARTS
rxn = AllChem.ReactionFromSmarts('[C:1](=[O:2])O.[N:3]>>[C:1](=[O:2])[N:3]')
products = []
for acid_smi in acids:
for amine_smi in amines:
acid = Chem.MolFromSmiles(acid_smi)
amine = Chem.MolFromSmiles(amine_smi)
if acid is None or amine is None:
continue
# Run reaction
products_tuple = rxn.RunReactants((acid, amine))
if products_tuple:
product = products_tuple[0][0]
Chem.SanitizeMol(product)
product_smi = Chem.MolToSmiles(product)
products.append({
'acid': acid_smi,
'amine': amine_smi,
'product': product_smi
})
return products
# Usage
acids = ["CC(=O)O", "c1ccccc1C(=O)O"]
amines = ["CCN", "c1ccccc1N"]
products = enumerate_amide_coupling(acids, amines)
print(f"Generated {len(products)} amides")Performance Optimization
Bulk Operations
from rdkit import Chem
from rdkit.Chem import AllChem
import pandas as pd
# Instead of loop
# ❌ SLOW
# fps = []
# for smi in smiles_list:
# mol = Chem.MolFromSmiles(smi)
# fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
# fps.append(fp)
# ✅ FAST: Use PandasTools for bulk operations
from rdkit.Chem import PandasTools
df = pd.DataFrame({'SMILES': smiles_list})
PandasTools.AddMoleculeColumnToFrame(df, 'SMILES', 'Molecule')
# Calculate properties in bulk
df['MW'] = df['Molecule'].apply(lambda x: Descriptors.MolWt(x) if x else None)
df['LogP'] = df['Molecule'].apply(lambda x: Descriptors.MolLogP(x) if x else None)Parallel Processing
from rdkit import Chem
from rdkit.Chem import AllChem
from multiprocessing import Pool
import pandas as pd
def process_molecule(smiles):
"""Process single molecule."""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return {
'smiles': smiles,
'mw': Descriptors.MolWt(mol),
'logp': Descriptors.MolLogP(mol),
'fp': AllChem.GetMorganFingerprintAsBitVect(mol, 2)
}
def process_library_parallel(smiles_list, n_jobs=4):
"""Process library in parallel."""
with Pool(n_jobs) as pool:
results = pool.map(process_molecule, smiles_list)
# Filter None results
results = [r for r in results if r is not None]
return pd.DataFrame(results)
# Usage
# df = process_library_parallel(large_smiles_list, n_jobs=8)Caching Calculations
from functools import lru_cache
from rdkit import Chem
from rdkit.Chem import Descriptors
@lru_cache(maxsize=10000)
def get_mol_properties(smiles):
"""Calculate properties with caching."""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return {
'mw': Descriptors.MolWt(mol),
'logp': Descriptors.MolLogP(mol),
'tpsa': Descriptors.TPSA(mol)
}
# Repeated calls will use cache
props1 = get_mol_properties("CCO") # Calculated
props2 = get_mol_properties("CCO") # Cached (instant)Common Pitfalls and Solutions
The "Invalid SMILES" Problem
Not all SMILES strings are valid.
# ❌ Problem: Assuming all SMILES are valid
smiles_list = ["CCO", "INVALID", "c1ccccc1"]
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
# Contains None!
# ✅ Solution: Filter invalid molecules
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
valid_mols = [m for m in mols if m is not None]
# ✅ Better: Track which failed
results = []
for smi in smiles_list:
mol = Chem.MolFromSmiles(smi)
if mol is not None:
results.append({'smiles': smi, 'mol': mol, 'valid': True})
else:
results.append({'smiles': smi, 'mol': None, 'valid': False})The "Stereochemistry Loss" Problem
SMILES generation can lose stereochemistry if not careful.
from rdkit import Chem
# Molecule with stereochemistry
chiral_smiles = "C[C@H](O)CC" # (S)-2-butanol
mol = Chem.MolFromSmiles(chiral_smiles)
# ❌ BAD: Lose stereochemistry
non_iso = Chem.MolToSmiles(mol, isomericSmiles=False)
print(non_iso) # "CC(O)CC" - lost chirality!
# ✅ GOOD: Preserve stereochemistry
iso = Chem.MolToSmiles(mol, isomericSmiles=True)
print(iso) # "C[C@H](O)CC" - preserved!The "3D Without Hydrogens" Problem
3D operations require explicit hydrogens.
from rdkit import Chem
from rdkit.Chem import AllChem
mol = Chem.MolFromSmiles("CCO")
# ❌ BAD: Generate 3D without hydrogens
result = AllChem.EmbedMolecule(mol)
# Poor quality or failure
# ✅ GOOD: Add hydrogens first
mol_h = Chem.AddHs(mol)
result = AllChem.EmbedMolecule(mol_h)
AllChem.MMFFOptimizeMolecule(mol_h)The "Fingerprint Type Mismatch" Problem
Comparing different fingerprint types gives meaningless results.
from rdkit import Chem
from rdkit.Chem import AllChem, MACCSkeys
from rdkit import DataStructs
mol1 = Chem.MolFromSmiles("CCO")
mol2 = Chem.MolFromSmiles("CCCO")
# ❌ BAD: Comparing different fingerprint types
fp1_morgan = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2_maccs = MACCSkeys.GenMACCSKeys(mol2)
# This will error or give nonsense!
# similarity = DataStructs.TanimotoSimilarity(fp1_morgan, fp2_maccs)
# ✅ GOOD: Use same fingerprint type
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, 2)
similarity = DataStructs.TanimotoSimilarity(fp1, fp2)The "Memory Explosion" Problem
Processing millions of molecules can exhaust memory.
# ❌ BAD: Load entire library into memory
suppl = Chem.SDMolSupplier('huge_library.sdf')
mols = [mol for mol in suppl] # Out of memory!
# ✅ GOOD: Process in batches
def process_in_batches(sdf_file, batch_size=10000):
suppl = Chem.SDMolSupplier(sdf_file)
batch = []
for mol in suppl:
if mol is not None:
batch.append(mol)
if len(batch) >= batch_size:
# Process batch
yield batch
batch = []
# Process remaining
if batch:
yield batch
# Usage
for batch in process_in_batches('huge_library.sdf'):
# Process each batch
passRDKit is the cornerstone of computational drug discovery and cheminformatics. Its comprehensive toolkit for molecular manipulation, descriptor calculation, and similarity searching makes it indispensable for pharmaceutical research, virtual screening, and chemical data analysis. Master RDKit, and you'll have the power to computationally explore vast chemical spaces and accelerate drug discovery.