Qsar modeling

Version Compatibility

Reference examples tested with: chemprop 2.0+ (major API change from 1.x), RDKit 2024.09+, scikit-learn 1.4+, MAPIE 0.8+ (conformal prediction), shap 0.44+, pytorch 2.1+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: chemprop train --help (chemprop 2.x); chemprop_train --help (1.x legacy)

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

QSAR Modeling

Build quantitative structure-activity relationship models from molecular structure inputs. The choice of model + featurization + split strategy determines whether the model captures real chemical signal or memorizes the training data. chemprop D-MPNN with optional Morgan / RDKit descriptors is the modern open-source standard; transformer-based methods (MolFormer, Uni-Mol, ChemBERTa) compete on benchmarks but offer minimal practical gain when domain-specific data is sparse. The OECD 5 principles structure the model for regulatory acceptance: defined endpoint, unambiguous algorithm, defined applicability domain, validation, and mechanistic interpretation.

For descriptor/fingerprint choices, see chemoinformatics/molecular-descriptors. For ADMET-specific QSAR, see chemoinformatics/admet-prediction. For molecular standardization (critical upstream), see chemoinformatics/molecular-standardization.

Model Taxonomy

Model	Architecture	Use case	Fails when
Random Forest + ECFP4	Classical baseline	Small data (<200 compounds), interpretability	Saturates at ~AUC 0.85 for complex endpoints
chemprop D-MPNN	Directed message passing	Modern default; 100-10k compounds	Very small datasets (<100) overfits
chemprop D-MPNN + RDKit 2D	Hybrid graph + descriptors	hERG SOTA (Cai 2024 AUC 0.956)	Diminishing returns at large data
MolFormer	Transformer (87M params)	Large public training data benefit	Compute overhead; OOD risk
Uni-Mol	3D-aware transformer	3D-relevant endpoints (binding)	Requires 3D conformers
ChemBERTa-2	Transformer (77M params)	SMILES language model	Large training data needed
Gaussian Process + ECFP4	Probabilistic	Active learning; uncertainty	O(N^3) scaling
MultiTask DNN	Joint training	Multiple endpoints	Data must overlap
AttentiveFP	GNN with attention	Pre-chemprop SOTA	Less actively maintained
KPGT / Graphormer / GROVER	Pretrained GNN	Pretraining lift	Diminishing returns

Decision: For 200-10k compounds, chemprop 2.0 D-MPNN + RDKit 2D descriptors is the modern standard. For <200 compounds, Random Forest + ECFP4 is competitive and more interpretable.

Decision Tree by Scenario

Dataset size	Endpoint type	Model
<50	Anything	Don't model; use as test set for literature models
50-200	Regression / classification	RF + ECFP4 + cross-validation
200-1k	Regression	chemprop D-MPNN + RDKit 2D, 5-fold
1k-10k	Anything	chemprop D-MPNN + RDKit 2D ensemble
10k-100k	Classification	chemprop ensemble OR MolFormer fine-tuning
>100k	Anything	MolFormer / Uni-Mol pretrained + LoRA fine-tune
Multi-task	Related endpoints (CYP3A4, CYP2D6, etc.)	chemprop MultiTask
3D-relevant	Binding, conformer-dependent	Uni-Mol with conformer ensemble

OECD 5 Principles

For regulatory use (REACH, ECHA, FDA submissions):

1. Defined endpoint: specific bioassay, units, threshold definitions
2. Unambiguous algorithm: reproducible code, fixed random seeds, version-pinned dependencies
3. Defined applicability domain (AD): where the model is valid
4. Appropriate statistical validation: external test set, cross-validation
5. Mechanistic interpretation: biological/chemical rationale where possible

For non-regulatory QSAR, all 5 still good practice; especially AD definition is critical.

Applicability Domain Methods

