Biomarker discovery analysis: stability-scored signatures from nested cross-validation

Key facts

What is biomarker discovery?

Biomarker discovery uses supervised machine learning on high-dimensional feature matrices to find a small multivariate panel that separates outcome groups or ranks risk [1]. Univariate differential testing plus a top-N cutoff ranks features one at a time. Discovery instead optimizes correlated feature combinations under regularization or wrapper selection inside cross-validation [2]. Mirza et al. [3] applied seven ML methods to 701 breast cancer samples across 11 GEO datasets, yielding nine-gene diagnostic and eight-gene prognostic signatures validated by qRT-PCR. Pepkio starts from client-supplied or Pepkio-processed matrices. Custom workflows are agreed at kickoff. See the biomarker discovery glossary.

When should you use biomarker discovery?

Use biomarker discovery when labeled samples provide a high-dimensional feature matrix and the question calls for a compact multivariate signature. Univariate ranking alone, or a pre-specified analyte panel, usually does not need this workflow [1].

• Neoadjuvant breast cancer response: Sammut et al. [7] built a multi-omic ML predictor from 168 neoadjuvant patients; the integrated model achieved AUC 0.87 on a 75-patient external cohort.
• Breast cancer transcriptomic signatures: Mirza et al. [3] derived nine-gene diagnostic and eight-gene prognostic signatures from 701 samples across 11 GEO datasets, with qRT-PCR validation.
• HCC prognostic panel: Feng et al. [8] used 77 algorithms on scRNA-seq and bulk RNA-seq to derive an 11-gene NK-related signature validated across TCGA, GEO, and ICGC cohorts.

How the analysis works — step by step

1. 1. Scope outcome, signature target, and validation plan

   Pepkio confirms the prediction target, endpoint type (binary, multiclass, or survival when scoped), signature size range, and whether an independent holdout or external cohort exists [1, 6]. Confounded batch-by-outcome designs and duplicate patients are flagged before modeling.

   Tools and outputs

   Tools used: Custom scoping templates; DOME checklist [6]

   Output: Signed scope memo with primary metrics (AUROC, AUPRC, balanced accuracy, or C-index when survival is scoped)

2. 2. Validate feature matrix and harmonize metadata

   Sample IDs are matched between the feature matrix and outcome table. Orientation, missingness, duplicate features, and class balance are audited. Feature identifiers are mapped to GENCODE v44 gene IDs or UniProt accessions when annotation is in scope.

   Tools and outputs

   Tools used: Custom validation scripts; clusterProfiler 4.20.0 for ID mapping when scoped

   Output: sample_manifest.csv; feature_qc_summary.csv

3. 3. Profile missingness and unsupervised variance

   Missing values, near-zero-variance features, and outlier samples are summarized. Unsupervised variance filtering is documented but not applied using outcome labels. Supervised filtering before cross-validation inflates performance estimates [4].

   Tools and outputs

   Tools used: scikit-learn 1.6.1 (VarianceThreshold); custom R diagnostics

   Output: missingness_profile.csv; variance_filter_log.csv; sample and feature QC plots (PDF/SVG)

4. 4. Audit data leakage

   Pepkio applies the seven guiding questions from Bernett et al. [9] before splitting to confirm preprocessing, feature selection, and tuning occur only on training folds.

   Tools and outputs

   Tools used: Documented leakage checklist [9]

   Output: leakage_audit.md

5. 5. Define patient-level nested cross-validation splits

   Patient-level nested CV separates inner-loop feature selection and tuning from outer-loop performance estimation [4, 6]. Stratification preserves class proportions when sample size permits.

   Tools and outputs

   Tools used: scikit-learn 1.6.1 (StratifiedKFold, GridSearchCV); nestedcv (R) when client prefers R-native nested CV [5]

   Output: cv_fold_assignments.csv with outer and inner fold IDs per sample

