Single Cell Multi Omics Data Harmonization

Multi-Omics Data Harmonization

Prepare your RNA-seq, proteomics, methylation, and other omics datasets for joint integration by applying per-assay normalization, cross-assay batch correction, feature ID alignment, and missing value handling.

This is Step 1 of the multi-omics integration pipeline — all downstream skills (MOFA, DIABLO, SNF) depend on clean, consistently scaled output from this step.

What it does

1. Loads each omics layer into a unified container (MultiAssayExperiment in R; MuData in Python)
2. Applies the correct normalization strategy per data type:
   • RNA-seq counts → VST (DESeq2)
   • Proteomics LFQ intensity → log2 + median centering
   • Methylation β-values → M-value transformation: log2(β / 1−β)
   • ATAC-seq peaks → log1p(CPM); miRNA → log2(CPM + 1)
3. Generates PCA plots before batch correction to visualize batch structure
4. Applies ComBat batch correction across assays, preserving biological condition signal
5. Generates PCA plots after correction to confirm batch removal
6. Maps protein UniProt IDs and methylation probe IDs to HGNC gene symbols via Ensembl BioMart
7. Filters features with > 30% missing values, then imputes remaining NAs with MinProb
8. Z-scores all features and exports as .rds, .h5mu, and .csv for all downstream tools

Why this exists

If you ask a general AI to "prepare my multi-omics data for integration," it will:

• Apply the same normalization to all data types (wrong — RNA counts need VST; methylation needs M-value transformation, not log2)
• Run ComBat without checking for batch–condition confounding, silently removing biological signal
• Skip feature ID alignment, leaving protein UniProt IDs unmapped to gene symbols
• Not generate PCA plots to verify batch correction worked
• Export data in a format incompatible with MOFA, DIABLO, or SNF

This skill encodes the correct methodological decisions:

• Uses VST for RNA, log2+median centering for protein, and M-value for methylation — each chosen for statistical properties of that data type
• Checks table(Batch, Condition) before ComBat to detect confounding
• Aligns all features to HGNC gene symbols via Ensembl BioMart for cross-omics compatibility
• Filters high-missingness features before imputation to avoid noise amplification
• Exports in both R (.rds) and Python (.h5mu) formats for full downstream flexibility

Reference Methods

Normalization:

Batch correction: ComBat (Johnson et al., 2007) — parametric empirical Bayes, robust for small sample sizes. Requires ≥ 2 samples per batch per condition.

Feature alignment: Ensembl BioMart REST API — maps UniProt accessions and Illumina probe IDs to HGNC gene symbols for cross-omics feature matching.

Missing value imputation: MinProb — samples replacement values from a low-intensity distribution, appropriate for proteomics data Missing Not At Random (MNAR) patterns.

Usage (R)

library(MultiAssayExperiment)
library(DESeq2)
library(sva)
library(biomaRt)
library(ggplot2)

# ── Step 1: Build unified container ──────────────────────────────────────────
rna         <- as.matrix(read.csv("rnaseq_counts.csv",    row.names = 1))
protein     <- as.matrix(read.csv("proteomics.csv",        row.names = 1))
methylation <- as.matrix(read.csv("methylation_beta.csv", row.names = 1))
sample_info <- read.csv("sample_metadata.csv",             row.names = 1)

mae <- MultiAssayExperiment(
  experiments = ExperimentList(
    RNA         = SummarizedExperiment(assays  = list(counts    = rna),
                                        colData = sample_info),
    Protein     = SummarizedExperiment(assays  = list(intensity = protein),
                                        colData = sample_info),
    Methylation = SummarizedExperiment(assays  = list(beta      = methylation),
                                        colData = sample_info)
  ),
  colData = sample_info
)
mae <- intersectColumns(mae)   # keep only samples present in ALL assays
cat("Samples in all assays:", ncol(mae[[1]]), "\n")

# ── Step 2: Normalize per assay ───────────────────────────────────────────────
dds     <- DESeqDataSet(mae[["RNA"]], design = ~1)
vst_rna <- assay(vst(dds, blind = TRUE))

log2_prot <- log2(assay(mae[["Protein"]]))
log2_prot[is.infinite(log2_prot)] <- NA
norm_prot <- sweep(log2_prot, 2,
                   apply(log2_prot, 2, median, na.rm = TRUE) -
                   median(log2_prot, na.rm = TRUE))

beta   <- pmin(pmax(assay(mae[["Methylation"]]), 0.001), 0.999)
m_vals <- log2(beta / (1 - beta))

# ── Step 3: PCA before correction ────────────────────────────────────────────
pca_plot <- function(mat, meta, title) {
  df <- as.data.frame(prcomp(t(mat))$x[, 1:2])
  df$Batch <- meta$Batch; df$Condition <- meta$Condition
  ggplot(df, aes(PC1, PC2, color = Batch, shape = Condition)) +
    geom_point(size = 3) + ggtitle(title) + theme_bw()
}
p_before <- pca_plot(vst_rna, colData(mae), "RNA — before correction")

# ── Step 4: Batch correction ──────────────────────────────────────────────────
# Always check confounding first
print(table(colData(mae)$Batch, colData(mae)$Condition))

batch          <- colData(mae)$Batch
mod            <- model.matrix(~ Condition, data = as.data.frame(colData(mae)))
corrected_rna  <- ComBat(vst_rna,   batch = batch, mod = mod)
corrected_prot <- ComBat(norm_prot, batch = batch, mod = mod)
corrected_meth <- ComBat(m_vals,    batch = batch, mod = mod)

