Version Compatibility
Reference examples tested with: numpy 1.26+, pandas 2.2+, scikit-learn 1.4+, shap 0.47+ (the feature_perturbation='auto' estimand and per-class 3-D .values behavior the code relies on; xgboost 2.0+ optional).
Before using code patterns, verify installed versions match. If versions differ:
- Python:
pip show <package> then help(module.function) to check signatures
scikit-learn drift: CalibratedClassifierCV(cv='prefit') deprecated in 1.6 (use FrozenEstimator); LogisticRegression(penalty=) deprecated in 1.8, and LogisticRegressionCV(penalty='l1') too -- the 1.8+ migration drops penalty= entirely and passes l1_ratios=(1.0,) alone (leave penalty at its default; penalty='elasticnet' still emits the FutureWarning). XGBoost moved early_stopping_rounds to the constructor in 2.x. If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Biomarker Discovery Pipeline
"Build a validated biomarker panel from my omics data" -> Orchestrate group-aware splitting, feature selection, leakage-safe cross-validation, calibration, and SHAP interpretation to produce a robust, honestly-validated biomarker signature.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.
The governing principle
The whole pipeline stands or falls on four commitments made at the seams; each one, if broken, inflates the reported performance and a held-out set cannot detect the leak because it was already contaminated.
- The independent unit of splitting is the highest biological unit — patient/donor/site, NOT the sample — and it is committed first. Multiple biopsies, longitudinal samples, or technical replicates from one subject in both train and test is group leakage; the model memorizes the subject, not the biology. Split with
GroupKFold/StratifiedGroupKFold on a subject key. For single-cell-derived features the unit is the donor, not the cell.
- Every data-dependent transform is fit INSIDE the CV fold — scaling, library-size/quantile normalization, ComBat/SVA, PCA/UMAP, imputation, AND feature selection. The discovery panel may be selected on all training data (that IS the deliverable), but the performance NUMBER must come from a pipeline that re-runs selection per fold. Selection is the dominant overfitting capacity in p>>n and gives near-perfect apparent accuracy on pure noise (Ambroise & McLachlan 2002).
- The locked test set is touched exactly once. Every threshold, feature count, hyperparameter, and "best epoch" chosen on it leaks; when hyperparameters are tuned, use nested CV to report performance (Varma & Simon 2006).
- The metric is matched to the data regime, and calibration is separate from discrimination. AUC for discrimination, AUPRC/MCC when imbalanced, and Brier + a reliability curve whenever risk estimates will be used — AUC is invariant to any monotone score transform, so it says nothing about calibration.
Workflow Overview
Expression matrix + Metadata
|
v
[1. Data Preparation] -----> StandardScaler, train/test split
|
v
[2. Feature Selection] ----> Boruta or LASSO stability selection
|
v
[3. Model Training] -------> Pipeline with selection inside CV (leakage-safe)
|
v
[4. Model Interpretation] -> SHAP values, feature importance
|
v
[5. Validation] -----------> Hold-out test, bootstrap CI
|
v
Validated biomarker panel + classifier
Step 1: Data Preparation
Goal: Load the matrix and hold out a GROUP-aware test set before anything is fit.
Approach: Split by the subject key so no subject appears in both train and test, then fit the scaler on training only; per-fold scaling is re-applied inside the CV pipeline in Step 3.
import pandas as pd
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.preprocessing import StandardScaler
expr = pd.read_csv('expression.csv', index_col=0)
meta = pd.read_csv('metadata.csv', index_col=0)
X = expr.T
POSITIVE_CLASS = 'disease'
y = (meta.loc[X.index, 'condition'].values == POSITIVE_CLASS).astype(int)
groups = meta.loc[X.index, 'subject_id'].values
sgkf = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=42)
train_idx, test_idx = next(sgkf.split(X, y, groups))
X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
y_train, y_test = y[train_idx], y[test_idx]
groups_train = groups[train_idx]
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
QC Checkpoint 1: Check class balance, sample counts, and group separation
- Minimum 10 samples per class recommended; classes reasonably balanced (ratio <3:1)
- Confirm NO subject id appears in both train and test (
set(groups[train_idx]) & set(groups[test_idx]) is empty)
Step 2: Feature Selection
Goal: Produce the discovery panel (all-relevant with Boruta, or a stable minimal set with LASSO).
Approach: Optionally pre-filter, then run the selector and map the mask back to the full feature space for downstream indexing.
Option A: Boruta (All-Relevant Selection)
import numpy as np
from boruta import BorutaPy
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectKBest, f_classif
if X_train_scaled.shape[1] > 10000:
selector = SelectKBest(f_classif, k=5000)
selector.fit(X_train_scaled, y_train)
prefilter_idx = np.where(selector.get_support())[0]
X_train_filt = X_train_scaled[:, prefilter_idx]
else:
prefilter_idx = None
X_train_filt = X_train_scaled
rf = RandomForestClassifier(n_estimators=100, max_depth=5, n_jobs=-1, random_state=42)
boruta = BorutaPy(rf, n_estimators='auto', max_iter=100, random_state=42, verbose=0)
boruta.fit(X_train_filt, y_train)
selected_idx = np.zeros(X_train.shape[1], dtype=bool)
selected_idx[prefilter_idx[boruta.support_] if prefilter_idx is not None else boruta.support_] = True
print(f'Selected {selected_idx.sum()} features')
Option B: LASSO Stability Selection
from sklearn.linear_model import LogisticRegressionCV
import numpy as np
n_bootstrap = 100
stability_scores = np.zeros(X_train_scaled.shape[1])
for i in range(n_bootstrap):
idx = np.random.choice(len(y_train), size=len(y_train), replace=True)
model = LogisticRegressionCV(penalty='l1', solver='saga', Cs=10, cv=3, random_state=i, max_iter=1000)
model.fit(X_train_scaled[idx], y_train[idx])
stability_scores += (model.coef_[0] != 0).astype(int)
stability_scores /= n_bootstrap
selected_idx = stability_scores > 0.6
print(f'Selected {selected_idx.sum()} features (stability >0.6)')
QC Checkpoint 2:
- Selected features: 5-200 range
- Too few (<5): lower threshold, increase iterations
- Too many (>200): increase threshold, add pre-filtering
Step 3: Leakage-Safe Performance Estimation
Goal: Estimate performance without the selection-before-CV leakage that inflates AUC toward 1.0 even on noise.
Approach: The Step 2 selection produced the discovery panel (fit on all training data) -- that is fine for the final panel, but it must NOT be the data the performance number is computed on. Estimate performance with scaling and selection wrapped in a Pipeline so they re-fit inside each fold; for raw RNA-seq, do per-sample normalization outside the fold and gene scaling/selection inside it.
from sklearn.model_selection import StratifiedGroupKFold, cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
pipe = Pipeline([
('scaler', StandardScaler()),
('select', SelectKBest(f_classif, k=min(50, X_train.shape[1]))),
('clf', LogisticRegression(max_iter=5000, class_weight='balanced')),
])
outer_cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=42)
cv_scores = cross_val_score(pipe, X_train, y_train, groups=groups_train, cv=outer_cv, scoring='roc_auc')
print(f'Leakage-safe CV AUC: {cv_scores.mean():.3f} +/- {cv_scores.std():.3f}')
QC Checkpoint 3:
- AUC reported with its fold spread, not a bare number (small-n CV is high-variance)
- Confirm selection is inside the pipeline and folds are group-aware; selection-before-CV inflates AUC toward 1.0 even on noise
- For imbalanced data report AUPRC/MCC, not accuracy; check the model predicts biology not batch (machine-learning/omics-classifiers)
Step 4: Model Interpretation
Goal: Audit what the final model keys on, not select biomarkers.
Approach: Fit the final model on the discovery panel, then compute interventional SHAP against a background and aggregate over modules to catch shortcut/batch learning.
import shap
import numpy as np
from sklearn.ensemble import RandomForestClassifier
sel = X_train.columns[selected_idx]
clf = RandomForestClassifier(n_estimators=300, random_state=42, n_jobs=-1).fit(X_train[sel], y_train)
background = shap.utils.sample(X_train[sel], 100)
explainer = shap.TreeExplainer(clf, data=background, feature_perturbation='interventional')
shap_values = explainer(X_test[sel])
if shap_values.values.ndim == 3:
shap_values = shap_values[:, :, 1]
mean_shap = np.abs(shap_values.values).mean(axis=0)
QC Checkpoint 4:
- SHAP is an audit, not a selection method: use it to confirm the model is not keying on batch/housekeeping shortcuts (machine-learning/prediction-explanation)
- Aggregate SHAP over co-expression modules before ranking; within-module order is not a finding
- SHAP directions should be biologically plausible; treat top-SHAP genes as hypotheses, not a validated panel
Step 5: Final Validation -- Discrimination AND Calibration
Goal: Report honest held-out performance, including calibration when risks will be used.
Approach: Report discrimination with an interval, but if the panel will produce risk estimates, also check calibration: AUC is invariant to any monotone transform of the score, so a high AUC says nothing about whether the probabilities are honest (machine-learning/model-validation). External validation on an independent cohort is the real bar.
from sklearn.metrics import roc_auc_score
from sklearn.metrics import brier_score_loss
import numpy as np
y_prob = clf.predict_proba(X_test[sel])[:, 1]
test_auc = roc_auc_score(y_test, y_prob)
boot = []
for _ in range(1000):
i = np.random.choice(len(y_test), len(y_test), replace=True)
if len(np.unique(y_test[i])) == 2:
boot.append(roc_auc_score(y_test[i], y_prob[i]))
ci_lower, ci_upper = np.percentile(boot, [2.5, 97.5])
print(f'Hold-out AUC: {test_auc:.3f} 95% CI [{ci_lower:.3f}, {ci_upper:.3f}]')
print(f'Brier score (calibration + refinement): {brier_score_loss(y_test, y_prob):.3f}')
Parameter Recommendations
| Step | Parameter | Recommendation |
|---|
| Split | n_splits (StratifiedGroupKFold) | 5 -> ~0.2 held out; lower n_splits for a larger test fraction |
| Boruta | max_iter | 100 (sufficient), 200 if tentative features |
| LASSO | n_bootstrap | 100 (quick), 500 for publication |
| LASSO | stability_threshold | 0.6 (standard), 0.8 for strict |
| Leakage-safe CV | folds | 5 (standard), 10 for small datasets; selection inside each fold |
| RF | n_estimators | 100-500 |
| XGBoost | learning_rate | 0.1 (conservative) |
Common Errors
The leakage seams first (each silently inflates performance and a held-out set cannot detect it), then operational issues.
| Symptom | Cause | Fix |
|---|
| Near-perfect CV AUC that collapses on external data | Features selected on the full dataset, then only the classifier CV'd | Wrap selection INSIDE the pipeline so it re-fits per fold (Ambroise & McLachlan 2002) |
| Optimistic AUC despite in-fold selection | Scaler/ComBat/PCA/imputation fit on all data before the split | Fit every data-dependent transform inside the fold (Pipeline) |
| Great CV, poor real-world performance | Repeated subjects (biopsies/longitudinal/replicates) split across train/test | Split by subject with StratifiedGroupKFold; the unit is the donor, not the sample |
| Reported AUC higher than any real fold | Same CV used to tune hyperparameters AND report | Nest: inner CV tunes, outer CV reports (Varma & Simon 2006) |
| Good AUC but risk estimates are miscalibrated | Resampling (SMOTE/undersampling) for imbalance, or AUC used as the only metric | Report AUPRC/MCC + Brier; recalibrate on a disjoint fold; do not resample-then-report calibration |
| No features selected | Too strict threshold | Lower stability threshold, increase iterations |
| Too many features (>200) | Noisy data | Add pre-filtering, increase regularization |
| Low CV AUC (<0.6) | No signal, low power | Check data quality, add samples |
| High variance across folds | Small sample size | Repeated stratified k-fold with an interval (LOOCV is degenerate for AUC) |
| SHAP features differ from selected | Correlated features split credit; attribution describes the model | Aggregate over modules; do not expect SHAP to match selection |
Export Results
import pandas as pd
import joblib
feature_names = X_train.columns[selected_idx].tolist()
pd.DataFrame({'feature': feature_names}).to_csv('biomarker_panel.csv', index=False)
joblib.dump(clf, 'biomarker_classifier.joblib')
joblib.dump(scaler, 'feature_scaler.joblib')
Related Skills
- database-access/geo-data - Public expression cohorts for validation sets
- database-access/sra-data - Pull raw FASTQ for re-quantified validation cohorts
- database-access/uniprot-access - Protein-level features (sequence, GO terms, PTMs) for protein biomarkers
- machine-learning/biomarker-discovery - Detailed feature selection methods
- machine-learning/model-validation - Nested CV implementation details
- machine-learning/omics-classifiers - Classifier options and tuning
- machine-learning/prediction-explanation - SHAP and LIME interpretation
- differential-expression/de-results - Pre-filter with DE genes
- pathway-analysis/go-enrichment - Functional enrichment of biomarkers
References
- Ambroise C, McLachlan GJ (2002) Selection bias in gene extraction on the basis of microarray gene-expression data. PNAS 99:6562-6566. DOI 10.1073/pnas.102102699. (feature selection must be inside the CV fold.)
- Varma S, Simon R (2006) Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics 7:91. DOI 10.1186/1471-2105-7-91. (nested CV for unbiased performance.)
- Whalen S, Schreiber J, Noble WS, Pollard KS (2022) Navigating the pitfalls of applying machine learning in genomics. Nature Reviews Genetics 23:169-181. DOI 10.1038/s41576-021-00434-9. (genomics-specific leakage and distribution-shift pitfalls.)
- Kapoor S, Narayanan A (2023) Leakage and the reproducibility crisis in machine-learning-based science. Patterns 4:100804. DOI 10.1016/j.patter.2023.100804. (a taxonomy of leakage, including group leakage.)
