Metadata-Version: 2.4
Name: cyp-substrate-predictor
Version: 0.2.0
Summary: Pharmacophore-aware fingerprint classifier for CYP substrate prediction (TDC ADMET)
License: Apache-2.0
Project-URL: Homepage, https://github.com/seqeralabs/cyp-substrate-predictor
Project-URL: Source, https://github.com/seqeralabs/cyp-substrate-predictor
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: rdkit
Requires-Dist: scikit-learn
Requires-Dist: numpy
Requires-Dist: pandas
Provides-Extra: eval
Requires-Dist: PyTDC; extra == "eval"
Requires-Dist: setuptools<81; extra == "eval"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# cyp-substrate-predictor

A small, CPU-only, interpretable classifier for predicting whether a small
molecule is a **substrate of a cytochrome-P450 enzyme** from its SMILES string.
It pairs two complementary molecular representations — **count Morgan
fingerprints** and the **Extended Reduced Graph (ErG)** pharmacophore descriptor —
with a **random-subspace bagged logistic regression**.

On the Therapeutics Data Commons (TDC) ADMET benchmark
`CYP2D6_Substrate_CarbonMangels` it reaches **AUPRC 0.780 ± 0.007** under the
official 5-seed protocol, with ~0.14M parameters and no GPU.

---

## Why it works (the science)

**The recognition problem.** CYP2D6 substrate space is
governed by a well-characterised pharmacophore: substrates are typically
**lipophilic bases** carrying a **protonatable nitrogen** that forms a salt
bridge with an acidic residue in the enzyme's active site, with the site of
oxidation located roughly **5–7 Å** from that basic centre. Good discrimination
therefore depends on capturing both *local substructure* (what chemical groups
are present) and *pharmacophore geometry* (which functional centres sit at which
topological distances). No single fingerprint captures both well.

**Two complementary feature blocks.**

- **Count Morgan (radius 1, 512 bits).** Morgan/ECFP fingerprints encode circular
  atom environments; using **counts** (not bits) preserves how *often* each local
  environment occurs (ring counts, repeated substitution patterns), which carries
  real signal for lipophilicity and aromatic content. A small radius keeps the
  representation compact and collision-tolerant on a few-hundred-molecule dataset.
- **ErG (Extended Reduced Graph, 315 float features).** ErG reduces a molecule to
  pharmacophore-typed nodes (H-bond donor/acceptor, positive/negative ionisable,
  aromatic, hydrophobic) and encodes **pairs of node types at each topological
  distance**. This is almost exactly the representation the CYP substrate
  pharmacophore calls for — e.g. *"positively ionisable centre at distance d from
  an aromatic/hydrophobic centre"* — and it captures relationships that circular
  fingerprints miss.

Concatenating the two (827 features) gives the model both a local-substructure
view and a pharmacophore-geometry view of every molecule.

**Why random-subspace bagged logistic regression.** The datasets are small
(~500 labelled molecules) and the feature space is high-dimensional, sparse, and
heterogeneously scaled (integer Morgan counts alongside continuous ErG values).
That is a regime where large non-linear models overfit and are hard to trust. The
model here is deliberately simple and heavily regularised:

- a **linear base learner** (L2-regularised logistic regression) that cannot
  overfit a few hundred points in the way a deep or boosted model can;
- wrapped in a **BaggingClassifier of 200 estimators**, each trained on a
  bootstrap of the rows *and* a **random 85% subspace of the features**
  (`bootstrap_features=True`). The random subspace method decorrelates the base
  learners across the correlated fingerprint dimensions, turning bagging into an
  effective variance-reduction and implicit feature-selection mechanism;
- with **StandardScaler inside each bag's pipeline**, so every base learner
  standardises only its own feature subset — reconciling the count/ErG scale
  mismatch without leaking global statistics.

The result is a model that behaves like a smoothly regularised linear ensemble:
strong on small data, fast on CPU, and transparent about what it uses.

---

## Results

Official TDC ADMET `admet_group` protocol (5 seeds, `evaluate_many`), metric AUPRC:

