Bioactivity Classification

Bioactivity classification is the task of predicting whether a molecule will be active or inactive against a specific biological target, typically a protein receptor, enzyme, or ion channel. It is one of the most common ML tasks in early-stage drug discovery, used to prioritize compounds for synthesis and testing, and to train predictive models that guide medicinal chemistry.

This guide walks through the complete pipeline: retrieving and cleaning ChEMBL data, defining activity thresholds, choosing data splits, training and evaluating models with ECFP4 fingerprints and XGBoost, and understanding the unique challenges of bioactivity prediction including activity cliffs and analog bias.

What Bioactivity Classification Is

In biochemical assays, compound activity against a target is measured as a potency value such as IC50 (half-maximal inhibitory concentration), Ki (inhibition constant), or EC50 (half-maximal effective concentration). These continuous measurements are converted to binary labels for classification models using an activity cutoff, typically 1 or 10 micromolar (uM).

The classification problem is then: given the molecular structure of a compound, predict whether it will have IC50 below the threshold (active) or above it (inactive). This is a hit identification problem: the model is used to score large virtual libraries and identify the most promising candidates for experimental testing.

A compound with IC50 = 0.1 nM is clearly active; IC50 = 100 uM is clearly inactive. The gray zone around the cutoff (e.g., 0.5-20 uM) creates label noise that no model can fully overcome.

Data Sources

For most target-specific bioactivity classification work, ChEMBL is the primary source. Its Python client library allows programmatic data retrieval by target protein identifier (UniProt accession or ChEMBL target ID).

Activity Cutoffs: IC50, Ki, EC50, and pIC50

The choice of activity cutoff significantly affects the positive-to-negative ratio and the model's utility:

• IC50: Concentration at which 50% inhibition is observed. Standard for enzyme and receptor assays. Common cutoffs: 1 uM (nanomolar range = active, micromolar range = inactive) or 10 uM for broader datasets.
• Ki: Thermodynamic binding constant. Directly related to binding affinity, preferred over IC50 when available. Cutoff typically at 1-10 uM.
• EC50: Used for functional assays (agonist potency). Context-dependent cutoff.
• pIC50: Negative log10 of IC50 in molar (so IC50 = 1 nM gives pIC50 = 9). Converting to pIC50 scale linearizes the distribution and is preferred for regression tasks. For classification, a common cutoff is pIC50 >= 6 (i.e., IC50 <= 1 uM = active).

When building a dataset from ChEMBL, use only a single activity type (do not mix IC50 and Ki without unit conversion), standardize to the same assay format where possible, and document your cutoff clearly.

Featurization: ECFP4 and Graph Representations

The two main featurization approaches for bioactivity classification:

• ECFP4 fingerprints (2048 bits): Fast, interpretable, strong baseline. Captures the circular atom neighborhoods that determine binding complementarity with the target. For typical ChEMBL target datasets (500-10,000 compounds), ECFP4 + XGBoost is highly competitive.
• Molecular graphs: Atoms as nodes, bonds as edges. Input for Graph Neural Networks. Allows end-to-end learning of task-specific representations. Benefits emerge with 5,000+ compounds or when the SAR has complex graph-level patterns that fingerprints cannot capture.

Data Splits in Drug Discovery: Scaffold Split vs. Random Split

The choice of data split has a dramatic effect on reported performance and its real-world relevance:

Always use scaffold split for published benchmarks. When comparing your model to the literature, use the same split. Random splits can make a simple 1-NN look like it achieves 0.99 AUROC on bioactivity data, which is meaningless.

Code Example: ChEMBL Data and XGBoost Classifier

import pandas as pd
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from chembl_webresource_client.new_client import new_client
from xgboost import XGBClassifier
from sklearn.metrics import roc_auc_score, average_precision_score, classification_report

# -------------------------------------------------------
# 1. Retrieve ChEMBL data for a specific target
# -------------------------------------------------------
# Example: EGFR kinase (ChEMBL ID: CHEMBL203)
activity_client = new_client.activity
target_id = 'CHEMBL203'  # EGFR

print(f"Fetching activities for {target_id}...")
activities = activity_client.filter(
    target_chembl_id=target_id,
    standard_type='IC50',
    relation='=',
    assay_type='B'  # Binding assay
).only(['molecule_chembl_id', 'canonical_smiles', 'standard_value', 'standard_units'])

df = pd.DataFrame(list(activities))
print(f"Retrieved {len(df)} activity records")

# -------------------------------------------------------
# 2. Data cleaning and label assignment
# -------------------------------------------------------
# Keep only nM measurements
df = df[df['standard_units'] == 'nM'].copy()

# Convert to numeric, drop missing
df['ic50_nm'] = pd.to_numeric(df['standard_value'], errors='coerce')
df = df.dropna(subset=['ic50_nm', 'canonical_smiles'])

# Remove duplicates: keep median IC50 per compound
df = df.groupby('canonical_smiles').agg({'ic50_nm': 'median'}).reset_index()

# Binary label: active if IC50 < 1000 nM (1 uM)
CUTOFF_NM = 1000.0
df['active'] = (df['ic50_nm'] < CUTOFF_NM).astype(int)

print(f"Dataset: {len(df)} compounds")
print(f"Active: {df['active'].sum()} ({df['active'].mean():.1%})")

