Cheminformatics Drug Discovery

Cheminformatics and Drug Discovery

When to Use

• Parsing/validating SMILES, SDF, or PDB ligand files
• Computing molecular descriptors (MW, LogP, TPSA) or Morgan/MACCS fingerprints, or running Tanimoto similarity search
• Building a ChEMBL-derived QSAR classifier for IC50/pIC50 bioactivity prediction
• Running structure-based virtual screening or AutoDock Vina docking
• Predicting ADMET properties or filtering candidates with Lipinski's Rule of Five
• Training a graph neural network (PyTorch Geometric) for molecular property prediction

Version Compatibility

RDKit ≥2024.03, Python ≥3.10, chembl_webresource_client ≥0.10, scikit-learn ≥1.4, pandas ≥2.0, AutoDock Vina 1.2.x, Meeko ≥0.5 (ligand prep), torch-geometric ≥2.5 (GNN section, optional).

Prerequisites

pip install rdkit chembl_webresource_client scikit-learn pandas numpy. For docking: AutoDock Vina binary, Meeko's mk_prepare_ligand.py, and pdbfixer for receptor prep. For GNN work: torch + torch_geometric. Assumes familiarity with SMILES notation and basic ML evaluation (train/test split, ROC-AUC).

Molecular Descriptors and Lipinski Filtering

Goal: parse a compound library and flag Rule-of-Five (drug-likeness) violations. Approach: validate each SMILES with RDKit, compute MW/LogP/HBD/HBA/TPSA, and apply the Ro5 cutoffs.

from rdkit import Chem
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd

def lipinski_filter(smiles_list):
    """Return a DataFrame with Ro5 descriptors and pass/fail flag per SMILES."""
    results = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            results.append({'smiles': smi, 'valid': False})
            continue
        props = {
            'smiles': smi,
            'valid': True,
            'MW': Descriptors.MolWt(mol),
            'LogP': Descriptors.MolLogP(mol),
            'HBD': rdMolDescriptors.CalcNumHBD(mol),
            'HBA': rdMolDescriptors.CalcNumHBA(mol),
            'TPSA': Descriptors.TPSA(mol),
        }
        props['ro5_pass'] = (
            props['MW'] <= 500 and props['LogP'] <= 5 and
            props['HBD'] <= 5 and props['HBA'] <= 10
        )
        results.append(props)
    return pd.DataFrame(results)

df = lipinski_filter(['CC(=O)Oc1ccccc1C(=O)O', 'Cn1cnc2c1c(=O)n(C)c(=O)n2C'])  # aspirin, caffeine
print(f"Ro5 pass rate: {df['ro5_pass'].mean():.1%}")

Fingerprints and Similarity Search

Goal: rank a compound library by similarity to a query molecule. Approach: compute ECFP4 (Morgan, radius=2) bit vectors and bulk Tanimoto similarity.

from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem

def similarity_search(query_smiles, library_smiles, top_n=10):
    """Rank library_smiles by Tanimoto similarity (ECFP4) to query_smiles."""
    query_mol = Chem.MolFromSmiles(query_smiles)
    query_fp = AllChem.GetMorganFingerprintAsBitVect(query_mol, radius=2, nBits=2048)

    library_mols = [Chem.MolFromSmiles(s) for s in library_smiles]
    valid = [(s, m) for s, m in zip(library_smiles, library_mols) if m is not None]
    library_fps = [AllChem.GetMorganFingerprintAsBitVect(m, 2, 2048) for _, m in valid]

    sims = DataStructs.BulkTanimotoSimilarity(query_fp, library_fps)
    ranked = sorted(zip([s for s, _ in valid], sims), key=lambda x: -x[1])
    return ranked[:top_n]

hits = similarity_search('CC(=O)Oc1ccccc1C(=O)O', library_smiles)

ChEMBL QSAR Pipeline

Goal: train a bioactivity classifier (active/inactive) from ChEMBL IC50 data. Approach: pull assay data for a target, convert IC50 to pIC50, featurize with Morgan fingerprints, and train a Random Forest with a scaffold split to avoid leakage.

from chembl_webresource_client.new_client import new_client
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Scaffolds import MurckoScaffold
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GroupShuffleSplit
from sklearn.metrics import roc_auc_score
import pandas as pd
import numpy as np