Method	Definition	Pro	Con
Ensemble variance (RECOMMENDED for chemprop)	Std across N-model ensemble predictions	Built-in to chemprop predict; no extra code	Assumes ensemble diversity (random seeds give independent fits)
kNN distance	Mean Tanimoto to k nearest in training	Easy to interpret	Doesn't account for label distribution
Leverage	Hat matrix diagonal	Statistical	Linear assumptions
KDE on PCA	Density in feature space	Captures multivariate structure	Density choice subjective
Mahalanobis distance	Covariance-aware distance	Theoretically motivated	High-dim instability
Conformal prediction	Per-prediction confidence interval	Distribution-free guarantee	Requires calibration set + sklearn-compatible estimator
Bayesian / MC-dropout	Posterior or dropout variance	Direct uncertainty	Computational cost
Tanimoto coverage	At least 1 NN within threshold	Practical	Threshold subjective

Operational rule: For chemprop, set --ensemble-size 5 at training; at predict time, flag predictions with ensemble std > P95 of the training-set ensemble std distribution as out-of-AD. This is the standard chemprop-native AD measure and avoids the overhead of wrapping with MAPIE/conformal.

chemprop 2.0 Training (CLI)

Goal: Train a 5-fold x 5-model chemprop D-MPNN ensemble with RDKit 2D descriptor features and scaffold-balanced split.

Approach: Invoke chemprop train with --molecule-featurizers rdkit_2d_normalized, --num-folds 5, --ensemble-size 5, and --split scaffold_balanced to produce 25 models for ensemble mean prediction + variance-based applicability domain.

# chemprop 2.x CLI (current): use 'chemprop train' (space; dashes not underscores)
chemprop train \
    --data-path data.csv \
    --task-type classification \
    --save-dir model_dir \
    --molecule-featurizers rdkit_2d_normalized \
    --num-folds 5 \
    --ensemble-size 5 \
    --epochs 50 \
    --batch-size 128 \
    --split scaffold_balanced \
    --split-sizes 0.8 0.1 0.1 \
    --metric roc

# chemprop 1.x legacy CLI (for backwards reference):
# chemprop_train --data_path data.csv --dataset_type classification ...

Key flags (chemprop 2.x):

• --molecule-featurizers rdkit_2d_normalized: include physchem descriptors (was --features_generator in 1.x)
• --num-folds 5: 5-fold cross-validation
• --ensemble-size 5: 5 models per fold for uncertainty
• --split scaffold_balanced: prevent scaffold leakage (was --split_type in 1.x)
• --split-sizes 0.8 0.1 0.1: 80/10/10 train/val/test

Total models: 25 (5 folds x 5 ensemble). Use ensemble mean as prediction; ensemble std as uncertainty.

Scaffold-Balanced Split (chemprop 2.0 default)

Goal: Partition a SMILES dataset into train/val/test such that no Bemis-Murcko scaffold appears in more than one split (prevents chemotype leakage).

Approach: Group compounds by scaffold and greedily assign whole scaffolds to train until target fraction, then val, then test. chemprop's --split scaffold_balanced implements this with optional class-stratification.

scaffold_balanced split assigns entire scaffolds to one of train / val / test, preventing chemotype leakage. Modern QSAR uses scaffold split exclusively.

Split type	Random	Scaffold balanced
Train AUC	0.99	0.95
Test AUC	0.95	0.75-0.85
Generalization	Optimistic	Realistic

The gap between random and scaffold split is the true generalization gap; report both.

Conformal Prediction for Calibrated Uncertainty

Use conformal prediction when calibrated coverage guarantees matter (regulatory submission, decision thresholds with cost asymmetry). For routine QSAR, chemprop ensemble variance is simpler and sufficient.

# MAPIE expects a scikit-learn-compatible estimator (.fit / .predict / .predict_proba).
# chemprop 2.x is NOT scikit-learn-compatible out of the box -- either wrap chemprop
# in a thin sklearn estimator class or use MAPIE only with the sklearn baseline.
from mapie.regression import MapieRegressor
from sklearn.ensemble import RandomForestRegressor

base = RandomForestRegressor(n_estimators=500, random_state=42)
mapie = MapieRegressor(estimator=base, method='plus', cv=5)
mapie.fit(X_train, y_train)
y_pred, y_intervals = mapie.predict(X_test, alpha=0.1)  # alpha=0.1 -> 90% coverage

