GPTomics/bio-workflows-biomarker-pipeline

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.

1. 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.
2. 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).
3. 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).
4. 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  # samples x genes
# y must be 0/1: brier_score_loss and calibration_curve raise on string labels, and sklearn orders
# classes alphabetically -- for a case/control column that makes 'control' the positive class, so
# predict_proba[:, 1], the SHAP [:, :, 1] slice, and Brier all silently describe the wrong class.
# AUC is symmetric and will not expose the flip. Encode the disease class as 1 explicitly.
POSITIVE_CLASS = 'disease'
y = (meta.loc[X.index, 'condition'].values == POSITIVE_CLASS).astype(int)
# The critical key: the SUBJECT (patient/donor/site), not the sample. If truly one
# sample per subject, groups = X.index; otherwise it MUST be the subject id.
groups = meta.loc[X.index, 'subject_id'].values

# Group- AND class-aware hold-out: take one StratifiedGroupKFold fold as the test set so no
# subject spans train/test (train_test_split(stratify=y) alone would leak repeated subjects).
sgkf = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=42)   # 1/5 held out (~0.2)
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]

# Fit scaler on training only to prevent data leakage
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

# Pre-filter if >10k features. selected_idx is a positional boolean mask aligned to X_train.columns.
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

# max_depth=5: Shallow trees for stable importances
rf = RandomForestClassifier(n_estimators=100, max_depth=5, n_jobs=-1, random_state=42)
# max_iter=100: Usually sufficient; 200 if many tentative
boruta = BorutaPy(rf, n_estimators='auto', max_iter=100, random_state=42, verbose=0)
boruta.fit(X_train_filt, y_train)

# Map the (possibly pre-filtered) Boruta mask back onto the FULL feature space.
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: Quick; use 500 for publication
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)
    # Cs=10: 10 regularization values to search
    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
# stability_threshold=0.6: Standard; 0.8 for strict
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

# Selection lives INSIDE the pipeline -> re-fit per fold, no leakage. Use the unscaled X_train.
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('select', SelectKBest(f_classif, k=min(50, X_train.shape[1]))),
    ('clf', LogisticRegression(max_iter=5000, class_weight='balanced')),
])
# Group-aware outer CV: pass groups_train so no subject spans a fold boundary.
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}')
# If hyperparameters are tuned, wrap a GridSearchCV (inner group CV) as the pipeline's estimator
# and report the OUTER cross_val_score -- flat CV that both tunes and reports is optimistic (Varma & Simon 2006).

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

# Fit the FINAL model on the discovery panel for interpretation and deployment.
sel = X_train.columns[selected_idx]
clf = RandomForestClassifier(n_estimators=300, random_state=42, n_jobs=-1).fit(X_train[sel], y_train)

# Interventional SHAP ('what the model uses') needs a background; set feature_perturbation
# explicitly because the 0.47+ 'auto' default flips the estimand on whether data= is given.
background = shap.utils.sample(X_train[sel], 100)
explainer = shap.TreeExplainer(clf, data=background, feature_perturbation='interventional')
shap_values = explainer(X_test[sel])
# RF returns one output per class in shap 0.47+ (n_samples, n_features, n_classes); keep the positive class.
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)

# Bootstrap CI for AUC (1000 resamples). Skip single-class resamples (roc_auc_score is nan there).
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}')   # y must be 0/1-encoded; string labels raise unless pos_label is passed
# If risks will be used, recalibrate on a disjoint fold and report a reliability curve
# (machine-learning/model-validation); do not resample for imbalance -- it breaks calibration.

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

# Save biomarker panel
feature_names = X_train.columns[selected_idx].tolist()
pd.DataFrame({'feature': feature_names}).to_csv('biomarker_panel.csv', index=False)

# Save model and scaler for deployment
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.)