Bio Applied Mixomics

mixOmics: Supervised Multi-Omics Integration

When to Use

• Classifying samples (e.g. cancer subtypes) from one or more omics layers with known group labels
• Selecting a sparse, stable biomarker panel from high-dimensional omics (p >> n) via sPLS-DA
• Integrating 2+ omics blocks (RNA-seq, proteomics, methylation, metabolomics) supervised by a shared outcome with DIABLO
• Comparing single-omics vs multi-block predictive power, or reporting cross-validated BER/AUC for a multi-omics classifier
• Needing a Python prototype of PLS-DA/DIABLO logic before running the real R mixOmics package

Version Compatibility

• R mixOmics ≥ 6.28 (Bioconductor 3.18+), R ≥ 4.3
• Python prototypes: scikit-learn ≥ 1.3 (PLSRegression), numpy ≥ 1.24, pandas ≥ 2.0

Prerequisites

• R: BiocManager::install("mixOmics")
• Python: pip install scikit-learn pandas numpy matplotlib seaborn
• Data already normalized/batch-corrected per omics layer (see bio-multi-omics-integration-data-harmonization) and matched sample IDs across blocks
• Familiarity with cross-validation and one-hot encoding of class labels

Method Family

PLS-DA: Single-Omics Classification

Goal: classify samples (e.g. 3 cancer subtypes) from a single standardized omics matrix. Approach: one-hot encode labels as Y, fit PLS regression maximizing cov(X, Y), classify by argmax of predicted Y. Tune n_components (2-5) by CV to minimize balanced error rate (BER).

import numpy as np
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import balanced_accuracy_score

class PLSDA:
    """PLS-DA: PLS regression against a one-hot-encoded class matrix."""
    def __init__(self, n_components=2):
        self.pls = PLSRegression(n_components=n_components, scale=True)

    def fit(self, X, Y_onehot):
        self.pls.fit(X, Y_onehot)
        return self

    def transform(self, X):
        return self.pls.transform(X)[0]

    def predict_class(self, X):
        return np.argmax(self.pls.predict(X), axis=1)

def cv_plsda(X, y_onehot, y, n_comp=3, cv=5):
    """5-fold stratified CV; returns mean/std balanced accuracy."""
    skf = StratifiedKFold(n_splits=cv, shuffle=True, random_state=42)
    accs = []
    for train, test in skf.split(X, y):
        model = PLSDA(n_components=n_comp).fit(X[train], y_onehot[train])
        y_pred = model.predict_class(X[test])
        accs.append(balanced_accuracy_score(y[test], y_pred))
    return np.mean(accs), np.std(accs)

sPLS-DA: Sparse Feature Selection

Goal: classify AND select a small, stable biomarker panel when most features are noise. Approach: in real mixOmics, keepX sets features retained per component, tuned by CV to minimize BER. Bootstrap the fit to score feature selection stability (>50% = reliable, >80% = high-confidence biomarker).

# R mixOmics — real tuning + fit
library(mixOmics)
test.keepX <- c(5, 10, 20, 50, 100)
tune.res <- tune.splsda(X, Y, ncomp = 3, validation = "Mfold",
                         folds = 5, dist = "max.dist",
                         test.keepX = test.keepX, nrepeat = 50)
optimal.keepX <- tune.res$choice.keepX

final.splsda <- splsda(X, Y, ncomp = 3, keepX = optimal.keepX)
selectVar(final.splsda, comp = 1)$name   # selected feature names, comp 1

import numpy as np
from sklearn.cross_decomposition import PLSRegression

