Biomarker Discovery with ML

Biomarkers are measurable molecular indicators (genes, proteins, metabolites, methylation marks) that associate with a biological state of clinical interest — most commonly response to a drug or overall survival. In drug development, predictive biomarkers guide patient selection for clinical trials and companion diagnostic development. Prognostic biomarkers inform treatment intensity decisions independent of therapy.

The fundamental challenge: omics datasets are high-dimensional (thousands of features) but low-sample (hundreds of patients), creating severe overfitting risk. Rigorous cross-validation and independent validation cohorts are non-negotiable.

Omics Data Sources: TCGA and GEO

import pandas as pd
import numpy as np

# TCGA data access via RTCGA bioconductor (use Python wrapper) or firehose
# GEO data access via GEOparse library

def load_geo_dataset(geo_accession: str, output_dir: str = "/tmp/geo") -> dict:
    """
    Download and parse a GEO dataset into expression matrix and clinical data.

    Args:
        geo_accession: GEO accession e.g. 'GSE12345'
        output_dir: Local cache directory

    Returns:
        Dict with 'expression' (samples × features) and 'clinical' DataFrames.
    """
    import GEOparse
    gse = GEOparse.get_GEO(geo=geo_accession, destdir=output_dir, silent=True)

    # Build expression matrix: rows=genes, cols=samples
    expression_frames = []
    for sample_name, sample in gse.gsms.items():
        sample_df = sample.table[["ID_REF", "VALUE"]].set_index("ID_REF")
        sample_df.columns = [sample_name]
        expression_frames.append(sample_df)

    expression = pd.concat(expression_frames, axis=1).T  # samples × genes
    expression.index.name = "sample_id"

    # Clinical metadata
    clinical_records = []
    for sample_name, sample in gse.gsms.items():
        meta = {k: v[0] if isinstance(v, list) and v else v
                for k, v in sample.metadata.items()}
        meta["sample_id"] = sample_name
        clinical_records.append(meta)
    clinical = pd.DataFrame(clinical_records).set_index("sample_id")

    return {"expression": expression, "clinical": clinical}

Feature Selection for Biomarker Identification

from sklearn.linear_model import LassoCV, ElasticNetCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import shap
from sklearn.ensemble import GradientBoostingClassifier


def lasso_feature_selection(
    X: pd.DataFrame,
    y: np.ndarray,
    cv: int = 5,
    alphas: np.ndarray = None,
) -> pd.Series:
    """
    LASSO with cross-validated alpha selection for omics feature selection.

    Particularly effective for p >> n (many genes, few samples).

    Args:
        X: Feature matrix (samples × genes/proteins), z-score normalised
        y: Binary response labels (0/1) or continuous outcome
        cv: Number of CV folds
        alphas: Alpha grid for cross-validation

    Returns:
        Series of non-zero LASSO coefficients (selected biomarkers) sorted by |coef|.
    """
    if alphas is None:
        alphas = np.logspace(-4, 1, 50)

    pipeline = Pipeline([
        ("scaler", StandardScaler()),
        ("lasso", LassoCV(alphas=alphas, cv=cv, max_iter=10000, n_jobs=-1)),
    ])
    pipeline.fit(X, y)

    lasso = pipeline.named_steps["lasso"]
    coefs = pd.Series(lasso.coef_, index=X.columns)
    selected = coefs[coefs != 0].sort_values(key=abs, ascending=False)

    return selected


def shap_feature_importance(
    X: pd.DataFrame,
    y: np.ndarray,
    n_estimators: int = 300,
) -> pd.DataFrame:
    """
    SHAP-based feature importance using gradient boosting classifier.

    SHAP values capture both main effects and interaction effects,
    providing global importance and per-sample explanations.

    Args:
        X: Feature matrix (samples × features), z-score normalised
        y: Binary class labels

    Returns:
        DataFrame with features ranked by mean |SHAP value|.
    """
    model = GradientBoostingClassifier(
        n_estimators=n_estimators,
        max_depth=3,
        learning_rate=0.05,
        subsample=0.8,
        random_state=42,
    )
    model.fit(X, y)

    explainer = shap.TreeExplainer(model)
    shap_values = explainer.shap_values(X)

    mean_abs_shap = np.abs(shap_values).mean(axis=0)
    importance_df = pd.DataFrame({
        "feature": X.columns,
        "mean_abs_shap": mean_abs_shap,
    }).sort_values("mean_abs_shap", ascending=False).reset_index(drop=True)

    return importance_df