Alpha choice: 0.05 (95% coverage) for safety endpoints; 0.10 (90%) for activity ranking. For chemprop, the simpler ensemble-variance route (set --ensemble-size 5 and use prediction std as AD signal) is preferred unless conformal coverage guarantees are specifically required.

SHAP / Atomic Attribution

For mechanistic interpretation:

For a scikit-learn-style model (e.g., Random Forest baseline on ECFP4), SHAP integrates directly:

import shap
from sklearn.ensemble import RandomForestClassifier

# X_train / X_test are Morgan fingerprint arrays (n_samples, n_bits)
model = RandomForestClassifier(n_estimators=500, random_state=42).fit(X_train, y_train)
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# Per-bit contribution; for atomic interpretation, map bits back to
# generating atoms via AllChem.GetMorganFingerprintAsBitVect(mol, ..., bitInfo=bi)
# and aggregate SHAP across all bits triggered by each atom.

For chemprop D-MPNN, SHAP requires a custom wrapper (chemprop is not sklearn-compatible). Use shap.GradientExplainer on the underlying PyTorch module after exposing it via a wrapper class. In practice, chemprop's built-in atom attribution (chemprop predict --uncertainty-method classification with atom-level outputs) is the supported route for atomic interpretation. Use SHAP only when comparing chemprop to a classical baseline.

Bayesian Optimization for Active Learning

from sklearn.gaussian_process import GaussianProcessRegressor

gp = GaussianProcessRegressor(kernel='rbf')
gp.fit(X_train, y_train)
mu, sigma = gp.predict(X_pool, return_std=True)

# Expected Improvement
def expected_improvement(mu, sigma, y_best, xi=0.01):
    z = (mu - y_best - xi) / sigma
    return (mu - y_best - xi) * norm.cdf(z) + sigma * norm.pdf(z)

ei = expected_improvement(mu, sigma, y_train.max())
next_to_test = X_pool[ei.argmax()]

For chemprop + active learning, replace GP with chemprop ensemble + ensemble variance.

Calibration (Platt / Isotonic)

Deep learning native probabilities rarely well-calibrated. Apply Platt scaling on validation set:

from sklearn.calibration import CalibratedClassifierCV
from sklearn.ensemble import RandomForestClassifier

# For sklearn-compatible models, calibration is direct:
base = RandomForestClassifier(n_estimators=500, random_state=42).fit(X_train, y_train)
calibrated = CalibratedClassifierCV(base, method='isotonic', cv='prefit')
calibrated.fit(X_val, y_val)
y_proba_calibrated = calibrated.predict_proba(X_test)

Use Platt (logistic) for binary; isotonic for general. For chemprop, native calibration is applied automatically when training with classification + --metric roc; for further re-calibration, save validation probabilities and fit isotonic regression externally:

from sklearn.isotonic import IsotonicRegression
iso = IsotonicRegression(out_of_bounds='clip').fit(val_chemprop_probs, val_true)
test_calibrated = iso.predict(test_chemprop_probs)

Multi-Task QSAR

Train multiple related endpoints jointly:

df = pd.DataFrame({
    'smiles': [...],
    'CYP1A2_inhibition': [...],
    'CYP2D6_inhibition': [...],
    'CYP3A4_inhibition': [...],
})
df.to_csv('multitask.csv', index=False)

chemprop train --data-path multitask.csv --task-type classification \
               --target-columns CYP1A2_inhibition CYP2D6_inhibition CYP3A4_inhibition \
               --save-dir multitask_model

Joint training improves performance when endpoints are correlated (Wu 2018 MoleculeNet). For independent endpoints, multitask hurts.

Per-Tool Failure Modes

Random split for QSAR

Trigger: Default sklearn train_test_split.

Mechanism: Compounds from same scaffold scatter across train/test; performance optimistic.

Symptom: Cross-validation AUC > 0.95; real-world performance much lower.

Fix: Use --split scaffold_balanced in chemprop 2.x (or --split_type scaffold_balanced in chemprop 1.x legacy); or scaffold_split from chemoinformatics/scaffold-analysis.

Class imbalance not handled

Trigger: 10:1 negative:positive ratio in dataset.

Mechanism: Default loss treats classes equally; model learns majority class.

Symptom: High accuracy but precision/recall on minority class poor.

