bio-machine-learning-biomarker-discovery

What this skill does

## Version Compatibility

Reference examples tested with: numpy 1.26+, pandas 2.2+, scikit-learn 1.4+, boruta 0.4+, mrmr-selection 0.2+.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures

BorutaPy expects numpy arrays and breaks on newer numpy where the `np.float`/`np.int` aliases were removed -- pin a compatible numpy or use a maintained fork. On scikit-learn 1.8+ the `LogisticRegression(penalty=)` argument is deprecated (removed in 1.10) in favor of `l1_ratio`+`C`; the examples show the 1.4-1.7 form. If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

# Feature Selection for Biomarker Discovery

**"Find the biomarkers in my omics data"** -> First decide which question is being answered (all-relevant vs minimal-optimal), then select features INSIDE a resampling loop, then quantify stability -- because a selected list means little without it.
- All-relevant (which genes carry signal?): `BorutaPy(rf)`
- Minimal-optimal (smallest predictive set?): `ElasticNetCV`, `LogisticRegressionCV(penalty='elasticnet')`
- Stability (does the list reproduce?): bootstrap selection frequencies + a stability index

## The Single Most Important Modern Insight -- Most Gene Signatures Do Not Replicate, and Significance Is the Wrong Bar

A signature being "significantly associated with outcome" is near-worthless evidence: *random* gene sets -- and signatures of biologically irrelevant phenomena -- are significantly associated with breast-cancer survival, often matching published prognostic signatures, because the transcriptome is dominated by a few axes (proliferation) that almost any large gene set captures (Venet 2011). The correct null is not "no association" but *random gene sets of equal size* plus a proliferation meta-gene. Two further hard facts complete the picture: many disjoint gene lists predict equally well (Ein-Dor 2005), so non-overlap with a prior list is the *expected* result, not a contradiction; and obtaining a *stable* list (as opposed to an accurate predictor) needs on the order of thousands of samples (Ein-Dor 2006), far more than typical omics n.

The operational consequences run through every section below: report a stability index next to accuracy; benchmark against a random-signature and proliferation-meta-gene null; never interpret the specific genes a minimal-optimal selector kept as "the biomarkers"; and keep selection inside the cross-validation loop or the reported performance is fiction.

## All-Relevant vs Minimal-Optimal (the distinction usually conflated)

This axis matters more than filter/wrapper/embedded. Choosing the wrong one is the most common conceptual error in applied biomarker papers.

- **Minimal-optimal** = the *smallest* subset giving optimal prediction (LASSO, RFE, forward selection). If two genes are correlated and both informative, it keeps **one and drops the other**; the dropped gene is still biologically relevant. Minimal-optimal sets are non-unique, unstable, and systematically exclude redundant-but-real features. *Absence from a minimal-optimal set is not evidence of irrelevance.*
- **All-relevant** = *every* feature carrying information, redundant or not (Boruta: keep anything beating the best "shadow" permuted feature). This is the right framing for *biological interpretation* -- the whole co-expression module is wanted, not one representative.

Decision rule: parsimonious assay with few measurements -> minimal-optimal; understand biology / enumerate implicated genes / pathway analysis -> all-relevant; stable deployable signature -> elastic net or stability selection.

## Methods Taxonomy

| Family | Method | Optimizes | Redundancy handling | Output | Key trap |
|--------|--------|-----------|---------------------|--------|----------|
| Filter (univariate) | t-test / `SelectKBest(f_classif)` | Marginal association, one gene at a time | None (keeps correlated blocks) | Ranked list | Ignores multivariate structure; huge multiplicity |
| Filter (multivariate) | mRMR (Peng 2005) | Relevance minus redundancy | Explicit penalty | Ranked K | Greedy/first-order; K must still be chosen |
| Wrapper | RFE / RFECV; SVM-RFE | A specific model's accuracy | Indirect | Ranked subset | Expensive; **must be inside CV**; SVM-RFE needs a linear kernel |
| Embedded | LASSO (Tibshirani 1996) | Prediction + L1 sparsity | **None** -- arbitrarily keeps one of a correlated group | Sparse coefs | Unstable under collinearity; caps at n features when p>n |
| Embedded | Elastic net (Zou-Hastie 2005) | Prediction + L1+L2 grouping | Keeps correlated groups together | Sparse coefs | Two hyperparameters; still not "causal" |
| All-relevant | Boruta (Kursa 2010) | Every feature beating shadow features | Keeps all relevant (redundant included) | Confirmed/Tentative/Rejected | Slow; returns redundant sets by design |
| Meta / stability | Stability selection (Meinshausen 2010; Shah-Samworth 2013) | Selection probability under subsampling | Inherits base learner | Selection frequencies + threshold | Error bounds assume exchangeability omics violates |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Want every implicated gene for pathway/biology interpretation | Boruta (all-relevant), or stability-based consensus | Keeps whole correlated modules, not one representative |
| Want a small deployable assay/signature | Elastic-net (not bare LASSO); report stability | L2 grouping keeps correlated genes together and resamples more stably |
| p is huge (>20k); selection is slow | Univariate pre-filter to a few thousand, then Boruta/elastic-net, all inside the CV fold | Cheap dimensionality cut; never pre-filter on the full dataset |
| Need to report model performance | Wrap selection in a `Pipeline`, estimate by nested CV | Selection outside CV inflates AUC to ~perfect on pure noise |
| Single-cell biomarker across conditions | Pseudobulk per donor, then select at the donor level | The unit is the donor, not the cell (Squair 2021); cells are pseudoreplicates |
| Want to know which genes "drive" a trained model | -> machine-learning/prediction-explanation | SHAP ranking is not validated selection |
| Want unbiased accuracy/calibration of the selected model | -> machine-learning/model-validation | Selection is one step; validation is its own discipline |

## Leakage-Safe Selection (the single most damaging error to avoid)

**Goal:** Estimate the performance of a selection-plus-model pipeline without optimistic bias.

**Approach:** Put selection in a `Pipeline` so it is re-fit on each training fold only; the held-out fold never informs which features are kept. Selecting the top-k features on the *whole* dataset before cross-validating the classifier produces near-zero apparent error even on pure noise (Ambroise-McLachlan 2002). Selection is where almost all overfitting capacity lives when p>>n.

```python
from sklearn.pipeline import Pipeline
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold

pipe = Pipeline([
    ('select', SelectKBest(f_classif, k=20)),               # re-fit per fold -> no leakage
    ('clf', LogisticRegression(penalty='l2', max_iter=5000)),
])
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=0)
auc = cross_val_score(pipe, X, y, cv=cv, scoring='roc_auc')   # honest estimate
print(f'Nested-safe AUC: {auc.mean():.3f} +/- {auc.std():.3f}')
```

The standalone Boruta/LASSO blocks below select features on a full matrix to *discover* candidates; that is fine for discovery, but any performance number must come from the Pipeline pattern above, with selection inside the fold.

## All-Relevant: Boruta

**Goal:** Enumerate every fe