# Bio Applied Qsar Modeling

> Build QSAR classifiers from ChEMBL IC50 data with RDKit Morgan fingerprints, Random Forest, scaffold splits, and k-NN applicability domain. Use when predicting activity/pIC50 from SMILES or building structure-activity relationship models.

- Skill: `pavel-kravchenko/bio-applied-qsar-modeling` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-qsar-modeling`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-qsar-modeling/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-qsar-modeling

---


# QSAR Modeling

## When to Use

- Predicting bioactivity (active/inactive, pIC50, Ki) for compounds from SMILES using a trained model.
- Building a structure-activity relationship (SAR) model from ChEMBL/PubChem bioactivity assay data.
- Scoring a virtual screening library and flagging predictions that fall outside the training chemical space.
- Comparing scaffold-based vs. random train/test splits to get a realistic estimate of generalization.
- Prioritizing synthesis candidates by combining a QSAR activity score with an applicability domain (AD) flag.

## Version Compatibility

- Python ≥3.10, RDKit ≥2023.09, scikit-learn ≥1.4, pandas ≥2.0, chembl_webresource_client ≥0.10.

## Prerequisites

- `pip install rdkit chembl_webresource_client scikit-learn pandas numpy`
- Familiarity with SMILES and basic ML classification metrics (AUC-ROC, MCC).
- Related skills: `bio-chemoinformatics-molecular-descriptors`, `bio-chemoinformatics-similarity-searching`, `bio-chemoinformatics-virtual-screening`.

## 1. Data Retrieval from ChEMBL

**Goal:** pull bioactivity measurements (IC50) for a target protein.
**Approach:** query `new_client.activity` filtered by `target_chembl_id` and `standard_type`, restricted to the fields needed downstream.

```python
from chembl_webresource_client.new_client import new_client
import pandas as pd

def fetch_chembl_ic50(target_chembl_id: str = "CHEMBL203") -> pd.DataFrame:
    """Fetch IC50 bioactivity records for a ChEMBL target (default: EGFR).

    Returns a DataFrame with molecule_chembl_id, standard_value (nM),
    standard_units, and canonical_smiles.
    """
    activity = new_client.activity
    records = activity.filter(
        target_chembl_id=target_chembl_id,
        standard_type="IC50",
    ).only(
        ["molecule_chembl_id", "standard_value", "standard_units", "canonical_smiles"]
    )
    return pd.DataFrame(list(records))
```

## 2. Data Curation

**Goal:** turn raw IC50 values into a clean binary activity label.
**Approach:** drop missing/invalid rows, convert nM IC50 to pIC50 (`-log10(IC50_M)`), collapse duplicate compounds to their median pIC50, and threshold at pIC50 ≥ 6 (IC50 ≤ 1 µM) for "active".

```python
import numpy as np

def curate_activity(df: pd.DataFrame, active_threshold: float = 6.0) -> pd.DataFrame:
    """Convert IC50 (nM) to pIC50, dedupe by median, and label active/inactive."""
    df = df.dropna(subset=["standard_value", "canonical_smiles"]).copy()
    df["standard_value"] = pd.to_numeric(df["standard_value"], errors="coerce")
    df = df[df["standard_value"] > 0]
    df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)

    # Keep median pIC50 per compound to reduce assay noise
    df = (
        df.groupby(["molecule_chembl_id", "canonical_smiles"], as_index=False)["pIC50"]
        .median()
    )
    df["active"] = (df["pIC50"] >= active_threshold).astype(int)
    return df
```

## 3. Feature Generation

**Goal:** turn each SMILES into a numeric feature vector.
**Approach:** compute a 2048-bit Morgan fingerprint (radius 2 ≈ ECFP4) plus the full set of RDKit 2D descriptors; concatenate them. Invalid SMILES are skipped and their rows dropped.

```python
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors

_DESC_NAMES = [name for name, _ in Descriptors.descList]
_CALC = MoleculeDescriptors.MolecularDescriptorCalculator(_DESC_NAMES)

def featurize(smiles_list, n_bits: int = 2048, radius: int = 2):
    """Compute Morgan fingerprint + RDKit 2D descriptors for each SMILES.

    Returns (X, valid_mask) where X is an (n_valid, n_bits + n_descriptors)
    float array and valid_mask marks which input rows parsed successfully.
    """
    rows, valid_mask = [], []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            valid_mask.append(False)
            continue
        fp = list(AllChem.GetMorganFingerprintAsBitVect(mol, radius, n_bits))
        desc = list(_CALC.CalcDescriptors(mol))
        rows.append(fp + desc)
        valid_mask.append(True)
    X = np.array(rows, dtype=float)
    return X, np.array(valid_mask)
```

## 4. Scaffold Split, Training, and Evaluation

**Goal:** train a classifier and evaluate it honestly (scaffold split avoids leaking near-duplicate analogs between train/test).
**Approach:** bin molecules by Murcko scaffold, assign whole scaffold groups to train/test, then fit `RandomForestClassifier` and score with AUC-ROC and MCC.

```python
from collections import defaultdict
from rdkit.Chem.Scaffolds import MurckoScaffold
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, matthews_corrcoef