Survival Analysis with Cox Proportional Hazards

from lifelines import CoxPHFitter, KaplanMeierFitter
from lifelines.statistics import logrank_test
import matplotlib.pyplot as plt


def fit_cox_biomarker(
    clinical_df: pd.DataFrame,
    biomarker_col: str,
    duration_col: str = "os_days",
    event_col: str = "os_event",
    covariates: list[str] = None,
) -> dict:
    """
    Fit a Cox PH model for a single biomarker with optional covariates.

    Args:
        clinical_df: DataFrame with survival time, event indicator, and biomarker
        biomarker_col: Column name of the candidate biomarker (continuous or binary)
        duration_col: Follow-up time column
        event_col: Event/death indicator (1 = event observed)
        covariates: List of additional covariate column names to adjust for

    Returns:
        Dict with hazard ratio, 95% CI, p-value, and concordance index.
    """
    cols = [duration_col, event_col, biomarker_col] + (covariates or [])
    df = clinical_df[cols].dropna()

    cph = CoxPHFitter(penalizer=0.1)  # L2 regularisation for stability
    cph.fit(df, duration_col=duration_col, event_col=event_col)

    summary = cph.summary
    biomarker_row = summary.loc[biomarker_col] if biomarker_col in summary.index else None

    return {
        "hazard_ratio": float(biomarker_row["exp(coef)"]) if biomarker_row is not None else None,
        "hr_lower_95": float(biomarker_row["exp(coef) lower 95%"]) if biomarker_row is not None else None,
        "hr_upper_95": float(biomarker_row["exp(coef) upper 95%"]) if biomarker_row is not None else None,
        "p_value": float(biomarker_row["p"]) if biomarker_row is not None else None,
        "concordance_index": float(cph.concordance_index_),
        "n_samples": len(df),
        "n_events": int(df[event_col].sum()),
    }


def kaplan_meier_stratified(
    clinical_df: pd.DataFrame,
    group_col: str,
    duration_col: str = "os_days",
    event_col: str = "os_event",
) -> dict:
    """
    Kaplan-Meier curves for two patient groups with log-rank test.

    Args:
        clinical_df: DataFrame with survival data and group assignment
        group_col: Binary group column (0/1 or 'high'/'low' expression)
        duration_col: Follow-up time in days
        event_col: Event indicator

    Returns:
        Dict with log-rank test p-value and group median survivals.
    """
    groups = clinical_df[group_col].unique()
    if len(groups) != 2:
        raise ValueError(f"Expected 2 groups, got {len(groups)}: {groups}")

    g0 = clinical_df[clinical_df[group_col] == groups[0]]
    g1 = clinical_df[clinical_df[group_col] == groups[1]]

    result = logrank_test(
        g0[duration_col], g1[duration_col],
        event_observed_A=g0[event_col], event_observed_B=g1[event_col],
    )

    kmf0 = KaplanMeierFitter().fit(g0[duration_col], g0[event_col])
    kmf1 = KaplanMeierFitter().fit(g1[duration_col], g1[event_col])

    return {
        "group_labels": [str(g) for g in groups],
        "p_value": float(result.p_value),
        "log_rank_statistic": float(result.test_statistic),
        "median_survival_group0_days": float(kmf0.median_survival_time_),
        "median_survival_group1_days": float(kmf1.median_survival_time_),
        "significant": result.p_value < 0.05,
    }

Cross-Validation Strategy for Omics Biomarkers

from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import warnings


