Bio Qsar Modeling

Version Compatibility

Reference examples target: chemprop 2.2.x (major API change from 1.x), RDKit 2024.09+, scikit-learn >=1.4,<1.6, MAPIE >=0.8,<1.0 for the MapieRegressor example, shap 0.44+, and pytorch 2.1+. Recheck examples before widening these bounds because Chemprop, scikit-learn calibration, and MAPIE interfaces evolve independently.

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, and split strategy determines whether the model captures transferable chemical signal or memorizes the training data. chemprop D-MPNN with optional Morgan / RDKit descriptors is a useful open-source approach; transformer-based methods (MolFormer, Uni-Mol, ChemBERTa) should be compared on the same split and endpoint. The OECD validation principles support transparent documentation and evaluation of (Q)SAR models, but following them does not by itself confer regulatory acceptance.

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

Decision: Compare a fingerprint-based baseline with chemprop under the same split and endpoint. Add a pretrained transformer or 3D model only when its representation, compute cost, and validation design fit the deployment question; dataset size alone does not determine the winner.

Decision Tree by Scenario

OECD 5 Principles

The OECD principles were agreed in 2004; the 2007 guidance explains their application:

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 measures of goodness-of-fit, robustness, and predictivity: external test set and suitable validation
5. Mechanistic interpretation, if possible: biological/chemical rationale where available

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

Applicability Domain Methods

Ensemble disagreement is one useful uncertainty diagnostic, not a formally defined applicability domain or calibrated coverage guarantee. If using a threshold such as a training-distribution percentile, label it as a project-defined heuristic and validate it prospectively.

chemprop 2.x Training (CLI)

Goal: Train five replicated chemprop runs, each containing a five-model D-MPNN ensemble with RDKit 2D descriptor features and a scaffold-balanced train/validation/test split.

Approach: For current chemprop 2.x, invoke chemprop train with --molecule-featurizers rdkit_2d, --num-replicates 5, --ensemble-size 5, and --split scaffold_balanced. Replicates repeat splitting/training with incremented seeds; they are not five-fold cross-validation. Confirm the exact flags with chemprop train --help because the v2 CLI continues to evolve.

# 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 \
    --num-replicates 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: include current v2 RDKit descriptors, which are scaled by default (the legacy v1-normalized generator is v1_rdkit_2d_normalized)
• --num-replicates 5: repeat the split/training workflow with successive seeds; this replaced --num-folds in chemprop 2.1
• --ensemble-size 5: train five models per replicate for an ensemble prediction
• --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 replicates x 5 ensemble members). Report which predictions are being aggregated and treat ensemble standard deviation as an uncertainty diagnostic, not a calibrated guarantee.

At prediction time, uncertainty output is opt-in and requires the actual saved model paths:

chemprop predict --test-path test.csv \
    --model-paths path/to/model_1.ckpt path/to/model_2.ckpt \
    --uncertainty-method ensemble \
    --preds-path predictions.csv

Scaffold-Balanced Split

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 assign whole scaffold groups to train, validation, or test. chemprop's --split scaffold_balanced implements a scaffold-based allocation; --class-balance is a separate training option and does not make this split outcome-stratified. The chemprop default split is random, so request scaffold-balanced explicitly when it matches the deployment question.

scaffold_balanced assigns each scaffold group to one of train / validation / test, reducing direct scaffold leakage. It is not universally the correct validation design: time splits, externally defined series, grouped cross-validation, and prospective tests may better represent a particular deployment setting.

Choose and document the primary split before model selection. A random split can answer an interpolation question but often shares close analogues across partitions; a scaffold split tests transfer across scaffold groups; a time or prospective split tests the historical deployment process. If several splits are reported, interpret their differences as split-specific sensitivity rather than a universal "true generalization gap."

Conformal Prediction for Calibrated Uncertainty

Use conformal prediction when calibrated marginal coverage under the stated exchangeability assumptions matters. Ensemble variance is simpler, but it is not a substitute for a conformal guarantee.

# 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 0.05 targets 95% marginal coverage and alpha 0.10 targets 90%, subject to the conformal method's assumptions. MAPIE supports the sklearn baseline directly; integrating chemprop requires a separately implemented and tested compatible wrapper.

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). A PyTorch attribution method must be adapted to the model's graph inputs and validated; the chemprop 2.x CLI does not provide the --uncertainty-method classification atom-attribution interface. Use directly supported fingerprint SHAP for the classical baseline unless a tested graph-attribution implementation is available.

Bayesian Optimization for Active Learning

import numpy as np
from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

gp = GaussianProcessRegressor(kernel=RBF(length_scale=1.0), random_state=42)
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):
    improvement = mu - y_best - xi
    ei = np.zeros_like(mu, dtype=float)
    nonzero = sigma > 0
    z = improvement[nonzero] / sigma[nonzero]
    ei[nonzero] = (
        improvement[nonzero] * norm.cdf(z)
        + sigma[nonzero] * norm.pdf(z)
    )
    return ei

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 probabilities are not guaranteed to be calibrated. Use Platt (logistic) or isotonic calibration for binary probabilities, choosing the method with a held-out calibration set. --metric roc evaluates ranking and does not automatically calibrate chemprop probabilities; export validation probabilities and fit the calibrator 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

Multitask learning can help when endpoints share predictive signal or data, but negative transfer is also possible. Compare single-task and multitask models under identical splits rather than assuming improvement from endpoint relatedness.

Per-Tool Failure Modes