Cheminformatics Pipelines for Drug Discovery

Before any molecule enters an ML model, it must be standardised. Inconsistent representations are the leading cause of irreproducible QSAR models. A compound submitted as a sodium salt will have different descriptors than its neutral parent — even though they are the same therapeutic entity. Tautomers can inflate apparent dataset diversity. Duplicate structures inflate the training set and poison generalisation metrics.

This guide builds a complete standardisation pipeline using RDKit and the MolVS/ChEMBL standardiser conventions, with a Dask backend for processing >10M compounds.

Molecule Standardisation with RDKit

from rdkit import Chem
from rdkit.Chem.MolStandardize import rdMolStandardize
from rdkit.Chem import Descriptors
from typing import Optional
import logging

logger = logging.getLogger(__name__)

# Initialise standardisation components once (thread-safe singletons)
_LARGEST_FRAGMENT_CHOOSER = rdMolStandardize.LargestFragmentChooser()
_UNCHARGER = rdMolStandardize.Uncharger()
_NORMALIZER = rdMolStandardize.Normalizer()
_TAUTOMER_ENUMERATOR = rdMolStandardize.TautomerEnumerator()
_TAUTOMER_ENUMERATOR.SetMaxTautomers(64)
_TAUTOMER_ENUMERATOR.SetRemoveBondStereo(False)  # Preserve E/Z on tautomer selection


def standardise_smiles(smiles: str, neutralise: bool = True) -> Optional[str]:
    """
    Full standardisation pipeline for a SMILES string.

    Steps:
        1. Parse SMILES
        2. Remove stereo-chemistry on parses-fail atoms
        3. Normalise (nitro groups, etc.)
        4. Remove salts / choose largest fragment
        5. Neutralise charges (optional)
        6. Canonicalise tautomer
        7. Return canonical SMILES

    Args:
        smiles: Input SMILES (may be non-canonical, charged salt, etc.)
        neutralise: Whether to neutralise charges after desalting.

    Returns:
        Canonical standardised SMILES, or None if parsing fails.
    """
    if not smiles or not isinstance(smiles, str):
        return None

    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        logger.debug("SMILES parse failed: %s", smiles[:80])
        return None

    try:
        mol = _NORMALIZER.normalize(mol)
        mol = _LARGEST_FRAGMENT_CHOOSER.choose(mol)

        if neutralise:
            mol = _UNCHARGER.uncharge(mol)

        mol = _TAUTOMER_ENUMERATOR.Canonicalize(mol)

        canonical = Chem.MolToSmiles(mol, isomericSmiles=True, canonical=True)
        return canonical

    except Exception as exc:
        logger.warning("Standardisation error for %s: %s", smiles[:80], exc)
        return None