Fix: Class-weighted loss; SMOTE; or report AUC/F1 not accuracy.

Over-engineered features

Trigger: Including hundreds of descriptors (e.g., rdkit_2d not normalized).

Mechanism: Some descriptors dominate scaling; model overfits.

Symptom: Validation performance differs widely across runs; high feature importance noise.

Fix: Use rdkit_2d_normalized; or restrict to <50 standard descriptors.

Missing AD assessment

Trigger: Predicting on novel chemotypes without AD check.

Mechanism: Model extrapolates; predictions unreliable.

Symptom: Confident predictions but actual values different.

Fix: Always compute kNN Tanimoto distance + ensemble variance; flag low-confidence predictions.

chemprop 1.x vs 2.x confusion

Trigger: Code/tutorial from before late 2024.

Mechanism: Major API change: chemprop_train → chemprop train; Python API redesigned.

Symptom: ImportError or different keyword arguments.

Fix: Use chemprop --version; check 2.x documentation; migrate APIs.

Pretrained Transformer overhead without data benefit

Trigger: Using MolFormer on <500 compounds.

Mechanism: Pretrained representation needs sufficient fine-tuning data to specialize.

Symptom: No improvement over chemprop; slower training.

Fix: For <500 compounds, stick with chemprop or RF. Use pretrained Transformers for >10k.

Validation leakage via standardization

Trigger: Different standardization for train vs test.

Mechanism: Train compounds with R-isomer; test compounds with S-isomer not in train.

Symptom: Test performance better than realistic.

Fix: Standardize unified preprocessing; same canonicalization train + test.

Reconciliation: Classical RF vs chemprop vs Transformer

Aspect	RF + ECFP4	chemprop D-MPNN	MolFormer
Data size sweet spot	50-1k	200-10k	>10k
Interpretability	High (Gini importance)	Medium (SHAP)	Low
Uncertainty	Bootstrap	Ensemble	MC-dropout
Hardware	CPU	GPU recommended	GPU mandatory
OOD performance	Most degraded	Better with ensemble	Worst (long-range memorization)
Production deployment	Easy (pickle)	OK (ONNX export)	Heavy

Common Errors

Symptom	Cause	Fix
chemprop hangs at start	GPU OOM	Reduce batch_size; check CUDA
All predictions same value	Constant target	Standardize labels
AUC mismatched across folds	Random seed not set	--seed 42
Test AUC = train AUC	No held-out data	Use scaffold_balanced split
Ensemble variance always small	All models identical	Set --seed per fold; check randomness
SHAP fails on D-MPNN	Custom architecture	Use GradientExplainer not TreeExplainer
MolFormer fine-tune slow	All parameters trained	Use LoRA or freeze early layers
Calibration broken	Pre-calibrated already	Skip calibration step

References

• Yang et al., J. Chem. Inf. Model. 59:3370 (2019) -- chemprop D-MPNN.
• Heid et al., J. Chem. Inf. Model. 64:9 (2024) -- chemprop 2.0 redesign.
• Wu et al., Chem. Sci. 9:513 (2018) -- MoleculeNet benchmark.
• Ross et al. (2022) -- MoLFormer architecture.
• Zhou et al., Nat. Commun. 14:3849 (2023) -- Uni-Mol 3D.
• OECD, "OECD Principles for the Validation of QSAR Models" (2007).
• OECD, "(Q)SAR Assessment Framework" (2023).
• Cortés-Ciriano & Bender, J. Chem. Inf. Model. 60:1184 (2020) -- conformal prediction for QSAR applicability domain.
• Alperstein et al., J. Cheminformatics 15:9 (2023) -- MolBERT chemistry-aware tokenization vs ChemBERTa-2.

Related Skills

• chemoinformatics/molecular-descriptors - Featurization choices
• chemoinformatics/molecular-standardization - Mandatory upstream
• chemoinformatics/scaffold-analysis - Bemis-Murcko split implementation
• chemoinformatics/admet-prediction - ADMET-specific QSAR
• chemoinformatics/generative-design - QSAR as scoring component
• machine-learning/model-validation - General ML validation principles
• machine-learning/biomarker-discovery - Adjacent ML approaches