Virtual Screening with ML

Virtual screening is the computational counterpart to high-throughput screening (HTS). Instead of physically testing millions of compounds in a robot-driven assay, you computationally score a large chemical library and identify the most promising subset for actual testing. Modern virtual screening campaigns combine rule-based filters, structure-based docking, and ML-based scoring in a hierarchical pipeline that dramatically reduces the cost and time of hit identification.

This guide covers the conceptual framework of virtual screening, traditional vs. ML-based approaches, a practical screening pipeline using RDKit and Python, evaluation metrics, and the common pitfalls that inflate reported enrichment rates.

What Virtual Screening Is in Drug Discovery Context

Drug discovery begins with a disease target (usually a protein) and the question: what small molecule might bind to this target and modulate its function? The universe of drug-like small molecules is estimated at 10^33 compounds. Even the largest virtual libraries contain 10^9 to 10^12 molecules. Virtual screening reduces this space to a manageable shortlist (typically 100-1000 compounds) for experimental testing.

A screening campaign typically follows three phases:

1. Hit identification: From billions down to thousands using property filters and fast computational methods
2. Hit validation: Experimental testing of the shortlist to confirm binding (dose-response curves, orthogonal assays)
3. Hit-to-lead: Medicinal chemistry optimization of confirmed hits into more potent, selective, drug-like lead compounds

Virtual screening is used in phase 1. Its quality is measured by how efficiently it enriches the shortlist with confirmed binders relative to random selection.

Traditional vs. ML-Based Screening

In practice, modern virtual screening campaigns use a hierarchical funnel combining multiple methods: property filters first (fastest, cheapest), then ML scoring to reduce to 10,000-100,000 compounds, then docking to reduce to 500-5,000, then manual review and experimental testing of the final 100-500.

The Screening Pipeline

Step 1: Library Preparation

A commercial or virtual screening library needs preparation before any scoring:

• SMILES standardization: Canonicalize SMILES, neutralize charges where appropriate, remove salts
• Deduplication: Remove duplicate structures (same canonical SMILES)
• Tautomer enumeration: Optionally enumerate common tautomers
• 3D conformer generation: Required for docking and 3D-based methods (RDKit ETKDG algorithm)

Step 2: Property Filtering

Remove compounds unlikely to be drug-like or likely to cause assay interference. Common filters:

• Lipinski Rule of Five: MW below 500, LogP below 5, HBD 5 or fewer, HBA 10 or fewer. Predicts oral bioavailability.
• Veber Rule: Rotatable bonds 10 or fewer, TPSA 140 or below. Additional oral bioavailability criteria.
• PAINS: Remove pan-assay interference compounds (aggregators, fluorescent compounds, redox cyclers).
• Reactive groups: Remove compounds with highly reactive functional groups (Michael acceptors, epoxides, aldehydes) unless specifically seeking covalent binders.

Step 3: ML-Based Scoring (QSAR)

Train a QSAR model on known actives and decoys for the target, then score the filtered library. The top-ranked compounds proceed to docking.

Step 4: Docking and Rescoring

Docking programs (AutoDock Vina, Glide, GOLD) predict the binding pose of a ligand in a protein pocket and estimate binding affinity. The docking score provides an orthogonal view to QSAR. ML-based rescoring functions trained on experimental binding data (EquiBind, Smina, GNINA) can improve enrichment over classical docking scores.

Tools

Code Example: Library Filtering and ML Scoring with RDKit

import pandas as pd
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams
from xgboost import XGBClassifier
from sklearn.metrics import roc_auc_score
import logging

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

# -------------------------------------------------------
# 1. Property filters
# -------------------------------------------------------
def lipinski_pass(mol) -> bool:
    """Return True if molecule passes Lipinski Rule of Five."""
    if mol is None:
        return False
    mw   = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd  = rdMolDescriptors.CalcNumHBD(mol)
    hba  = rdMolDescriptors.CalcNumHBA(mol)
    return mw <= 500 and logp <= 5 and hbd <= 5 and hba <= 10