def scaffold_split(smiles_list, test_frac: float = 0.2, seed: int = 42):
    """Split molecule indices into train/test by Bemis-Murcko scaffold.

    Whole scaffold groups go entirely to train or test, so structurally
    near-identical analogs never straddle the split (unlike a random split).
    """
    scaffold_to_idx = defaultdict(list)
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        scaffold = MurckoScaffold.MurckoScaffoldSmiles(mol=mol) if mol else ""
        scaffold_to_idx[scaffold].append(i)

    rng = np.random.RandomState(seed)
    groups = list(scaffold_to_idx.values())
    rng.shuffle(groups)

    n_total = len(smiles_list)
    test_idx, train_idx = [], []
    for group in groups:
        if len(test_idx) < test_frac * n_total:
            test_idx.extend(group)
        else:
            train_idx.extend(group)
    return np.array(train_idx), np.array(test_idx)


def train_and_evaluate(X, y, train_idx, test_idx, seed: int = 42):
    """Fit a Random Forest QSAR classifier and report AUC-ROC and MCC."""
    X_train, X_test = X[train_idx], X[test_idx]
    y_train, y_test = y[train_idx], y[test_idx]

    rf = RandomForestClassifier(n_estimators=500, random_state=seed, n_jobs=-1)
    rf.fit(X_train, y_train)

    y_pred = rf.predict(X_test)
    y_prob = rf.predict_proba(X_test)[:, 1]
    metrics = {
        "auc_roc": roc_auc_score(y_test, y_prob),
        "mcc": matthews_corrcoef(y_test, y_pred),
    }
    return rf, metrics
```

## 5. Applicability Domain (AD)

**Goal:** flag test/new compounds too dissimilar from the training set for the model's predictions to be trusted.
**Approach:** compute mean distance to the k=5 nearest training-set neighbors; anything beyond `mean + 3*std` of the training self-distances is out of domain.

```python
from sklearn.neighbors import NearestNeighbors

def applicability_domain(X_train, X_query, n_neighbors: int = 5, z: float = 3.0):
    """Flag query compounds outside the k-NN applicability domain of X_train.

    Returns (in_ad, avg_distances) — in_ad is a boolean array over X_query.
    """
    knn = NearestNeighbors(n_neighbors=n_neighbors)
    knn.fit(X_train)

    # Threshold from training set's own leave-one-out neighbor distances
    train_dist, _ = knn.kneighbors(X_train, n_neighbors=n_neighbors + 1)
    train_avg = train_dist[:, 1:].mean(axis=1)  # drop self-match at index 0
    threshold = train_avg.mean() + z * train_avg.std()

    query_dist, _ = knn.kneighbors(X_query, n_neighbors=n_neighbors)
    query_avg = query_dist.mean(axis=1)
    return query_avg <= threshold, query_avg
```

```python
def demo():
    """Self-check: featurize toy SMILES, scaffold-split, train, AD-check."""
    smiles = [
        "CCOc1ccc(cc1)C(=O)Nc1ccccc1", "CCOc1ccc(cc1)C(=O)Nc1ccccn1",
        "c1ccccc1C(=O)O", "c1ccccc1C(=O)N", "CCN(CC)CCOC(=O)c1ccccc1N",
        "CCN(CC)CCOC(=O)c1ccccc1", "COc1ccccc1", "COc1ccccc1O",
    ]
    y = np.array([1, 1, 0, 0, 1, 1, 0, 0])

    X, mask = featurize(smiles)
    assert X.shape[0] == len(smiles) and mask.all()

    train_idx, test_idx = scaffold_split(smiles, test_frac=0.25, seed=0)
    assert len(train_idx) + len(test_idx) == len(smiles)

    rf, metrics = train_and_evaluate(X, y, train_idx, test_idx)
    assert 0.0 <= metrics["auc_roc"] <= 1.0

    in_ad, dist = applicability_domain(X[train_idx], X[test_idx])
    assert in_ad.shape[0] == len(test_idx)
    print("demo OK:", metrics, "in_ad:", in_ad.tolist())


if __name__ == "__main__":
    demo()
```

## Pitfalls

- **SMILES canonicalization**: different SMILES can represent the same molecule; canonicalize (`Chem.MolToSmiles`) before deduplicating.
- **Stereochemistry**: default Morgan fingerprints ignore chirality unless `useChirality=True` is passed, which can merge enantiomers with different bioactivity.
- **Descriptor scaling**: RDKit 2D descriptors span wildly different ranges (MW vs. fraction Csp3); standardize (`StandardScaler`) before distance-based steps like AD or SVM/kNN.
- **Random vs. scaffold split**: a random split lets near-duplicate analogs leak between train/test, inflating AUC/MCC — always report scaffold-split metrics for realistic generalization estimates.
- **Applicability domain threshold**: the `mean + 3*std` k-NN rule is a heuristic, not a guarantee; tune `n_neighbors`/`z` against a held-out validation set for your chemical series.
- **Class imbalance**: bioactivity datasets are often skewed toward actives or inactives depending on the assay; stratify splits and prefer MCC/AUC-PR over raw accuracy.

## See Also

- `bio-chemoinformatics-molecular-descriptors`
- `bio-chemoinformatics-similarity-searching`
- `bio-chemoinformatics-virtual-screening`
- `bio-chemoinformatics-admet-prediction`