p_after <- pca_plot(corrected_rna, colData(mae), "RNA — after correction")
ggsave("figures/pca_before_after.pdf",
       gridExtra::grid.arrange(p_before, p_after, ncol = 2),
       width = 12, height = 5)

# ── Step 5: Feature ID alignment ──────────────────────────────────────────────
ensembl <- useEnsembl("genes", dataset = "hsapiens_gene_ensembl")
id_map  <- getBM(attributes = c("uniprotswissprot", "hgnc_symbol"),
                  filters    = "uniprotswissprot",
                  values     = rownames(corrected_prot), mart = ensembl)
rownames(corrected_prot) <- id_map$hgnc_symbol[
  match(rownames(corrected_prot), id_map$uniprotswissprot)]
corrected_prot <- corrected_prot[!is.na(rownames(corrected_prot)), ]

# ── Step 6: Missing value handling ───────────────────────────────────────────
corrected_prot <- corrected_prot[rowMeans(is.na(corrected_prot)) < 0.30, ]

impute_minprob <- function(mat) {
  for (j in seq_len(ncol(mat))) {
    nas <- is.na(mat[, j])
    if (any(nas)) {
      q01 <- quantile(mat[, j], 0.01, na.rm = TRUE)
      mat[nas, j] <- rnorm(sum(nas), mean = q01, sd = abs(q01) * 0.1)
    }
  }; mat
}
corrected_prot <- impute_minprob(corrected_prot)

# ── Step 7: Scale and export ──────────────────────────────────────────────────
scale_mat  <- function(mat) t(scale(t(mat)))
harmonized <- list(
  RNA         = scale_mat(corrected_rna),
  Protein     = scale_mat(corrected_prot),
  Methylation = scale_mat(corrected_meth),
  sample_meta = as.data.frame(colData(mae))
)
saveRDS(harmonized, "harmonized_multiomics.rds")
write.csv(harmonized$RNA,         "harmonized_rna.csv")
write.csv(harmonized$Protein,     "harmonized_protein.csv")
write.csv(harmonized$Methylation, "harmonized_methylation.csv")

Usage (Python — muon)

import muon as mu
import scanpy as sc
import pandas as pd
import numpy as np
from combat.pycombat import pycombat

adata_rna  = sc.read_csv("rnaseq_counts.csv").T
adata_prot = sc.read_csv("proteomics.csv").T
meta       = pd.read_csv("sample_metadata.csv", index_col=0)

for a in [adata_rna, adata_prot]:
    a.obs = meta.loc[a.obs_names]

sc.pp.normalize_total(adata_rna, target_sum=1e6)
sc.pp.log1p(adata_rna)

X = np.log2(adata_prot.X + 1)
X -= np.median(X, axis=0) - np.median(X)
adata_prot.X = X

batch = meta.loc[adata_rna.obs_names, "Batch"]
for adata in [adata_rna, adata_prot]:
    df = pd.DataFrame(adata.X.T, columns=adata.obs_names,
                      index=adata.var_names)
    adata.X = pycombat(df, batch).T.values

mdata = mu.MuData({"rna": adata_rna, "prot": adata_prot})
mu.pp.intersect_obs(mdata)
mdata.write("harmonized_multiomics.h5mu")

Example Output

Multi-Omics Data Harmonization
===============================
Input assays:   RNA | Protein | Methylation
Common samples: 92 (of 97 total — 5 missing at least one assay)

Normalization applied:
  RNA:         VST (DESeq2)         — 18,432 features
  Protein:     log2 + median center — 3,847 features
  Methylation: M-value              — 11,209 features

Batch correction (ComBat, 2 batches):
  RNA  — PC1 batch variance: 34.2% → 2.1%   ✓ corrected
  Protein — PC1 batch variance: 28.7% → 3.4% ✓ corrected

Feature ID alignment:
  Protein: 3,847 UniProt IDs → 3,612 HGNC symbols
  (235 removed: unmapped or duplicate)

Missing value summary (protein):
  Filtered (>30% missing): 3,612 → 3,589 features retained
  Imputed (MinProb):        0% missing after imputation

Exported:
  harmonized_multiomics.rds
  harmonized_rna.csv         (18,432 × 92)
  harmonized_protein.csv     ( 3,589 × 92)
  harmonized_methylation.csv (11,209 × 92)
  figures/pca_before_after.pdf

Interpretation Guide

• PCA before correction: Samples should separate by Batch on PC1/PC2 if a batch effect is present — this confirms ComBat is necessary
• PCA after correction: Batch structure should disappear; samples should now separate by Condition instead
• If batch and condition are confounded (e.g., all disease samples in Batch 1): skip ComBat — it will remove the biological signal you want to keep
• Protein missingness > 30%: features exceeding this threshold are filtered, not imputed — imputing very sparse features introduces more noise than signal
• M-value vs β-value: M-values are statistically preferable for downstream modelling (homoscedastic); β-values (0–1 range) are more interpretable for visualization — keep both if needed

Citation

If you use this skill in a publication, please cite:

• Johnson, W.E. et al. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics, 8(1), 118–127.
• Ramos, M. et al. (2017). Software for the integration of multiomics experiments in Bioconductor. Cancer Research, 77(21), e39–e42. (MultiAssayExperiment)
• Love, M.I. et al. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550.