def nested_cv_biomarker(
    X: pd.DataFrame,
    y: np.ndarray,
    k_features: int = 20,
    outer_folds: int = 5,
    inner_folds: int = 3,
) -> dict:
    """
    Nested cross-validation for unbiased biomarker model performance estimation.

    CRITICAL: Feature selection MUST happen inside the inner CV loop, not before.
    Pre-selecting features on the full dataset causes massive data leakage in
    high-dimensional omics settings (p >> n).

    Args:
        X: Expression matrix (samples × all_features) — NOT pre-filtered
        y: Binary outcome labels
        k_features: Number of features to select per fold
        outer_folds: Outer CV folds (for unbiased performance estimate)
        inner_folds: Inner CV folds (for hyperparameter selection)

    Returns:
        Dict with outer CV AUC scores and description of pipeline.
    """
    # Pipeline: scale → select → classify (ALL steps inside CV)
    pipeline = Pipeline([
        ("scaler", StandardScaler()),
        ("selector", SelectKBest(f_classif, k=k_features)),
        ("classifier", LogisticRegression(
            penalty="l2", C=1.0, solver="lbfgs", max_iter=1000, random_state=42
        )),
    ])

    outer_cv = StratifiedKFold(n_splits=outer_folds, shuffle=True, random_state=42)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        cv_results = cross_validate(
            pipeline, X, y,
            cv=outer_cv,
            scoring=["roc_auc", "average_precision"],
            return_train_score=True,
            n_jobs=-1,
        )

    return {
        "outer_cv_auc_mean": float(cv_results["test_roc_auc"].mean()),
        "outer_cv_auc_std": float(cv_results["test_roc_auc"].std()),
        "outer_cv_ap_mean": float(cv_results["test_average_precision"].mean()),
        "pipeline_note": "Feature selection is inside CV loop — no leakage",
        "n_features_selected_per_fold": k_features,
    }

Patient Stratification with NMF

from sklearn.decomposition import NMF
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score


def nmf_patient_stratification(
    expression_matrix: np.ndarray,
    n_subtypes: int = 3,
    n_init: int = 10,
    max_iter: int = 1000,
) -> dict:
    """
    Non-negative Matrix Factorisation for unsupervised patient subtype discovery.

    NMF decomposes X ≈ W·H where:
    - W (samples × k): Patient membership weights for each subtype
    - H (k × genes): Metagene signature for each subtype

    Patient subtype assigned by argmax(W, axis=1) — strongest subtype membership.

    Args:
        expression_matrix: Non-negative expression values (samples × genes)
        n_subtypes: Number of molecular subtypes to identify
        n_init: Number of random initialisations (best result selected)
        max_iter: Maximum NMF iterations per run

    Returns:
        Dict with subtype assignments, metagene signatures, and reconstruction error.
    """
    best_model = None
    best_error = np.inf

    for seed in range(n_init):
        model = NMF(
            n_components=n_subtypes,
            init="nndsvda",
            random_state=seed,
            max_iter=max_iter,
            alpha_W=0.1,
        )
        W = model.fit_transform(expression_matrix)
        if model.reconstruction_err_ < best_error:
            best_error = model.reconstruction_err_
            best_model = model
            best_W = W

    subtype_assignments = np.argmax(best_W, axis=1)
    sil_score = silhouette_score(
        expression_matrix, subtype_assignments, metric="euclidean", sample_size=500
    ) if len(np.unique(subtype_assignments)) > 1 else None

    return {
        "subtype_assignments": subtype_assignments.tolist(),
        "n_subtypes": n_subtypes,
        "subtype_counts": {i: int((subtype_assignments == i).sum()) for i in range(n_subtypes)},
        "reconstruction_error": float(best_error),
        "silhouette_score": float(sil_score) if sil_score is not None else None,
        "metagene_matrix_H": best_model.components_,  # Shape: (n_subtypes, n_genes)
    }

Key Takeaways

• Feature selection must happen inside the cross-validation loop. Pre-filtering features on the full dataset is the most common and severe form of data leakage in omics biomarker studies. Use nested CV with feature selection as a pipeline step.
• Cox PH model concordance index (C-index) is the primary metric for survival biomarkers. C-index >0.70 is generally considered clinically meaningful; >0.75 is strong for a single omics biomarker.
• NMF is preferable to k-means for omics patient stratification because the non-negativity constraint on expression values produces interpretable, additive metagene programmes rather than arbitrary cluster centroids.
• LASSO alpha selection via cross-validation is sample-efficient in the p>>n regime. For very small cohorts (<50 samples), elastic net (mixing L1+L2) is more stable than pure LASSO.
• Always validate discovered biomarkers in an independent cohort before reporting clinical utility. TCGA and GEO provide matched discovery/validation pairs for many tumour types.

Try it live

Explore a related model hands-on: open Molecule Screening Engine — Similarity Search Live and change its parameters while it runs. Everything stays in your browser.