# Cheminformatics Drug Discovery

> Parse SMILES/SDF with RDKit, build ChEMBL QSAR models, dock ligands with AutoDock Vina, predict ADMET/Lipinski properties. Use for compound screening, Tanimoto similarity, Rule of Five filtering.

- Skill: `pavel-kravchenko/cheminformatics-drug-discovery` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/cheminformatics-drug-discovery`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/cheminformatics-drug-discovery/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Product & Planning
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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.

```python
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.

```python
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.

```bash
# 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
```

```python
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