# -------------------------------------------------------
# 3. Scaffold split using RDKit
# -------------------------------------------------------
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import defaultdict

def get_scaffold(smiles: str) -> str:
    """Return Bemis-Murcko scaffold SMILES for a molecule."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return ''
    scaffold = MurckoScaffold.GetScaffoldForMol(mol)
    return Chem.MolToSmiles(scaffold)

df['scaffold'] = df['canonical_smiles'].apply(get_scaffold)

# Group compounds by scaffold
scaffold_groups = defaultdict(list)
for idx, row in df.iterrows():
    scaffold_groups[row['scaffold']].append(idx)

# Sort scaffolds by size (descending) and assign train/test
scaffolds = sorted(scaffold_groups.items(), key=lambda x: len(x[1]), reverse=True)
train_idx, test_idx = [], []
test_cutoff = int(0.8 * len(df))

total_assigned = 0
for scaffold, indices in scaffolds:
    if total_assigned < test_cutoff:
        train_idx.extend(indices)
    else:
        test_idx.extend(indices)
    total_assigned += len(indices)

train_df = df.loc[train_idx].reset_index(drop=True)
test_df  = df.loc[test_idx].reset_index(drop=True)

print(f"Train: {len(train_df)} compounds, Test: {len(test_df)} compounds")

# -------------------------------------------------------
# 4. Featurization: ECFP4 fingerprints
# -------------------------------------------------------
def compute_ecfp4(smiles_list: list, n_bits: int = 2048) -> np.ndarray:
    X = np.zeros((len(smiles_list), n_bits), dtype=np.uint8)
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        if mol is not None:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=n_bits)
            fp.ToNumpyArray(X[i])
    return X

X_train = compute_ecfp4(train_df['canonical_smiles'].tolist())
X_test  = compute_ecfp4(test_df['canonical_smiles'].tolist())
y_train = train_df['active'].values
y_test  = test_df['active'].values

# -------------------------------------------------------
# 5. XGBoost classifier
# -------------------------------------------------------
pos_weight = (y_train == 0).sum() / max((y_train == 1).sum(), 1)

clf = XGBClassifier(
    n_estimators=400,
    max_depth=7,
    learning_rate=0.03,
    subsample=0.8,
    colsample_bytree=0.8,
    scale_pos_weight=pos_weight,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42,
    n_jobs=-1
)
clf.fit(X_train, y_train)

y_pred_proba = clf.predict_proba(X_test)[:, 1]
y_pred_label = (y_pred_proba >= 0.5).astype(int)

auroc = roc_auc_score(y_test, y_pred_proba)
auprc = average_precision_score(y_test, y_pred_proba)

print(f"\nResults on scaffold test set:")
print(f"  AUROC: {auroc:.4f}")
print(f"  AUPRC: {auprc:.4f}")
print(classification_report(y_test, y_pred_label, target_names=['Inactive', 'Active']))

Evaluation and Common Pitfalls

Evaluation Metrics

• AUROC: Primary metric for ranking ability. Reports the probability that a randomly chosen active is ranked higher than a randomly chosen inactive. Target above 0.80 on scaffold split for a reliable model.
• AUPRC: Especially important when the active class is a small minority (as in large screening datasets with 1-5% hit rate).
• Enrichment Factor (EF): The ratio of actives in the top X% of predictions to the baseline active rate. EF at 1% or 5% is the primary metric for virtual screening applications.

Activity Cliffs

An activity cliff occurs when two structurally very similar compounds (high Tanimoto similarity) have dramatically different activity values. For example, adding a single methyl group can shift IC50 from 0.1 nM to 10,000 nM due to a steric clash in the binding pocket. Activity cliffs are the primary limitation of fingerprint-based models: two nearly identical fingerprints map to very different labels, and the model cannot learn this from structure alone. GNNs handle some activity cliffs better because they can learn more nuanced local structural features, but neither approach handles all activity cliffs reliably.

Analog Bias

If the training set contains many analogs of the active compounds (compounds that differ by only one or two atoms from a known active), the model learns a narrow region of chemical space very well but fails to generalize to structurally novel actives. Scaffold splits partially address this, but cluster-based splits with explicit similarity thresholds are more stringent.

Assay Interference Compounds

Pan-assay interference compounds (PAINS) are molecules that give false-positive readings in many assay types due to assay interference (fluorescence, aggregation, redox activity) rather than genuine binding. Including PAINS in your training data as "actives" trains a model to predict interference rather than binding. Filter PAINS using RDKit's FilterCatalog with the PAINS A/B/C catalog before training.

from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams

def is_pains(smiles: str) -> bool:
    """Return True if molecule matches any PAINS filter."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return False
    params = FilterCatalogParams()
    params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
    catalog = FilterCatalog(params)
    return catalog.HasMatch(mol)

# Filter PAINS from dataset before training
df['is_pains'] = df['canonical_smiles'].apply(is_pains)
df_clean = df[~df['is_pains']].copy()
print(f"Removed {df['is_pains'].sum()} PAINS compounds")

Further Reading

• ChEMBL Database
• ChEMBL Python Web Resource Client
• Wu et al. (2018): MoleculeNet Benchmark
• Baell and Holloway (2010): PAINS Filters
• DeepChem for molecular ML

Try it live

Explore a related model hands-on: open Molecule Screening Engine — Similarity Search Live and change its parameters while it runs. Everything stays in your browser.