# Bio Applied Mofa2

> Run MOFA2 (mofapy2/muon) to fuse RNA-seq, proteomics, methylation into latent factors; decompose per-view R2, interpret weights. Use for multi-omics integration, MOFA/MOFA2 analysis, or latent factor discovery.

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

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

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

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

| Parameter | Default | Notes |
|---|---|---|
| `n_factors` | 10-15 | Start at 15-20, prune factors with R^2 < 2% in every view |
| `convergence_mode` | "fast" | fast/medium/slow; use "slow" for final published results |
| `spikeslab_weights` | True | Sparse feature weights via ARD prior |
| `scale_views` | True | Normalize each view to unit variance so high-dimensional views don't dominate |

## 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

| Aspect | PCA | MOFA2 |
|---|---|---|
| Input | Single matrix | Multiple matrices (views) |
| Shared vs specific | Mixed together | Explicit per-view decomposition |
| Missing views/samples | Cannot handle | Native support |
| Sparse weights | No | Yes (ARD/spike-slab prior) |

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