Integrate multi-omics layers into interpretable factors with MOFA+

Hand Claude Code two or more omics matrices measured on the same samples or cells (e.g., RNA + ATAC + protein, or transcriptome + methylation + CNV); get back a trained MOFA+ model, a variance-explained breakdown per factor per layer, factor–metadata associations, and the top feature loadings that make each factor interpretable.

   
Problem class Data analysis
Subject areas Molecular and Cellular Biology, Immunology and Microbiology, Translational Medicine
Evidence level Proposed
Complexity One skill or MCP
Availability Fully open
Compute Laptop

Problem

Multi-omics studies measure several molecular layers on one set of samples — transcriptome, chromatin accessibility, methylation, surface protein, copy number — and the analyst has to find the axes of variation shared across layers versus those private to one. Concatenating the matrices and running PCA mixes scales and lets the highest-variance layer dominate; analyzing each layer separately misses the cross-layer structure that is usually the point. MOFA+ solves this by fitting an unsupervised factor model that, like PCA, returns a low-dimensional latent space, but decomposes the variance per factor per layer — so you can see that Factor 1 is driven by RNA+methylation while Factor 3 is ATAC-only, and correlate each factor with sample metadata (subtype, survival, treatment). Solved looks like: a trained model, a variance-decomposition heatmap, a table of factors ranked by association with your outcome, and the top-loading features per factor as the biological handle — all reproducible from a committed script.

This is a different problem from batch integration. scVI integration aligns the same modality across batches; MOFA+ integrates different modalities of the same samples into shared factors. Reach for MOFA+ when the layers differ, not when the batches do.

  1. Assemble one view per omics layer, aligned on the same samples/cells. Each view is a (samples × features) matrix; features can differ across views but the sample axis must match. For single-cell multi-modal data (10x Multiome, CITE-seq), build a MuData object with the muon skill; for bulk multi-omics, a dict of AnnData views is enough (see the AnnData skill). Feature-select each view first (e.g., highly variable genes/peaks) — MOFA+ scales poorly with tens of thousands of features per view and factors get noisier.

  2. Install the MOFA+ skill. From the SciAgent-Skills collection:

    git clone https://github.com/jaechang-hits/SciAgent-Skills

    Then inside Claude Code run /plugin install sciagent-skills and confirm it appears under /plugin → Installed. The skill runs mofapy2 locally via Bash/Python; install its declared dependencies when prompted.

  3. Have the skill build the model, train, and write the artifact. Drive it to a versioned script — a minimal prompt:

    Use the mofaplus-multi-omics skill. Load these views (samples on
the shared axis):
  - rna:    data/rna_hvg.h5ad
  - atac:   data/atac_features.h5ad
  - protein: data/adt.h5ad
Build a MOFA+ model, scale each view, train with 15 factors and
the default convergence settings, seed=1. Write the workflow to a
committed script mofa_run.py that:
  - trains the model and saves it to results/mofa_model.hdf5
  - writes results/variance_explained.csv (factor x view)
  - writes results/factor_metadata_assoc.csv (Pearson/ANOVA of
    each factor vs the columns in obs: subtype, treatment, ...)
  - writes results/top_loadings_<view>.csv (top 30 features per
    factor per view)
  - saves the variance-decomposition heatmap to figures/.
Pin the environment in requirements.txt.

  4. Read the model as a diagnostic, not a black box. Drop factors that explain negligible variance across all views; a factor loading on a single view at high variance often flags a technical or batch axis rather than biology. Ask the skill to correlate factors against known covariates (batch, library size, percent-mito) as a sanity check before interpreting any factor as a biological program.

  5. Ground the interpretation and hand off. Ask Claude to name only factors and features that appear in the saved CSVs, citing the variance-explained and loading values. Feed the top-loading gene list of an interesting factor into the functional-enrichment recipe to annotate its biology.

  6. Record provenance. Emit a provenance.json capturing the MOFA+/mofapy2 and muon versions, the number of factors and seed, per-view feature counts, input file sha256s, the run date, and the model id. mofapy2 training is stochastic in initialization; the seed plus the saved mofa_model.hdf5 make the run re-attributable. Keep mofa_run.py, requirements.txt, the CSV tables, and provenance.json under version control.

Why this assembly

Rung 2 of the simplicity ladder. Multi-omics factor analysis is a single statistical procedure with one canonical implementation (mofapy2), and the MOFA+ skill encapsulates the view-construction, scaling, training, and variance-decomposition idioms that are easy to get wrong by hand. Plain Claude Code (rung 1) can write mofapy2 code from docs but reliably fumbles the view-scaling and the per-view feature-selection footgun. A toolbelt (rung 3) buys nothing — the muon/AnnData helper is just the data container feeding the one core skill. Autonomous systems (rung 4) are the wrong tier for a bounded, single-model decomposition.

Availability

Fully open. The MOFA+ and muon skills are community OSS in jaechang-hits/SciAgent-Skills (CC BY 4.0 collection); mofapy2 is LGPL-3.0; muon/MuData is BSD-3-Clause. No subscription, no institutional access, no API key — all computation is local.

Compute requirements

Laptop-sufficient for typical designs. MOFA+ variational inference on a few thousand samples/cells × a few thousand features per view across 3–4 views trains in minutes to tens of minutes on CPU; no GPU. Memory is dominated by the input views held in memory. Very wide views (tens of thousands of features) or hundreds of thousands of cells push training to tens of minutes–hours and 16–32 GB RAM — feature-select aggressively (step 1) rather than throwing raw matrices at the model, and downsample cells per group if needed.

Evidence

Proposed. No documented end-to-end attempt of “Claude + the MOFA+ skill” on a real dataset, with quantitative pass/fail, is known to the curator. The evidence is component-level and strong:

Alternatives considered

  • scVI/scANVI batch integration (rung 2). Integrate multiple single-cell RNA-seq datasets is the answer when the layers are the same modality across batches. Use it for donor/technology harmonization; use MOFA+ for cross-modality factorization. They compose — integrate within modality first, then MOFA+ across modalities.
  • WNN / MultiVI joint embeddings. muon’s Weighted Nearest Neighbors and scvi-tools’ MultiVI produce a joint embedding for clustering paired single-cell modalities but do not give the interpretable per-factor per-view variance decomposition that makes MOFA+ useful for hypothesis generation. Reach for WNN/MultiVI when the deliverable is a joint UMAP; reach for MOFA+ when it is “which axis of variation is shared and what drives it.”
  • Plain Claude Code + mofapy2 (rung 1). Viable for a throwaway one-off if the package is already installed; the skill earns its place by pinning the view-scaling and feature-selection conventions and keeping the run reproducible.
  • An autonomous system (Biomni). Overkill for a single decomposition step; reach for it only when factor analysis is one node in a larger autonomous loop.

See also

Sources

Tried this recipe?

Share feedback — what worked, what didn’t, what you’d change. The form opens with this recipe pre-selected and a link back to this page.