🧬 Single-Cell Multi-Omics Integration

📖 What is Single-Cell Multi-Omics Integration?

Single-cell multi-omics integration combines measurements of different molecular layers (transcriptome, epigenome, proteome, spatial location) from the same or related cells to build comprehensive cellular maps. This integration is essential for understanding cell states, developmental trajectories, disease mechanisms, and therapeutic responses.

Why Multi-Omics Integration Matters

• Comprehensive Cell State Definition: Single modalities provide incomplete views; RNA tells what's transcribed, ATAC reveals accessible chromatin, proteins show functional output
• Regulatory Mechanism Discovery: Linking chromatin accessibility → transcription → protein abundance reveals gene regulatory networks and signaling cascades
• Batch Effect Correction: Harmonizing data across experiments, technologies, and labs enables atlas-scale analyses and meta-studies
• Missing Modality Imputation: Predicting unmeasured features (e.g., protein from RNA) reduces experimental costs while maintaining biological insights
• Spatial Context Integration: Combining molecular profiles with spatial locations reveals tissue architecture and cell-cell interactions
• Perturbation Response Modeling: Understanding how genetic or chemical perturbations affect multiple molecular layers simultaneously

Integration Paradigms (Fu, Shaliu, et al. Nature Methods, 2025), (Liu, Chunlei, et al. Nature Methods, 2025)

According to the benchmark papers, there are six major multi-omics integration paradigms:

Each paradigm addresses different data structures and analytical challenges in multi-omics analysis.

📅 Evolution Timeline: From Paired Measurements to Foundation Models

🔬 Method Taxonomy: Algorithmic Approaches

By Computational Framework

By Computational Profile

📄 Landmark Papers by Computational Framework (2015-2025)

🧬 Experimental Technologies (Foundation)

🔷 Variational Autoencoders & Probabilistic Models

🕸️ Graph-Based Methods

🔄 Optimal Transport Methods

🧠 Deep Learning & Neural Networks

🎯 Contrastive Learning

🤖 Foundation Models & Transformers

🌉 Mosaic & Bridge Integration

🔬 Perturbation & Response Prediction

📊 Benchmark & Review Papers

🔧 Methods Development & Innovation

📚 Comprehensive Method Comparison

By Integration Category & Performance

🎯 From Scientific Question to Method

34 landmark papers are a lot to hold at once. The table below compresses them back into the question that actually matters when starting a project: what am I asking, what object do I need as output, and which approach typically produces it?

💡 Practical Implementation Guide

Choosing the Right Method: Decision Framework

Common Pitfalls & Best Practices

• Data Quality First: Ensure appropriate cell calling, modality-specific QC, doublet handling, and contamination assessment before integration — avoid universal hard thresholds; QC cutoffs should reflect assay, sequencing depth, and dataset structure. Use standard pipelines (Scanpy, Seurat). Bad data → bad integration.
• Modality- and model-specific preprocessing: RNA, ATAC, and protein have different measurement distributions, and preprocessing must match both the assay and the downstream model. For example, WNN commonly uses log-normalized RNA, TF-IDF/LSI for ATAC, and CLR/DSB for ADT, whereas generative models such as totalVI/MultiVI expect their own specified count-level inputs rather than those pre-derived features.
• Evaluation Metrics Matter: Use multiple complementary metrics (ARI, NMI, ASW for clustering; kBET, iLISI for batch mixing). Single metrics can be misleading.
• Batch Effect vs Biological Variation: Over-correction removes biology; under-correction leaves technical artifacts. Check cell type separation alongside batch mixing.
• Held-Out Testing: Evaluate on held-out batches or cell types, not just held-out cells from training conditions. Tests true generalization.
• Baseline Comparisons: Always compare against simple baselines (PCA, naive concatenation, per-modality analysis). Complex methods should outperform simple ones.
• Computational Scalability: Test methods on your actual data size. Many methods fail silently or degrade performance on large datasets.
• Biological Validation: Computational integration should reveal biologically meaningful patterns. Validate key findings with marker genes, known cell types, or experimental follow-up.
• Hyperparameter Tuning: Default parameters often suboptimal. Use cross-validation or grid search for critical parameters (latent dimensions, regularization, learning rate).
• Reproducibility: Report random seeds, package versions, preprocessing steps, and hyperparameters. Provide code and (when possible) data.

Software Ecosystem & Tools

• Scanpy: Standard Python toolkit for preprocessing, analysis, visualization. Foundation for many integration methods.
• Seurat: R package with Seurat v3/v5 integration, WNN, UMAP. Industry standard for single-cell RNA-seq.
• scvi-tools: PyTorch-based framework for probabilistic models (scVI, totalVI, MultiVI, scANVI). GPU-accelerated.
• Muon: Python package for multimodal omics analysis. Integrates with Scanpy ecosystem.
• LIGER: R package for integrative non-negative matrix factorization. Good for cross-platform integration.
• Harmony: Fast batch correction algorithm. Works on PCA embeddings. Widely used for large-scale integration.
• ArchR: R package for scATAC-seq analysis (2021, Nat Genetics) — Seurat's epigenomic companion for peak calling, gene-score computation, and downstream integration with RNA.
• Pertpy: Python framework unifying perturbation analysis workflows (CPA, GEARS, and optimal-transport methods like CellOT behind one interface with standard evaluation metrics).
• MuData (2022, Genome Biology): The stable multi-modal AnnData container — paired modalities stored in one .h5mu file with linked obs/var, the format Muon (above) is built on.
• spatialdata (2024, Nat Methods): OME-Zarr-based container unifying points, polygons, images, and tables for spatial multi-omics across platforms.

