Bio Applied Mofa2

MOFA2: Multi-Omics Factor Analysis

When to Use

• Integrating 2+ omics layers (RNA-seq, proteomics, ATAC/methylation, metabolomics) measured on the same samples into a shared low-dimensional representation.
• Separating "shared" biological signal (appears in all views) from "view-specific" signal (e.g., epigenetic-only or proteomic-only variation).
• Handling missing views/samples (e.g., proteomics missing for a subset of patients) -- MOFA2 natively supports this, unlike PCA on a concatenated matrix.
• Finding candidate genes/proteins/CpGs driving a factor, then running gene set enrichment on the top weights.
• User says "MOFA", "MOFA2", "multi-omics factor analysis", "shared vs specific variance", or "factor analysis across omics layers".

Version Compatibility

• mofapy2 >= 0.7 (Python >= 3.8), muon >= 0.1.6, anndata >= 0.10, mudata >= 0.2
• R alternative: MOFA2 Bioconductor package >= 1.12 (same model, R interface)
• HDF5 output format is shared between languages -- a model trained in Python can be loaded in R with MOFA2::load_model() and vice versa.

Prerequisites

• pip install mofapy2 muon (or BiocManager::install("MOFA2") in R)
• Each omics view pre-processed to a samples x features matrix, ideally variance-stabilized/log-normalized and roughly Gaussian (count data should be transformed first; see bio-differential-expression-deseq2-basics or bio-single-cell-preprocessing)
• Views share the same sample IDs (partial overlap is fine)

Core Workflow

Goal: Fit a MOFA2 model on matched omics views and get per-factor, per-view variance explained.

Approach: Build a MuData object from per-view AnnData, run muon.tl.mofa, then load the trained model to extract Z (factors) and variance decomposition.

import numpy as np
import anndata as ad
import muon as mu


def build_mudata(rna_df, prot_df, meth_df):
    """Assemble a MuData object from samples x features DataFrames sharing an index.

    Each df: rows=samples, cols=features, already normalized/scaled per view.
    """
    mdata = mu.MuData({
        "rna": ad.AnnData(rna_df.values, obs=rna_df[[]], var=rna_df.T[[]]),
        "prot": ad.AnnData(prot_df.values, obs=prot_df[[]], var=prot_df.T[[]]),
        "meth": ad.AnnData(meth_df.values, obs=meth_df[[]], var=meth_df.T[[]]),
    })
    mdata.update()
    return mdata


def run_mofa(mdata, n_factors=15, outfile="mofa_model.hdf5", seed=42):
    """Train a MOFA2 model on a MuData object and save it to HDF5.

    n_factors: start high (15-20) and prune inactive factors afterward.
    """
    mu.tl.mofa(
        mdata,
        n_factors=n_factors,
        outfile=outfile,
        seed=seed,
        convergence_mode="fast",
    )
    return mdata  # factors written to mdata.obsm["X_mofa"], weights in mdata.varm

For direct control over the generative model (long-format data, custom priors), use mofapy2 directly:

from mofapy2.run.entry_point import entry_point


def run_mofa_low_level(data_df, n_factors=10, outfile="mofa_model.hdf5"):
    """Train MOFA2 via mofapy2 entry_point on long-format data.

    data_df columns: sample, group, feature, view, value.
    """
    ent = entry_point()
    ent.set_data_options(scale_groups=False, scale_views=True)
    ent.set_data_df(data_df)
    ent.set_model_options(factors=n_factors, spikeslab_weights=True, ard_weights=True)
    ent.set_train_options(iter=1000, convergence_mode="fast", seed=42)
    ent.build()
    ent.run()
    ent.save(outfile)
    return ent

Goal: Decompose variance explained (R^2) per factor per view to distinguish shared vs. view-specific factors.

Approach: Read the trained model back with mofax (or parse the HDF5 directly) and tabulate R^2.

import pandas as pd


def variance_decomposition(model_path):
    """Load a MOFA2 HDF5 model and return an R2 table (views x factors).

    Requires: pip install mofax
    """
    import mofax as mfx

    m = mfx.mofa_model(model_path)
    r2 = m.get_r2()  # columns: Factor, View, R2, (Group)
    table = r2.pivot_table(index="View", columns="Factor", values="R2")
    m.close()
    return table


def flag_shared_vs_specific(r2_table, threshold=2.0):
    """Classify factors as shared (high R2 in >=2 views) or view-specific.

    threshold: minimum R2 (%) to count a view as "active" for a factor.
    """
    active = r2_table >= threshold
    n_active_views = active.sum(axis=0)
    return pd.DataFrame({
        "n_active_views": n_active_views,
        "classification": np.where(n_active_views >= 2, "shared", "view-specific"),
    })

Model Parameters

Factor Interpretation Workflow

1. Correlate factor scores (mdata.obsm["X_mofa"]) with sample metadata (type, clinical variables, batch).
2. Rank features by |weight| per view/factor -> candidate genes/proteins/CpGs.
3. Run gene set enrichment (e.g., KS test or GSEA) on ranked RNA weights for a factor of interest -> pathway interpretation. See bio-pathway-analysis-gsea.

MOFA2 vs PCA

Pitfalls

• Factor 1 depth artifact: check whether Factor 1 correlates with library size/sequencing depth rather than biology before interpreting it.
• Batch effects: if batch confounds biology, MOFA will capture batch as its own factor -- always correlate factors against known technical covariates.
• View scaling: always set scale_views=True; otherwise a high-dimensional view (e.g., 20k genes vs. 150 proteins) dominates the shared factors purely by feature count.
• Too few factors: under-specifying n_factors merges distinct biological and technical signals into one factor -- start high and prune post hoc.
• Non-Gaussian input: raw counts violate the default Gaussian likelihood; log/VST-transform RNA-seq and normalize methylation M-values first.

See Also

• bio-multi-omics-integration-mofa-integration
• bio-multi-omics-integration-mixomics-analysis
• bio-multi-omics-integration-data-harmonization
• bio-pathway-analysis-gsea