Bio Applied Data Harmonization

Data Harmonization for Multi-Omics

When to Use

• Combining RNA-seq, proteomics, methylation, or metabolomics matrices from the same samples before joint modeling (MOFA2, DIABLO, iCluster, concatenation-based PCA).
• Proteomics or metabolomics data has 20-40% missing values and you need to decide MCAR/MAR/MNAR imputation strategy.
• Samples cluster by sequencing batch/run date on a PCA plot instead of by biological group.
• Not all samples have all omics layers (partial overlap) and you need to build a matched subset or handle it explicitly.
• Choosing a normalization scheme per data type (counts vs. intensities vs. beta values) before scaling for integration.

Version Compatibility

Python ≥3.10, scikit-learn ≥1.4, pandas ≥2.0, scipy ≥1.11, numpy ≥1.26. Concepts apply equally to MOFA2 ≥1.12 (R/Bioconductor) and mixOmics ≥6.28 workflows that consume the harmonized matrices produced here.

Prerequisites

• pip install scikit-learn pandas numpy scipy matplotlib seaborn
• Familiarity with PCA and basic statistics (variance, mean-centering).
• Each omics layer already QC'd and in samples x features orientation with a shared SampleID.

Integration Strategies (pick before harmonizing)

• Early (concatenation): merge all matrices column-wise -> single dimension reduction. Simple but lets high-dimensional layers (e.g., RNA) dominate.
• Late (meta-analysis): analyze each omics layer separately -> combine p-values/rankings. Loses cross-layer correlation structure.
• Intermediate (factor analysis): MOFA2, iCluster — learn shared latent factors across layers, robust to differing feature counts.
• Supervised (DIABLO/mixOmics): use sample labels to guide which cross-layer correlations are found.

Key principle: normalize and scale within each omics layer before any of the above — never z-score the concatenated matrix as a whole, since layers with more features would dominate the shared variance.

Normalization by Data Type

Goal: put each omics layer on a comparable scale so no single layer's feature count or variance magnitude dominates a joint PCA/factor model. Approach: apply the layer-appropriate transform first (log/logit/VST), then a generic scaler (z-score is the safe default across layers), and sanity-check with PCA that structure looks biological, not technical.

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler, QuantileTransformer
from sklearn.decomposition import PCA


def normalize_layer(matrix: np.ndarray, method: str = "zscore") -> np.ndarray:
    """Normalize one omics layer (samples x features).

    method: "zscore" (default, safe general choice) or
            "quantile" (forces identical distributions across features;
            can over-correct if batch is confounded with biology).
    """
    if method == "zscore":
        return StandardScaler().fit_transform(matrix)
    if method == "quantile":
        return QuantileTransformer(
            output_distribution="normal", random_state=42
        ).fit_transform(matrix)
    raise ValueError(f"unknown method: {method}")


def pca_variance_check(matrix: np.ndarray, n_components: int = 10) -> pd.Series:
    """Return % variance explained per PC — use to spot batch dominating PC1."""
    pca = PCA(n_components=n_components).fit(matrix)
    return pd.Series(
        pca.explained_variance_ratio_ * 100,
        index=[f"PC{i+1}" for i in range(n_components)],
    )

Missing Data Imputation

Proteomics missing values (typically 20-40%) are predominantly MNAR (low-abundance proteins fall below detection). Best practice:

1. Drop features (proteins/metabolites) missing in >50% of samples.
2. Half-minimum imputation for MNAR features; KNN for the remainder (MAR).
3. Flag every imputed cell (boolean mask) and carry it into downstream sensitivity analysis.

from sklearn.impute import KNNImputer


def impute_proteomics(matrix: np.ndarray, mnar_threshold: float = 0.5) -> np.ndarray:
    """Impute a proteomics/metabolomics matrix (samples x features) with NaNs.

    Drops features missing in more than `mnar_threshold` fraction of samples,
    then applies half-minimum imputation (MNAR-appropriate) to the rest.
    Use KNNImputer instead when missingness is closer to MAR (e.g. random
    dropout, not concentrated in low-intensity features).
    """
    miss_frac = np.isnan(matrix).mean(axis=0)
    keep = miss_frac <= mnar_threshold
    filtered = matrix[:, keep]

    col_min = np.nanmin(filtered, axis=0)
    imputed = np.where(np.isnan(filtered), col_min / 2, filtered)
    return imputed


def impute_knn(matrix: np.ndarray, n_neighbors: int = 5) -> np.ndarray:
    """KNN imputation for MAR missing data (e.g. random instrument dropouts)."""
    return KNNImputer(n_neighbors=n_neighbors).fit_transform(matrix)

Batch Effect Detection and Correction

Detection: PCA (samples cluster by batch, not biology), PVCA (variance-component decomposition), RLE plots.

Never correct batch when confounded with biology (e.g., all tumor samples in Batch1, all normal in Batch2 — correction would erase the biological signal). Always cross-tabulate pd.crosstab(batch, sample_type) first.

def combat_correct(data: np.ndarray, batch_labels: list[str]) -> np.ndarray:
    """Simplified mean-only ComBat-style batch correction.

    Shifts each batch's per-feature mean to the grand mean. This is a
    location-only approximation of ComBat's empirical-Bayes model — use the
    real `pycombat`/`sva::ComBat` for variance (scale) correction too.
    """
    corrected = data.copy()
    grand_mean = data.mean(axis=0)
    for b in set(batch_labels):
        mask = np.array(batch_labels) == b
        batch_mean = data[mask].mean(axis=0)
        corrected[mask] -= (batch_mean - grand_mean)
    return corrected

## ComBat-seq for RNA-seq raw counts (Bioconductor `sva` package, sva >= 3.50)
library(sva)

adjusted_counts <- ComBat_seq(
  counts = as.matrix(count_matrix),   # genes x samples, raw integer counts
  batch  = metadata$Batch,
  group  = metadata$SampleType         # preserve biological groups explicitly
)

Sample Matching Across Layers

When omics layers don't all have the same samples (common in real cohorts), decide up front whether to use only complete cases (all layers present) or a partial-overlap strategy per analysis:

def complete_case_samples(layer_sample_sets: dict[str, set]) -> set:
    """Intersect sample sets across omics layers to get complete cases.

    layer_sample_sets: e.g. {"rna": {...}, "prot": {...}, "meth": {...}}
    """
    all_sets = list(layer_sample_sets.values())
    return set.intersection(*all_sets)

Pitfalls

• Coordinate systems: BED is 0-based half-open; VCF/GFF are 1-based inclusive — mixing them causes off-by-one errors when joining interval-based annotations into a harmonized table.
• Batch effects: always check for batch/biology confounding (pd.crosstab) before correcting — correcting a confounded batch destroys the signal you're trying to detect.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when testing thousands of features simultaneously after harmonization.
• Normalizing the concatenated matrix: z-scoring after merging layers lets the layer with more features dominate variance — always normalize within each layer first.
• Quantile normalization over-correction: it forces identical distributions per feature, which can flatten real biological differences if batch and biology are entangled.

See Also

• bio-applied-mofa2 — factor-analysis integration (MOFA2) of harmonized layers
• bio-applied-mixomics — supervised integration (DIABLO) using sample labels
• bio-applied-sc-integration — single-cell-specific batch integration (Harmony, MNN)
• bio-applied-rna-seq-analysis — upstream RNA-seq normalization (TMM/DESeq2 VST)