# Bio Machine Learning Biomarker Discovery

> Selects biomarker features from high-dimensional omics data using Boruta all-relevant selection, mRMR, LASSO/elastic-net, and stability selection, while controlling the leakage, irreproducibility, and correlated-feature traps that make most published signatures fail to replicate. Use when identifying candidate biomarkers, deciding between an all-relevant and a minimal-optimal selector, or judging whether a selected gene set is reproducible. For unbiased performance estimation of the resulting model see machine-learning/model-validation; for interpreting a trained model see machine-learning/prediction-explanation.

- Skill: `pku-yuangroup/bio-machine-learning-biomarker-discovery` (Agent Skill, multi-file: 4 files)
- Install (CLI): `npx skillmds@latest add pku-yuangroup/bio-machine-learning-biomarker-discovery`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pku-yuangroup/bio-machine-learning-biomarker-discovery/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Product & Planning
- Author: pku-yuangroup (https://skillmd.com/u/pku-yuangroup)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pku-yuangroup/bio-machine-learning-biomarker-discovery

---


## 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 feature carrying signal, including redundant co-expressed genes.

**Approach:** Compare each real feature's importance to the maximum importance of permuted "shadow" features over many iterations; confirm features that consistently beat the best shadow.

```python
from boruta import BorutaPy
from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier(n_estimators=100, n_jobs=-1, class_weight='balanced', max_depth=5, random_state=42)
# perc=100 uses the max shadow importance (strict); two_step (default True) controls the multiple-testing correction.
boruta = BorutaPy(rf, n_estimators='auto', perc=100, two_step=True, max_iter=100, random_state=42)
boruta.fit(X.values, y.values)                              # numpy arrays, not pandas

confirmed = X.columns[boruta.support_]                      # all-relevant set (redundant by design)
tentative = X.columns[boruta.support_weak_]
```

## Minimal-Optimal: Elastic Net (prefer over bare LASSO)

**Goal:** A small, stable predictive signature from correlated omics features.

**Approach:** Use elastic net, whose L2 term induces a grouping effect so correlated genes enter or leave together; standardize first because the penalty is scale-sensitive. Bare LASSO keeps one arbitrary member of a correlated group and flips on tiny data perturbations.

```python
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler

X_scaled = StandardScaler().fit_transform(X)               # for real scoring, do this inside the Pipeline
# saga is the only solver supporting elasticnet; C = 1/lambda (opposite of alpha in Lasso/ElasticNet).
enet = LogisticRegressionCV(penalty='elasticnet', solver='saga',
                            l1_ratios=[0.1, 0.5, 0.9], Cs=20, cv=10, max_iter=10000)
enet.fit(X_scaled, y)
selected = X.columns[enet.coef_[0] != 0]
```

## Stability: Are the Selected Features Reproducible?

**Goal:** Distinguish a robust signature from a resampling accident, and report stability alongside accuracy.

**Approach:** Run the selector on many subsamples, count per-feature selection frequency, keep features above a threshold (0.6 is the common default), and compute a chance-corrected stability index. Use the Nogueira 2018 measure (handles variable-size selections, gives a confidence interval); the older Kuncheva index needs equal-size subsets and breaks for LASSO.

```python
import numpy as np
from sklearn.linear_model import LogisticRegression

n_subsample, p = 100, X.shape[1]
counts = np.zeros(p)
subsets = []
for _ in range(n_subsample):
    idx = np.random.choice(len(X), size=len(X) // 2, replace=False)   # n/2 subsampling
    fit = LogisticRegression(penalty='l1', solver='liblinear', C=0.1, max_iter=2000).fit(X.iloc[idx], y.iloc[idx])
    mask = fit.coef_[0] != 0
    counts += mask
    subsets.append(mask.astype(int))

stable = X.columns[counts / n_subsample > 0.6]             # pi_thr=0.6: Meinshausen-Buhlmann default
# Nogueira stability index (chance-corrected; 1 = identical selections, ~0 = random):
Z = np.array(subsets); pbar = Z.mean(axis=0); k = Z.sum(axis=1)
stability = 1 - (Z.var(axis=0, ddof=1).mean()) / ((k.mean() / p) * (1 - k.mean() / p))
print(f'{len(stable)} stable features; Nogueira stability = {stability:.2f}')
```

## Per-Method Failure Modes

### Interpreting minimal-optimal membership as biology
- **Trigger:** Reporting "LASSO selected gene X but not its co-expressed partner Y" as a biological finding.
- **Mechanism:** L1 geometry keeps one vertex of a correlated group arbitrarily; the choice flips across resamples.
- **Symptom:** Selected genes change completely on a different train/test split though accuracy is stable.
- **Fix:** Use elastic net (grouping effect) or report selection *frequencies*; never read membership as importance ordering.

### Selection-before-CV leakage
- **Trigger:** Pick top-k features on all samples, then cross-validate the classifier on those features.
- **Mechanism:** The held-out folds informed which genes were kept; selection is the dominant overfitting capacity in p>>n.
- **Symptom:** Near-perfect CV accuracy, even reproducible on label-permuted (null) data; collapse on an external cohort.
- **Fix:** Selection lives inside the CV fold (Pipeline); estimate by nested CV (machine-learning/model-validation).

### Significance against the wrong null
- **Trigger:** Concluding a signature is real because it significantly predicts outcome.
- **Mechanism:** Random gene sets clear that bar; the transcriptome's proliferation axis is captured by almost any large set (Venet 2011).
- **Symptom:** The signature does not beat a size-matched random signature or a proliferation meta-gene in independent data.
- **Fix:** Benchmark against random-signature and proliferation-meta-gene nulls; require added value over clinical covariates in an *independent* cohort.

### Winner's curse / inflated effect sizes
- **Trigger:** Estimating effect sizes or AUC on the same data used to select features.
- **Mechanism:** Selected features are disproportionately those whose noise inflated their apparent effect (Goring 2001); the inflation can be near-total for small true effects.
- **Symptom:** Discovery AUC much higher than replication; replication is under-powered because it was sized to the inflated effect.
- **Fix:** Estimate effects on an independent split (cross-fitting / data-splitting); size replication for the shrunken effect.

### Pseudoreplication in single-cell selection
- **Trigger:** Treating thousands of cells from a few donors as independent samples.
- **Mechanism:** Cells within a donor are correlated; the effective n is the number of donors.
- **Symptom:** Grossly inflated significance and false discoveries.
- **Fix:** Pseudobulk per donor, select at the donor level (Squair 2021); confront the small true n.

## Reconciliation: When Methods Disagree

| Pattern | Likely cause | Action |
|---------|--------------|--------|
| Boruta keeps 200 genes, LASSO keeps 12 | All-relevant vs minimal-optimal answering different questions | Both can be right; pick by goal, do not "average" them |
| A list barely overlaps a published signature | Many disjoint equally-predictive lists exist (Ein-Dor 2005) | Expected, not a contradiction; compare *performance* and stability, not membership |
| High accuracy, low stability index | Resampling accident exploiting a dominant axis | Distrust the specific genes; prefer the lower-accuracy higher-stability candidate |
| FDR-clean list still fails to replicate | FDR controls testing, not selection stability | They are orthogonal; add stability + independent validation |

## Quantitative Thresholds

| Threshold | Source | Rationale |
|-----------|--------|-----------|
| Samples for a *stable* gene list ~ thousands | Ein-Dor 2006 | Small effects need large n for reproducible membership (accuracy needs far fewer) |
| Selection inside every CV fold; nested CV for tuning | Ambroise 2002; Simon 2003 | Selection outside CV gives ~0% error on noise |
| Stability threshold pi_thr ~ 0.6-0.9 | Meinshausen-Buhlmann 2010 | Selection-frequency cutoff; tune to false-positive cost |
| Random-signature null | Venet 2011 | Benchmark against size-matched random sets + proliferation meta-gene |
| Single-cell unit = donor (pseudobulk) | Squair 2021 | Cells are pseudoreplicates |
| Biomarker clinical translation rate <1% | Kern 2012 | Sets expectations; failures follow a foreseeable taxonomy |

## Common Errors

| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| BorutaPy raises on `np.float`/pandas input | Newer numpy removed aliases; needs arrays | Pass `X.values`, `y.values`; pin numpy or use a fork |
| Regularization strength backwards | `C=1/lambda` (logistic) vs `alpha` (Lasso/ElasticNet) are opposite conventions | Verify which API; small C = strong shrinkage |
| `elasticnet` penalty errors | Only `solver='saga'` supports it | Set `solver='saga'`, pass `l1_ratio(s)` |
| `mrmr_classif` returns wrong type | Pandas backend needs a DataFrame X and Series y | Pass `X` DataFrame, `y=pd.Series(y)`; K must still be chosen |
| glmnet signature unstable across runs (R) | Used `lambda.min` | Use `lambda.1se` for a sparser, more reproducible set |

## References

- Tibshirani R. 1996. Regression shrinkage and selection via the lasso. *J R Stat Soc B* 58:267-288.
- Goring HHH, Terwilliger JD, Blangero J. 2001. Large upward bias in estimation of locus-specific effects from genomewide scans. *Am J Hum Genet* 69:1357-1369.
- Ambroise C, McLachlan GJ. 2002. Selection bias in gene extraction on the basis of microarray gene-expression data. *PNAS* 99:6562-6566.
- Simon R, Radmacher MD, Dobbin K, McShane LM. 2003. Pitfalls in the use of DNA microarray data for diagnostic and prognostic classification. *J Natl Cancer Inst* 95:14-18.
- Ein-Dor L, Kela I, Getz G, Givol D, Domany E. 2005. Outcome signature genes in breast cancer: is there a unique set? *Bioinformatics* 21:171-178.
- Peng H, Long F, Ding C. 2005. Feature selection based on mutual information: criteria of max-dependency, max-relevance, and min-redundancy. *IEEE Trans Pattern Anal Mach Intell* 27:1226-1238.
- Zou H, Hastie T. 2005. Regularization and variable selection via the elastic net. *J R Stat Soc B* 67:301-320.
- Ein-Dor L, Zuk O, Domany E. 2006. Thousands of samples are needed to generate a robust gene list for predicting outcome in cancer. *PNAS* 103:5923-5928.
- Kursa MB, Rudnicki WR. 2010. Feature selection with the Boruta package. *J Stat Softw* 36:1-13.
- Meinshausen N, Buhlmann P. 2010. Stability selection. *J R Stat Soc B* 72:417-473.
- Venet D, Dumont JE, Detours V. 2011. Most random gene expression signatures are significantly associated with breast cancer outcome. *PLoS Comput Biol* 7:e1002240.
- Kern SE. 2012. Why your new cancer biomarker may never work: recurrent patterns and remarkable diversity in biomarker failures. *Cancer Res* 72:6097-6101.
- Shah RD, Samworth RJ. 2013. Variable selection with error control: another look at stability selection. *J R Stat Soc B* 75:55-80.
- Nogueira S, Sechidis K, Brown G. 2018. On the stability of feature selection algorithms. *J Mach Learn Res* 18:1-54.
- Squair JW, Gautier M, Kathe C, et al. 2021. Confronting false discoveries in single-cell differential expression. *Nat Commun* 12:5692.

## Related Skills

- machine-learning/model-validation - Nested CV and leakage-safe estimation of the selected model
- machine-learning/prediction-explanation - Why SHAP rankings are not a validated selection method
- machine-learning/omics-classifiers - Build a classifier from the selected features
- differential-expression/de-results - Pre-filter candidates with differential expression
- experimental-design/multiple-testing - FDR control and why it is orthogonal to selection stability
- experimental-design/power-analysis - Sample size for a stable signature vs an accurate predictor
- pathway-analysis/go-enrichment - Functional enrichment of an all-relevant gene set