| Benchmark | This method | For reference (TDC leaderboard #1) |
|---|---|---|
| `CYP2D6_Substrate_CarbonMangels` | **0.780 ± 0.007** | ContextPred 0.736 ± 0.024 |

The released weights (a single model fitted on the full `train_val` split, rather
than the protocol's 5 refits) score **0.779** on the same test split.

---

## Pretrained weights

`cyp_substrate/weights/cyp2d6_substrate_carbonmangels.npz` (371 KB) holds a model fitted on all
532 `train_val` molecules. Scoring needs **only rdkit and numpy** — the weights are
plain arrays, not a pickled estimator:

```python
from cyp_substrate.pretrained import load

model = load()
model.predict_proba(["COc1ccc2c(c1)[C@@H]1CC3CCCC[C@]3(CCN1C)C2"])   # -> [0.836]
```

```bash
cyp-substrate-predict "CC(C)NCC(O)COc1cccc2ccccc12"    # propranolol -> 0.8028
cyp-substrate-predict --input molecules.csv --output scores.csv
python train.py                                        # regenerate the weights
pytest tests/                                          # verify them (no TDC needed)
```

The CLI reads a CSV whose header names a SMILES column (`smiles`,
`canonical_smiles`, `Drug`, `structure`) or a plain one-per-line list, carries any
`id` column through to the output, and writes
`id,smiles,p_substrate,prediction,valid`. SMILES that RDKit cannot parse are
flagged `valid=false` and left unscored rather than reported as ~0.23, the score
of an all-zero feature row.

Each of the 200 bags is a `StandardScaler → LogisticRegression`, which is affine
in the raw features, so the fitted ensemble collapses **exactly** into one dense
827-vector per bag (agreement with scikit-learn: 3.3e-16 in float64). That is what
the `.npz` stores. On the official 135-molecule test split these weights reach
**AUPRC 0.779 / AUROC 0.858**.

See [`docs/model-card.md`](docs/model-card.md) for the model card: threshold
selection, a sanity check against known CYP2D6 pharmacology, and limitations —
including benchmark label noise and why 0.78 is split-specific.

## Nextflow module

The predictor is packaged as a Nextflow module for
[registry.nextflow.io](https://registry.nextflow.io), written to the nf-core
module conventions:

```nextflow
include { CYPSUBSTRATEPREDICTOR } from 'seqera/cyp-substrate-predictor'

workflow {
    CYPSUBSTRATEPREDICTOR(
        Channel.of([[id: 'library'], file('molecules.csv')]),
        []                                  // [] = bundled weights
    )
}
```

Source lives in [`modules/seqera/cyp-substrate-predictor/`](modules/seqera/cyp-substrate-predictor/),
with its own [README](modules/seqera/cyp-substrate-predictor/README.md) covering
inputs, thresholds, and nf-core porting. Scoring is a matmul, so run **one task
per file** rather than one per molecule. Tests:

```bash
nf-test test modules/seqera/cyp-substrate-predictor/tests/main.nf.test   # needs Nextflow >= 26.04
```

Releases are automated: tag, cut a GitHub Release, and
[`.github/workflows/workflow.yml`](.github/workflows/workflow.yml) publishes to
PyPI via Trusted Publishing after verifying the built wheel. Full runbook,
including the container and module steps, is in
[`docs/releasing.md`](docs/releasing.md).

## Software

```bash
pip install cyp-substrate-predictor        # rdkit, scikit-learn, numpy, pandas
```

or from a checkout:

```bash
pip install -e .
```

Featurise + train + predict:

```python
from cyp_substrate.features import featurize
from cyp_substrate.model import build_model

X_train = featurize(train_smiles)          # count-Morgan(r1/512) ⊕ ErG(315)
clf = build_model()                        # RSM-bagged logistic, tuned defaults
clf.fit(X_train, y_train)
proba = clf.predict_proba(featurize(test_smiles))[:, 1]
```

Reproduce the TDC benchmark numbers (`pip install PyTDC`):

```bash
python evaluate.py --benchmark CYP2D6_Substrate_CarbonMangels   # -> 0.780 ± 0.007
```

Default hyperparameters (`cyp_substrate/model.py`): Morgan radius 1 / 512 bits,
`C=0.0793`, `max_features=0.85`, `n_estimators=200`, `max_samples=0.85`.

---

## License

Apache-2.0.