def smiles_to_inchikey(smiles: str) -> Optional[str]:
    """Convert SMILES to InChIKey for deduplication."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    from rdkit.Chem.inchi import MolToInchiKey, MolToInchi
    inchi = MolToInchi(mol)
    if inchi is None:
        return None
    return MolToInchiKey(inchi)

Pipeline with Error Handling and Statistics

import pandas as pd
from dataclasses import dataclass, field


@dataclass
class StandardisationResult:
    """Summary statistics from a standardisation run."""
    n_input: int = 0
    n_parse_failed: int = 0
    n_standardised: int = 0
    n_duplicates_removed: int = 0
    n_output: int = 0
    failure_reasons: dict = field(default_factory=dict)

    @property
    def success_rate(self) -> float:
        return self.n_standardised / self.n_input if self.n_input > 0 else 0.0


def standardise_dataframe(
    df: pd.DataFrame,
    smiles_col: str = "smiles",
    id_col: Optional[str] = None,
    neutralise: bool = True,
    deduplicate: bool = True,
    dedup_on: str = "inchikey",
) -> tuple[pd.DataFrame, StandardisationResult]:
    """
    Standardise all SMILES in a DataFrame and deduplicate.

    Args:
        df: Input DataFrame with a SMILES column.
        smiles_col: Column name containing SMILES.
        id_col: Optional compound ID column. Retained in output.
        neutralise: Neutralise charges during standardisation.
        deduplicate: Remove duplicate structures.
        dedup_on: 'inchikey' (recommended) or 'canonical_smiles'.

    Returns:
        Tuple of (cleaned DataFrame, StandardisationResult).
    """
    result = StandardisationResult(n_input=len(df))
    df = df.copy()

    # Step 1: Standardise SMILES
    df["canonical_smiles"] = df[smiles_col].apply(
        lambda s: standardise_smiles(s, neutralise=neutralise)
    )
    failed_mask = df["canonical_smiles"].isna()
    result.n_parse_failed = int(failed_mask.sum())
    df = df[~failed_mask].reset_index(drop=True)
    result.n_standardised = len(df)

    # Step 2: Compute InChIKeys for deduplication
    if deduplicate and dedup_on == "inchikey":
        df["inchikey"] = df["canonical_smiles"].apply(smiles_to_inchikey)
        n_before = len(df)
        df = df.drop_duplicates(subset="inchikey", keep="first").reset_index(drop=True)
        result.n_duplicates_removed = n_before - len(df)
    elif deduplicate:
        n_before = len(df)
        df = df.drop_duplicates(subset="canonical_smiles", keep="first").reset_index(drop=True)
        result.n_duplicates_removed = n_before - len(df)

    result.n_output = len(df)
    return df, result

Large-Scale Processing with Dask

import dask.dataframe as dd
from dask import delayed
import numpy as np

def standardise_partition(partition: pd.DataFrame, smiles_col: str, neutralise: bool) -> pd.DataFrame:
    """Applied to each Dask partition."""
    partition = partition.copy()
    partition["canonical_smiles"] = partition[smiles_col].apply(
        lambda s: standardise_smiles(s, neutralise=neutralise)
    )
    return partition


def process_large_library(
    input_parquet: str,
    output_parquet: str,
    smiles_col: str = "smiles",
    npartitions: int = 64,
    neutralise: bool = True,
) -> dict:
    """
    Standardise a large compound library stored as Parquet using Dask.

    Suitable for 10M+ compound libraries. Uses local processes scheduler
    by default; switch to dask.distributed for cluster execution.

    Args:
        input_parquet: Path to input Parquet file (or glob for multiple files).
        output_parquet: Output directory for partitioned Parquet.
        smiles_col: Column containing SMILES strings.
        npartitions: Dask partitioning (tune to ~100k rows/partition).
        neutralise: Whether to neutralise charges.

    Returns:
        Summary dict with processing statistics.
    """
    ddf = dd.read_parquet(input_parquet)
    n_input = len(ddf)  # Triggers metadata read only

    # Apply standardisation per partition
    meta = ddf._meta.copy()
    meta["canonical_smiles"] = pd.Series(dtype="string")

    ddf_std = ddf.map_partitions(
        standardise_partition,
        smiles_col=smiles_col,
        neutralise=neutralise,
        meta=meta,
    )

    # Remove failed parses
    ddf_clean = ddf_std[ddf_std["canonical_smiles"].notna()]

    # Write to Parquet
    ddf_clean.to_parquet(
        output_parquet,
        engine="pyarrow",
        compression="snappy",
        write_index=False,
    )

    return {
        "status": "complete",
        "output_path": output_parquet,
        "npartitions": npartitions,
        "note": "Run deduplication step on output using pandas groupby inchikey for final dedup.",
    }

Tautomer and Stereoisomer Enumeration

from rdkit.Chem.EnumerateStereoisomers import (
    EnumerateStereoisomers,
    StereoEnumerationOptions,
)


def enumerate_tautomers(smiles: str, max_tautomers: int = 32) -> list[str]:
    """
    Enumerate tautomers for a SMILES string.
    Returns canonical SMILES for each unique tautomer.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return []

    enumerator = rdMolStandardize.TautomerEnumerator()
    enumerator.SetMaxTautomers(max_tautomers)

    tautomers = enumerator.Enumerate(mol)
    smiles_set = {Chem.MolToSmiles(t) for t in tautomers if t is not None}
    return sorted(smiles_set)


def enumerate_stereoisomers(
    smiles: str,
    max_isomers: int = 16,
    try_embedding: bool = False,
) -> list[str]:
    """
    Enumerate undefined stereocenters in a SMILES.

    Args:
        smiles: Input SMILES with undefined stereocenters (marked as @? or missing)
        max_isomers: Maximum number of stereoisomers to generate
        try_embedding: If True, remove stereoisomers that fail 3D embedding

    Returns:
        List of SMILES for each stereoisomer.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return []

    opts = StereoEnumerationOptions(
        tryEmbedding=try_embedding,
        unique=True,
        maxIsomers=max_isomers,
        rand=0xDEADBEEF,  # Deterministic ordering
    )

    isomers = list(EnumerateStereoisomers(mol, options=opts))
    return [Chem.MolToSmiles(iso) for iso in isomers]

Common Standardisation Pitfalls

Key Takeaways

• Standardisation is the most impactful data preprocessing step in drug discovery ML. A rigorous pipeline: normalise → desalt → neutralise → canonicalise tautomer → canonical SMILES → InChIKey dedup.
• Always use InChIKey for deduplication — string-based SMILES deduplication misses equivalent structures with different atom orderings and tautomeric forms.
• For multi-million compound libraries, Dask map_partitions with ~100k rows per partition provides good parallelism. Deduplicate after writing Parquet to avoid cross-partition groupby overhead.
• For ADMET property prediction (logP, pKa, solubility), preserve the original ionisation state. For fingerprint-based ML models, neutralise first to avoid spurious descriptor differences.

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.