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

Decision Tree by Scenario

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.

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.

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.

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.

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

Quantitative Thresholds

Common Errors

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