def veber_pass(mol) -> bool:
    """Return True if molecule passes Veber oral bioavailability criteria."""
    if mol is None:
        return False
    rotbonds = rdMolDescriptors.CalcNumRotatableBonds(mol)
    tpsa     = rdMolDescriptors.CalcTPSA(mol)
    return rotbonds <= 10 and tpsa <= 140

# PAINS filter
_pains_params = FilterCatalogParams()
_pains_params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
_pains_catalog = FilterCatalog(_pains_params)

def pains_free(mol) -> bool:
    """Return True if molecule does NOT match any PAINS filter."""
    if mol is None:
        return False
    return not _pains_catalog.HasMatch(mol)

def filter_library(smiles_list: list) -> tuple:
    """
    Apply a sequential filter cascade to a compound library.

    Returns:
        kept_smiles: list of SMILES passing all filters
        stats: dict of filter pass counts
    """
    stats = {'total': len(smiles_list), 'valid_smiles': 0,
             'lipinski': 0, 'veber': 0, 'pains': 0}
    kept = []

    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        stats['valid_smiles'] += 1

        if not lipinski_pass(mol):
            continue
        stats['lipinski'] += 1

        if not veber_pass(mol):
            continue
        stats['veber'] += 1

        if not pains_free(mol):
            continue
        stats['pains'] += 1

        kept.append(smi)

    return kept, stats

# -------------------------------------------------------
# 2. QSAR scoring model training
# -------------------------------------------------------
def train_qsar_scorer(active_smiles: list, decoy_smiles: list) -> XGBClassifier:
    """
    Train a QSAR binary classifier on known actives and property-matched decoys.
    Returns a fitted XGBClassifier.
    """
    all_smiles = active_smiles + decoy_smiles
    labels = np.array([1] * len(active_smiles) + [0] * len(decoy_smiles))

    # ECFP4 featurization
    X = np.zeros((len(all_smiles), 2048), dtype=np.uint8)
    valid_mask = np.ones(len(all_smiles), dtype=bool)
    for i, smi in enumerate(all_smiles):
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            valid_mask[i] = False
            continue
        fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
        fp.ToNumpyArray(X[i])

    X = X[valid_mask]
    y = labels[valid_mask]

    pos_weight = (y == 0).sum() / max((y == 1).sum(), 1)
    clf = XGBClassifier(
        n_estimators=300,
        max_depth=6,
        learning_rate=0.05,
        scale_pos_weight=pos_weight,
        subsample=0.8,
        use_label_encoder=False,
        eval_metric='logloss',
        random_state=42
    )
    clf.fit(X, y)
    logger.info(f"Trained QSAR model on {y.sum()} actives and {(y==0).sum()} decoys")
    return clf

def score_library(clf: XGBClassifier, smiles_list: list) -> np.ndarray:
    """
    Score a filtered compound library with a trained QSAR model.
    Returns array of active probabilities.
    """
    X = np.zeros((len(smiles_list), 2048), 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=2048)
            fp.ToNumpyArray(X[i])
    return clf.predict_proba(X)[:, 1]

# -------------------------------------------------------
# 3. Full screening pipeline
# -------------------------------------------------------
def run_screening_pipeline(
    library_smiles: list,
    active_smiles: list,
    decoy_smiles: list,
    top_k: int = 500
) -> pd.DataFrame:
    """
    Complete screening pipeline: filter, score, rank, return top compounds.
    """
    # Step 1: Filter
    logger.info(f"Filtering {len(library_smiles)} compounds...")
    filtered_smiles, stats = filter_library(library_smiles)
    logger.info(f"Filter stats: {stats}")
    logger.info(f"Remaining after filters: {len(filtered_smiles)}")

    # Step 2: Train QSAR model
    logger.info("Training QSAR model...")
    clf = train_qsar_scorer(active_smiles, decoy_smiles)

    # Step 3: Score library
    logger.info("Scoring filtered library...")
    scores = score_library(clf, filtered_smiles)

    # Step 4: Rank and select top-k
    results = pd.DataFrame({
        'smiles': filtered_smiles,
        'qsar_score': scores
    }).sort_values('qsar_score', ascending=False).head(top_k)

    return results

