# Bio Applied Data Harmonization

> Harmonize multi-omics data (RNA-seq, proteomics, methylation, metabolomics) before integration — per-layer normalization, KNN/half-minimum missing-value imputation, PCA/PVCA batch-effect detection, and ComBat correction with pandas/scikit-learn. Use when prepping matrices for MOFA2/DIABLO/mixOmics, fixing missing values in proteomics MS data, or removing batch effects confounded with sequencing date/site before joint analysis.

- Skill: `pavel-kravchenko/bio-applied-data-harmonization` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-data-harmonization`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-data-harmonization/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Data & Analytics
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Data type | Recommended normalization | Rationale |
|---|---|---|
| RNA-seq counts | TMM/DESeq2 VST -> z-score | Library size + biological variation |
| Proteomics intensities | Median centering + log2 -> z-score | Inter-sample loading variation |
| Methylation beta values | logit to M-value -> z-score | Beta is bounded [0,1], not normal |
| Metabolomics | Probabilistic quotient normalization + log | Corrects for dilution effects |

**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.

```python
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

| Mechanism | Description | Strategy |
|---|---|---|
| MCAR | Random technical failures | KNN or mean imputation |
| MAR | Missingness depends on observed values | KNN imputation |
| MNAR | Below the instrument's detection limit | Half-minimum (min/2) 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.

```python
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.

| Tool | Method | Use case |
|---|---|---|
| ComBat | Empirical Bayes | Single covariate batch correction (bulk RNA/proteomics) |
| ComBat-seq | Negative binomial | RNA-seq raw counts |
| Harmony | Iterative correction | scRNA-seq integration |
| limma `removeBatchEffect` | Linear model | After normalization, bulk data |
| MNN (mutual nearest neighbors) | Graph-based | Single-cell, multiple batches |

**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.

```python
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
```

```r
## 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:

```python
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)