Bio Workflows Biomarker Pipeline

Version Compatibility

Reference examples tested with: matplotlib 3.8+, numpy 1.26+, pandas 2.2+, scanpy 1.10+, scikit-learn 1.4+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

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 feature selection (Boruta/LASSO), nested cross-validation classifier training, and SHAP interpretation to produce a robust, validated biomarker signature.

Complete pipeline from expression data to validated biomarker panels with classifier.

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] -------> RandomForest/XGBoost with nested CV
    |
    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

import pandas as pd
from sklearn.model_selection import train_test_split
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 = meta.loc[X.index, 'condition'].values

# test_size=0.2: Standard 80/20 split; use 0.3 for <100 samples
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)

# 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 per group

• Minimum 10 samples per class recommended
• Classes should be reasonably balanced (ratio <3:1)

Step 2: Feature Selection

Option A: Boruta (All-Relevant Selection)

from boruta import BorutaPy
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectKBest, f_classif

# Pre-filter if >10k features
if X_train_scaled.shape[1] > 10000:
    selector = SelectKBest(f_classif, k=5000)
    selector.fit(X_train_scaled, y_train)
    X_train_filt = X_train_scaled[:, selector.get_support()]
    feature_mask = selector.get_support()
else:
    X_train_filt = X_train_scaled
    feature_mask = None

# 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)

selected_idx = boruta.support_
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: Model Training with Nested CV

from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.ensemble import RandomForestClassifier

X_train_sel = X_train_scaled[:, selected_idx]
X_test_sel = X_test_scaled[:, selected_idx]

# outer_cv=5: Standard for performance estimation
outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# n_estimators=100: Sufficient for most omics
clf = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)
cv_scores = cross_val_score(clf, X_train_sel, y_train, cv=outer_cv, scoring='roc_auc')

print(f'Nested CV AUC: {cv_scores.mean():.3f} +/- {cv_scores.std():.3f}')

QC Checkpoint 3:

• AUC >0.7 acceptable, >0.8 good
• Fold variance <0.1 (stable performance)
• Check for overfitting: train AUC should not be >>test AUC

Step 4: Model Interpretation

import shap
import matplotlib.pyplot as plt

clf.fit(X_train_sel, y_train)

# SHAP v0.47+: call explainer directly
explainer = shap.TreeExplainer(clf)
shap_values = explainer(X_train_sel)

# Beeswarm: shows importance AND direction
shap.plots.beeswarm(shap_values, max_display=20, show=False)
plt.tight_layout()
plt.savefig('shap_beeswarm.png', dpi=150, bbox_inches='tight')
plt.close()

# Extract top features
import numpy as np
mean_shap = np.abs(shap_values.values).mean(axis=0)
top_shap_idx = np.argsort(mean_shap)[-20:]

QC Checkpoint 4:

• Top 20 SHAP features should have >50% overlap with selected features
• SHAP directions should be biologically plausible

Step 5: Final Validation

from sklearn.metrics import roc_auc_score, classification_report
import numpy as np

y_prob = clf.predict_proba(X_test_sel)[:, 1]
test_auc = roc_auc_score(y_test, y_prob)
print(f'Hold-out test AUC: {test_auc:.3f}')

# Bootstrap CI for AUC
# n_bootstrap=1000: Standard for publication-quality CI
n_bootstrap = 1000
boot_aucs = []
for i in range(n_bootstrap):
    idx = np.random.choice(len(y_test), size=len(y_test), replace=True)
    boot_aucs.append(roc_auc_score(y_test[idx], y_prob[idx]))

ci_lower, ci_upper = np.percentile(boot_aucs, [2.5, 97.5])
print(f'95% CI: [{ci_lower:.3f}, {ci_upper:.3f}]')

print(classification_report(y_test, clf.predict(X_test_sel)))

Parameter Recommendations

Troubleshooting

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

Biomarker Pipeline Usage Guide

Overview

End-to-end workflow for biomarker discovery combining feature selection, model training with nested cross-validation, interpretation, and validation. Produces a validated biomarker panel with an accompanying classifier.

Prerequisites

pip install scikit-learn boruta shap xgboost pandas numpy matplotlib joblib

Input data:

• Expression matrix (genes x samples) as CSV
• Metadata file with sample IDs and condition/label column

Quick Start

Tell your AI agent what you want to do:

• "Build a biomarker classifier from my expression data"
• "Run the full biomarker discovery pipeline with nested CV"
• "Select features and train a validated classifier"
• "Create a biomarker panel for disease vs control classification"

Example Prompts

Basic Biomarker Discovery

"I have expression.csv and metadata.csv. Build a biomarker classifier for my disease vs control samples."

"Run biomarker discovery with LASSO stability selection and nested cross-validation."

With Specific Methods

"Use Boruta for feature selection and train a Random Forest classifier with SHAP interpretation."

"Build a minimal biomarker signature using LASSO with strict stability threshold (0.8)."

Validation Focus

"Create a validated biomarker panel with bootstrap confidence intervals for AUC."

"Train a classifier with nested CV and export the model for external validation."

What the Agent Will Do

1. Load and prepare data with stratified train/test split 2. Scale features (fit on training only to prevent leakage) 3. Select features using Boruta or LASSO with stability selection 4. Train classifier with nested CV for unbiased performance estimation 5. Generate SHAP plots for model interpretation 6. Validate on held-out test set with bootstrap confidence intervals 7. Export biomarker panel and trained model

Tips

• Start with at least 20 samples per class for reasonable statistical power
• Use Boruta for comprehensive biomarker panels (finds all relevant features)
• Use LASSO for minimal signatures (finds sparse feature sets)
• Always use nested CV to avoid overfitting bias in performance estimates
• Check that SHAP top features align with selected features (sanity check)
• Pre-filter with differential expression if starting with >10k features
• Consider biological validation with independent dataset or orthogonal assay
• Class imbalance >3:1 may require stratification or SMOTE

Related Skills

• database-access/geo-data - Public expression cohorts for independent validation
• database-access/sra-data - Pull raw FASTQ to build re-quantified validation cohorts
• database-access/uniprot-access - Protein-level features (sequence, GO, 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