⚖️ Side-by-Side: Combining Modalities — Weighted-Graph (WNN) vs Joint Generative (totalVI / MultiVI)

When RNA and a second modality (protein in CITE-seq, chromatin in 10x Multiome) are measured in the same cell, the integration question is when to combine them: fuse per-modality neighbor graphs after embedding each separately (WNN), or learn a single joint latent space from the raw counts of both (a joint VAE like totalVI / MultiVI)?

🛠️ Hands-On Practice

The steps below walk through the joint-VAE route on CITE-seq with totalVI — one shared latent from RNA + protein — and note the WNN (graph) alternative and the MultiVI variant for RNA + ATAC.

Environment & packages

scvi-tools provides totalVI/MultiVI; scanpy/muon/mudata hold the modalities and run the WNN graph route in Python (Seurat's WNN is the R equivalent).

conda create -n multiomics python=3.10 -y
conda activate multiomics
pip install scvi-tools scanpy muon mudata
# WNN graph route is in muon (Python) or Seurat (R)

Hardware. WNN runs on a laptop; totalVI/MultiVI training wants a GPU (minutes) but works on CPU (slower) for datasets up to ~100k cells.

Data structures & formats

• CITE-seq AnnData — RNA counts in adata.X, protein (ADT) counts in adata.obsm["protein_counts"]
• 10x Multiome MuData — mdata["rna"] (genes) + mdata["atac"] (peaks), the same cells in both
• Pairing requirement varies by method — WNN and totalVI need the same cells measured in every modality (paired). MultiVI is designed for mosaic RNA/ATAC data — paired, unpaired, or a mix of RNA/ATAC availability across cells — and can infer/impute the unobserved RNA or ATAC modality within that supported RNA–ATAC setting.
• Joint latent — written to adata.obsm["X_totalVI"] / ["X_multivi"] for neighbors + UMAP
• Batch key — adata.obs["batch"] (donor / lane) for the model to correct

Minimal code walkthrough

Fit totalVI on CITE-seq to get one joint latent from RNA + protein.

import scvi, scanpy as sc

# CITE-seq: RNA counts in adata.X, protein (ADT) counts in adata.obsm["protein_counts"]
adata = sc.read_h5ad("cite_seq.h5ad")

# totalVI jointly models RNA (negative binomial) and protein (with ambient background)
scvi.model.TOTALVI.setup_anndata(
    adata,
    protein_expression_obsm_key="protein_counts",
    batch_key="batch",
)
model = scvi.model.TOTALVI(adata)
model.train()

# One joint latent from BOTH modalities; denoised protein is available too
adata.obsm["X_totalVI"] = model.get_latent_representation()
sc.pp.neighbors(adata, use_rep="X_totalVI"); sc.tl.umap(adata)

Swap in the graph route, or MultiVI for RNA + ATAC:

# Graph-level alternative — WNN via muon (Python port of Seurat's WNN):
# import muon as mu
# mu.pp.neighbors(mdata)        # per-modality graphs -> one weighted joint graph
# mu.tl.umap(mdata)

# RNA + ATAC (10x Multiome) instead of CITE-seq -> MultiVI:
# amvi = scvi.data.organize_multiome_anndatas(adata_paired, adata_rna, adata_atac)
# MultiVI's batch_key must hold the modality-availability label ("paired"/"expression"/"accessibility")
# added by organize_multiome_anndatas -- this is required by the model's architecture, not a
# stand-in for donor/lane batch. Real technical batch goes in categorical_covariate_keys instead:
# scvi.model.MULTIVI.setup_anndata(amvi, batch_key="modality", categorical_covariate_keys=["donor"])
# mvi = scvi.model.MULTIVI(amvi, n_genes=n_g, n_regions=n_r); mvi.train()

Common pitfalls & tips

• Paired vs mosaic vs diagonal. WNN and totalVI require the same cells in every modality. MultiVI is built for mosaic data — paired, unpaired, or mixed RNA/ATAC availability across cells — so it does not require full pairing. For fully diagonal data (no cell overlap, no shared modality at all) use Seurat anchor transfer, GLUE, or LIGER instead.
• Never concatenate raw features and PCA. RNA counts, protein ADT, and sparse binary ATAC peaks live on wildly different scales — naive concatenation lets one modality dominate the embedding.
• Protein has ambient background. ADT counts include background binding; totalVI models it explicitly, so don't preprocess protein like RNA (use CLR for WNN, raw counts for totalVI).
• ATAC preprocessing is model-specific. The classical WNN/LSI pipeline needs peak calling then TF-IDF/binarization before use. Generative models have their own input contract instead: MultiVI takes ATAC as its own count/binary-accessibility likelihood term and should not be fed pre-computed LSI/TF-IDF features meant for the WNN route — check what each model expects rather than assuming one preprocessing recipe fits all.
• Modality-availability label vs. technical batch are not the same. For MultiVI, batch_key must hold the modality-availability indicator (paired/expression/accessibility) — that's required by the model's architecture, not a stand-in for donor/lane batch. Pass real technical batch effects via categorical_covariate_keys instead, so both are corrected.
• Inspect the WNN weights. The per-cell modality weights are interpretable; a modality with no signal in some cells should get a low weight — check them to catch a broken or empty modality.