class SPLSDA:
    """Toy sPLS-DA: hard-thresholds PLS weights to keepX features per component,
    then deflates X before extracting the next component."""
    def __init__(self, n_components=3, keepX=50):
        self.n_comp = n_components
        self.keepX = keepX
        self.selected_features = {}

    def fit(self, X, Y):
        X_deflated = X.copy()
        for k in range(self.n_comp):
            pls = PLSRegression(n_components=1, scale=False)
            pls.fit(X_deflated, Y)
            weights = np.abs(pls.x_weights_[:, 0])
            top_idx = np.argsort(weights)[-self.keepX:]
            sparse_weights = np.zeros(X.shape[1])
            sparse_weights[top_idx] = weights[top_idx]
            self.selected_features[f"comp_{k+1}"] = top_idx
            t = X_deflated @ sparse_weights
            t = t / (t @ t + 1e-10)
            X_deflated = X_deflated - np.outer(t, X_deflated.T @ t)
        return self

Feature stability via bootstrap: refit on resampled rows, count how often each feature lands in selected_features, keep features selected in >50% of bootstraps as stable biomarkers.

DIABLO: Multi-Block Supervised Integration

Goal: integrate 2+ omics blocks (e.g. RNA-seq + proteomics + methylation) supervised by one outcome, exploiting complementary signal across blocks. Approach: define an (M×M) design matrix of expected between-block correlation (0 = independent, 0.1 = default weak coupling, 1 = maximize correlation), fit block-wise sPLS-DA jointly, and average/consensus predictions across blocks.

# R mixOmics — real DIABLO fit
library(mixOmics)
X <- list(rna = rna_scaled, prot = prot_scaled, meth = meth_scaled)
design <- matrix(0.1, ncol = 3, nrow = 3, dimnames = list(names(X), names(X)))
diag(design) <- 0

list.keepX <- list(rna = c(20, 20), prot = c(10, 10), meth = c(15, 15))
diablo.model <- block.splsda(X = X, Y = Y, ncomp = 2, design = design,
                              keepX = list.keepX)
perf.diablo <- perf(diablo.model, validation = "Mfold", folds = 5, nrepeat = 10)
perf.diablo$WeightedVote.error.rate   # per-block + consensus BER
selectVar(diablo.model, block = "rna", comp = 1)$rna$name

import numpy as np
from sklearn.cross_decomposition import PLSRegression

class DIABLO:
    """Simplified multi-block PLS-DA: fits one PLS-DA per omics block, then
    predicts by averaging per-block class probabilities (consensus vote)."""
    def __init__(self, n_components=2):
        self.n_comp = n_components
        self.block_models = {}

    def fit(self, X_list, view_names, Y):
        self.view_names = view_names
        for X, name in zip(X_list, view_names):
            pls = PLSRegression(n_components=self.n_comp, scale=True)
            pls.fit(X, Y)
            self.block_models[name] = pls
        return self

    def predict_class(self, X_list):
        preds = [self.block_models[name].predict(X)
                 for X, name in zip(X_list, self.view_names)]
        return np.argmax(np.mean(preds, axis=0), axis=1)

Model Evaluation

Use stratified Mfold CV (5-10 folds, 10-50 repeats) and report BER (1 - balanced accuracy) and macro ROC-AUC, not raw accuracy — omics cohorts are frequently class-imbalanced. Compare each single-omics PLS-DA against DIABLO to confirm multi-block integration actually improves over the best single view before reporting it as such.

Pitfalls

• Overfitting: PLS-DA overfits badly when p >> n; always use sPLS-DA with CV-tuned keepX, and prefer nested CV (outer loop for evaluation, inner for tuning)
• Batch effects: check for batch confounding before interpreting components as biological signal — a batch-driven component looks identical to a real one
• Balanced Error Rate: use BER/balanced accuracy, never raw accuracy, with imbalanced classes
• Data leakage: BER < 5% on a real omics cohort is suspicious — check that scaling/feature selection happens inside each CV fold, not before splitting
• Design matrix choice: setting DIABLO's off-diagonal design too high (near 1) forces blocks to correlate even when biologically unrelated, hurting classification; 0.1 is a safe default

See Also

• bio-multi-omics-integration-data-harmonization — normalize/align blocks before mixOmics
• bio-multi-omics-integration-mofa-integration — unsupervised alternative (no labels needed)
• bio-machine-learning-biomarker-discovery — general ML biomarker workflows
• bio-machine-learning-model-validation — nested CV and validation strategy details