6. 6. Preprocess within inner folds

   Scaling, imputation, and optional variance filtering are fit on inner training data only. Batch is regressed with limma removeBatchEffect on training folds when batch is not confounded with outcome [10].

   Tools and outputs

   Tools used: scikit-learn 1.6.1 (Pipeline, ColumnTransformer); limma 3.60

   Output: Fold-specific preprocessed matrices logged in pipeline parameters

7. 7. Run embedded feature selection in inner loops

   Feature selection runs inside inner CV folds, not on pooled train and test data [4]. Default: glmnet LASSO or elastic-net (glmnet 4.1-8); Boruta 11.0.0 when all-relevant identification is scoped; optional univariate prefilter inside inner folds only [2].

   Tools and outputs

   Tools used: glmnet 4.1-8; Boruta 11.0.0; scikit-learn 1.6.1

   Output: Per-fold selected feature lists; inner-CV tuning metrics

8. 8. Train classifiers and evaluate on outer folds

   A regularized logistic or elastic-net classifier is trained on inner-selected features and evaluated on outer-fold samples. XGBoost 2.1.3 is compared when sample size and scope permit [11]. Simpler baselines are reported alongside more complex models.

   Tools and outputs

   Tools used: scikit-learn 1.6.1; XGBoost 2.1.3; scikit-survival 0.24.1 when survival endpoints scoped

   Output: signature_performance_summary.csv (AUROC, AUPRC, balanced accuracy, or C-index per outer fold)

9. 9. Compute signature stability across folds

   Selection frequency, the proportion of outer folds in which each feature is selected, is computed for every candidate feature [2, 12]. A consensus signature includes features meeting a kickoff-agreed stability threshold (often ≥50-60% of outer folds).

   Tools and outputs

   Tools used: Custom R/Python stability scripts; scikit-learn 1.6.1

   Output: selected_features.csv (feature_id, selection_frequency, mean_coefficient, annotation); stability barplot and heatmap (PDF/SVG)

10. 10. Interpret signature biology and package deliverables

    Consensus features receive SHAP attribution and pathway enrichment (SHAP 0.46.0; clusterProfiler 4.20.0). Deliverables are assembled. External or holdout evaluation is reported separately when a cohort was reserved at kickoff [10].

    Tools and outputs

    Tools used: SHAP 0.46.0; clusterProfiler 4.20.0

    Output: pathway_enrichment_signature.csv; SHAP summary plots; final deliverable bundle; Methods draft

What Pepkio delivers

Technical decisions we make — and why

Validation design: nested CV vs single train/test split

Nested CV is the default because feature selection on pooled data inflates accuracy more than non-nested tuning [4], and DOME recommends nested CV for unbiased estimation [6]. A single holdout is used only when an external cohort is reserved at kickoff.

Feature selection: glmnet vs Boruta vs univariate filter

glmnet is the default for sparse signatures in p >> n settings [13]. Boruta identifies all relevant features when broader sets are needed [14]. Univariate prefilters run only inside inner folds. Haury et al. [2] found t-test filters competitive with wrappers on breast cancer data.

Stability threshold

Features meeting a kickoff-agreed selection frequency (often ≥50-60% of outer folds) enter the consensus panel [12].

Classifier: regularized logistic vs XGBoost

Regularized logistic or elastic-net is the default for small tabular cohorts [10, 11]. XGBoost is compared when sample size and scope support nonlinear boundaries.

Batch correction

limma removeBatchEffect on training folds when batch and outcome are separable [10]. Fully confounded designs are flagged at kickoff.

Common questions

Related services

• Predictive modeling — Full classifier or survival model validation, calibration assessment, and TRIPOD+AI-aligned reporting after a signature is defined.
• Multi-omics integration — Matched multi-layer feature matrices from MOFA2 or DIABLO as inputs for supervised signature construction.
• Bulk RNA-seq — Expression count matrices and differential-expression preprocessing as transcriptomic feature inputs.
• Olink proximity extension assay — Bridge-normalized NPX matrices and QC for targeted protein signature discovery.
• Bioinformatics consulting — Experimental design, power estimation, and covariate planning before committing to ML signature work.