name: bio-machine-learning-biomarker-discovery description: 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. origin: openai4s category: bioskills/machine-learning metadata: tool_type: python primary_tool: boruta third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT
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.
"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
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.
This axis matters more than filter/wrapper/embedded. Choosing the wrong one is the most common conceptual error in applied biomarker papers.
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.
| 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 |
| 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 |
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.
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_]
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]
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}')
| 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 |
| 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 |
| 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 |