def fetch_chembl_qsar_data(target_chembl_id, activity_threshold_pic50=6.0):
    """Fetch IC50 activities for a ChEMBL target and label active/inactive by pIC50."""
    activity = new_client.activity
    data = activity.filter(
        target_chembl_id=target_chembl_id, standard_type='IC50'
    ).only(['molecule_chembl_id', 'standard_value', 'canonical_smiles'])

    df = pd.DataFrame(list(data)).dropna(subset=['standard_value', 'canonical_smiles'])
    df = df[df['standard_value'].astype(float) > 0]
    df['pIC50'] = -np.log10(df['standard_value'].astype(float) * 1e-9)
    df['active'] = (df['pIC50'] >= activity_threshold_pic50).astype(int)
    return df

df = fetch_chembl_qsar_data('CHEMBL203')  # EGFR
X = np.array([list(AllChem.GetMorganFingerprintAsBitVect(Chem.MolFromSmiles(s), 2, 2048))
              for s in df['canonical_smiles']])
y = df['active'].values

scaffolds = [MurckoScaffold.MurckoScaffoldSmiles(mol=Chem.MolFromSmiles(s), includeChirality=False)
             for s in df['canonical_smiles']]
groups = np.array([hash(s) for s in scaffolds])
train_idx, test_idx = next(GroupShuffleSplit(test_size=0.2, random_state=42).split(X, y, groups))

rf = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)
rf.fit(X[train_idx], y[train_idx])
print(f"AUC: {roc_auc_score(y[test_idx], rf.predict_proba(X[test_idx])[:, 1]):.3f}")

Docking with AutoDock Vina

Goal: dock a ligand into a prepared receptor and rank binding poses. Approach: strip waters/add hydrogens with pdbfixer, convert ligand to PDBQT with Meeko, run Vina over a defined search box, then parse the log.

# Example target: Abl kinase (PDB 1IEP), ligand: imatinib
wget -q https://files.rcsb.org/download/1IEP.pdb
pdbfixer 1IEP.pdb --output 1IEP_fixed.pdb --add-hydrogens --remove-heterogens
mk_prepare_ligand.py -i imatinib.sdf -o imatinib.pdbqt

vina --receptor 1IEP_fixed.pdbqt --ligand imatinib.pdbqt \
    --center_x 22.5 --center_y 5.0 --center_z 18.0 \
    --size_x 20 --size_y 20 --size_z 20 \
    --exhaustiveness 8 --out imatinib_docked.pdbqt --log docking.log

import re
import pandas as pd

def parse_vina_log(log_path):
    """Parse an AutoDock Vina log into a DataFrame of pose affinities and RMSDs."""
    rows = []
    with open(log_path) as f:
        for line in f:
            m = re.match(r'\s+(\d+)\s+([-\d.]+)\s+([\d.]+)\s+([\d.]+)', line)
            if m:
                rows.append({
                    'mode': int(m.group(1)),
                    'affinity_kcal_mol': float(m.group(2)),
                    'rmsd_lb': float(m.group(3)),
                    'rmsd_ub': float(m.group(4)),
                })
    return pd.DataFrame(rows)

scores = parse_vina_log('docking.log')
print(f"Best pose: {scores.iloc[0]['affinity_kcal_mol']} kcal/mol")

Lipinski's Rule of Five

MW ≤ 500 Da, LogP ≤ 5, H-bond donors ≤ 5, H-bond acceptors ≤ 10, TPSA ≤ 140 Ų.

Pitfalls

• SMILES validation — always check Chem.MolFromSmiles() is not None before computing anything; malformed SMILES silently propagate as None otherwise.
• Scaffold split, not random split — random train/test splits leak near-duplicate scaffolds between sets and inflate QSAR AUC; use Murcko-scaffold grouping.
• pIC50 vs IC50 — convert nanomolar IC50 to pIC50 = -log10(IC50_nM × 1e-9) before regressing; mixing units silently breaks thresholds like pIC50 ≥ 6.
• Vina score is not free energy — it is a fast approximate scoring function; always confirm top hits with MD, MM-GBSA, or an experimental assay before committing resources.
• Filter before docking — apply Lipinski/Veber/PAINS filters first; docking every library member wastes compute on compounds that would fail ADMET anyway.
• Fingerprint bit collisions — Morgan fingerprints with too few bits (e.g. 512) cause hash collisions that hurt similarity search on large libraries; 2048 bits is the common default.

See Also

• bio-applied-docking — deeper AutoDock Vina / Glide docking workflows
• bio-applied-virtual-screening — library triage and composite scoring pipelines
• bio-applied-molecular-gnn — graph neural networks for molecular property prediction
• structural-bioinformatics — protein structure parsing and binding-site analysis