GPTomics/bio-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