bio-applied-qsar-modeling

name	bio-applied-qsar-modeling
description	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.
tool_type	python
primary_tool	RDKit

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.

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".

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.

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.

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=, random_state=seed, n_jobs=-)
    rf.fit(X_train, y_train)

    y_pred = rf.predict(X_test)
    y_prob = rf.predict_proba(X_test)[:, ]
    metrics = {
        : roc_auc_score(y_test, y_prob),
        : matthews_corrcoef(y_test, y_pred),
    }
     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.

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

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