bio-qsar-modeling

What this skill does

## 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.

```bash
# 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 vari