# Example usage
# library = pd.read_csv('zinc_drug_like.csv')['smiles'].tolist()
# actives = pd.read_csv('target_actives.csv')['smiles'].tolist()
# decoys  = pd.read_csv('dud_e_decoys.csv')['smiles'].tolist()
# hits = run_screening_pipeline(library, actives, decoys, top_k=1000)
# hits.to_csv('virtual_screening_hits.csv', index=False)

Evaluation: Enrichment Factor, ROC-AUC, and BEDROC

Virtual screening is evaluated using metrics that measure how efficiently actives are concentrated in the top of the ranked list:

import numpy as np
from sklearn.metrics import roc_auc_score

def enrichment_factor(y_true: np.ndarray, y_score: np.ndarray, fraction: float = 0.01) -> float:
    """
    Compute enrichment factor at given fraction of ranked library.

    EF = (actives in top fraction) / (total actives * fraction)
    EF = 1 is random; maximum EF = 1/fraction if all top hits are active.
    """
    n_total   = len(y_true)
    n_top     = max(1, int(n_total * fraction))
    n_actives = y_true.sum()

    if n_actives == 0:
        return 0.0

    sorted_idx = np.argsort(y_score)[::-1]  # Descending
    top_idx    = sorted_idx[:n_top]
    actives_in_top = y_true[top_idx].sum()

    baseline = n_actives * fraction
    return actives_in_top / baseline

# Example: 100 actives in 10,000 compound library (1% hit rate)
# After screening with ML model, 40 out of 100 top-ranked compounds are active
# EF@1% = (40 / 100) / 0.01 = 40  (40-fold enrichment over random)

# Using rdkit's implementation of BEDROC
from rdkit.ML.Scoring.Scoring import CalcBEDROC, CalcEnrichment, CalcAUC

def evaluate_screening(y_true: np.ndarray, y_score: np.ndarray) -> dict:
    """Comprehensive evaluation for virtual screening."""
    sorted_idx = np.argsort(y_score)[::-1]
    y_sorted = y_true[sorted_idx].reshape(-1, 1).tolist()  # rdkit format

    return {
        'auroc': roc_auc_score(y_true, y_score),
        'ef_1pct': enrichment_factor(y_true, y_score, 0.01),
        'ef_5pct': enrichment_factor(y_true, y_score, 0.05),
        'bedroc_20': CalcBEDROC(y_sorted, 0, 20.0),
    }

Practical Tips: Decoy Selection and Screening Bias

Decoy Selection

The quality of your decoy set is as important as your active set. Decoys should be:

• Property-matched: Similar molecular weight, LogP, HBD, HBA, and ring count to the actives. Otherwise, the model learns to discriminate on physicochemical properties rather than target-relevant features.
• Chemically diverse from actives: Low Tanimoto similarity to actives (below 0.4) to avoid introducing near-actives as negatives.
• Not assumed inactive: Randomly selected compounds from commercial databases are not confirmed inactives. True enrichment rates depend on the actual hit rate in the library.

The DUD-E (Directory of Useful Decoys Enhanced) provides standardized property-matched decoys for 102 drug targets. DEKOIS 2.0 provides an alternative with stricter similarity controls.

Bias in Docking Benchmarks

Reported enrichment factors for docking methods are often inflated due to:

• Analogue bias: Known actives used for benchmarking were co-crystallized with the target, making their binding poses easy for docking to reproduce. They are not representative of prospective screening scenarios.
• Decoy bias: Property-matching does not fully eliminate pharmacophore differences that make decoys trivially separable from actives.
• Retrospective vs. prospective: Retrospective benchmarks always overestimate prospective performance. When possible, evaluate on prospective screening data where the test compounds were not available during model development.

Further Reading

• RDKit Scoring Metrics (BEDROC, EF)
• DUD-E Docking Benchmark Dataset
• AutoDock Vina
• GNINA: ML-Based Docking
• Truchon and Bayly (2007): BEDROC Metric
• ZINC Virtual Compound Library

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.