Bio Multi Omics Mofa Integration

Version Compatibility

Reference examples tested with: MOFA2 1.12+ (Bioconductor), mofapy2 0.7+, muon 0.1+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('MOFA2') then ?function_name to verify parameters
• Python: pip show mofapy2 muon then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

MOFA2 is R/Bioconductor and trains through the Python mofapy2 backend via basilisk/reticulate; the trained model serializes to an .hdf5 file whose schema is version-coupled. The R defaults differ from the muon Python defaults (spikeslab_weights, ard_factors, seed 42 vs 1) - confirm per ecosystem.

MOFA2 Integration

"Find shared variation across my omics layers" -> Learn unsupervised latent factors that decompose variation into shared and view-specific axes - because the factors are found WITHOUT the phenotype, so a factor need not separate the groups of interest, and the variance-explained table is the read-out.

• R: create_mofa(data_list) -> prepare_mofa() -> run_mofa()
• Python: mofapy2 for training, muon / mofax for downstream

Scope: unsupervised cross-block factor modeling of bulk omics, variance decomposition, and signed-weight interpretation. Supervised discriminant integration -> mixomics-analysis. The method-selection decision -> integration-design. Per-omic normalization / HVG / transform-to-Gaussian -> data-harmonization. Enrichment mechanics -> pathway-analysis/gsea. Single-cell multimodal MOFA+ -> single-cell/multimodal-integration.

The Single Most Important Modern Insight -- A MOFA Factor Is an Unsupervised Axis of Variance, Not a Pathway, and a Factor That Correlates With Batch Is a Batch Factor

MOFA is PCA generalized to multiple omics: it returns factor scores (Z), per-view loadings (W), and a variance decomposition (R^2 per factor per view) that says which axis is shared across modalities and which is view-specific. That decomposition is the output. Three properties every misuse forgets:

1. Factors are unsupervised and blind to the phenotype. MOFA never saw the groups, so "factor 3 separates cases from controls" is genuinely strong evidence - but only if K factors were not scanned to find the one that splits the groups. A factor becomes biology only after its weights are annotated, it correlates with a known covariate, and it is shown NOT to be a technical-covariate factor.
2. MOFA greedily captures the largest variance. An unregressed batch effect or a depth gradient is often the largest variance, so it becomes a top factor. Always run correlate_factors_with_covariates against batch/depth/plate and exclude technical factors from biological interpretation. A factor that tracks run date is a batch factor, not confounded biology.
3. Factors are unordered and ARD-pruned. Factor 1 is not "the most important" the way PC1 is. The honest deliverable is "factor k explains X% of variance in views {A,B}, its top weights enrich for pathway P at sign s, it correlates with covariate C and not with batch, and it recurs across seeds" - anything less is an axis dressed up as a finding.

Tool Taxonomy (Unsupervised Integrators)

Decision Tree by Deliverable

Prepare the Views

Goal: Get each omic into the orientation and distribution MOFA assumes, so the factors reflect biology rather than measurement scale.

Approach: Each view must be features-by-samples (the transpose of the usual samples-by-features matrix), per-omic normalized and variance-stabilized upstream, and HVG-filtered so feature counts are within an order of magnitude across views. The per-omic transform and HVG selection are owned by data-harmonization.

library(MOFA2)

common <- Reduce(intersect, list(colnames(rna), colnames(prot), colnames(meth)))   # shared samples; mosaic samples may be kept, see below
data_list <- list(RNA=rna[, common], Protein=prot[, common], Methylation=meth[, common])   # each features x samples, already transformed
mofa <- create_mofa(data_list)
plot_data_overview(mofa)        # shows views, samples, and the missing-data pattern (grey = missing, tolerated)

MOFA tolerates missing samples in a view (it ignores missing entries in the likelihood), so a mosaic cohort can be passed directly rather than intersected to complete cases - this is a core reason to choose MOFA when data is incomplete.

Create and Train

Goal: Configure a factor model whose count and likelihoods match the data and the sample size, then train by variational inference.

Approach: Over-specify the factor count and let ARD prune, transform counts to a Gaussian likelihood rather than using Poisson, set a seed for reproducibility, and write the model to a versioned .hdf5.

data_opts  <- get_default_data_options(mofa)       # scale_views=FALSE, center_groups=TRUE
model_opts <- get_default_model_options(mofa)      # num_factors=10, likelihoods='gaussian'
train_opts <- get_default_training_options(mofa)   # convergence_mode='fast', drop_factor_threshold=-1, stochastic=FALSE, seed=42

model_opts$num_factors <- 15                       # over-specify; ARD prunes inactive factors
model_opts$likelihoods <- c(RNA='gaussian', Protein='gaussian', Methylation='gaussian')   # transform counts upstream, prefer gaussian
train_opts$drop_factor_threshold <- 0.01           # drop factors explaining <1% variance in ALL views
data_opts$scale_views <- TRUE                      # equalize per-view variance if feature counts cannot be balanced by filtering

mofa <- prepare_mofa(mofa, data_options=data_opts, model_options=model_opts, training_options=train_opts)
mofa <- run_mofa(mofa, outfile=file.path(tempdir(), 'model.hdf5'), use_basilisk=TRUE)

Read the Variance Decomposition (the Output)

Goal: Identify which factors are shared across views and which are view-specific before interpreting any of them.

Approach: The R^2 per factor per view is the central result: a factor active in two or more views is a shared axis, a factor active in one view is view-specific. Inspect this table first; factors with near-zero R^2 everywhere are noise.

var_exp <- get_variance_explained(mofa)            # $r2_total, $r2_per_factor  <- the central output
plot_variance_explained(mofa)                      # heatmap: factors x views
plot_variance_explained(mofa, plot_total=TRUE)     # total variance explained per view

Label the Factors (and Exclude Technical Ones)

Goal: Earn a biological label for a factor instead of asserting one from its existence.

Approach: Attach metadata after fitting (it is never used to train), correlate each factor with both biological and technical covariates, and exclude any factor that tracks batch/depth/plate from biological interpretation. Then run enrichment on the signed weights, treating the two poles of the axis separately.

md <- metadata[unlist(samples_names(mofa)), ]
md$sample <- rownames(md)                       # samples_metadata<- requires a literal 'sample' column
samples_metadata(mofa) <- md
correlate_factors_with_covariates(mofa, covariates=c('condition', 'batch', 'depth'))   # a factor that correlates with batch IS a batch factor

# enrichment per sign - the two poles are biological opposites along one axis
up   <- run_enrichment(mofa, view='RNA', feature.sets=msig_binary_matrix, factors=1:5, sign='positive')
down <- run_enrichment(mofa, view='RNA', feature.sets=msig_binary_matrix, factors=1:5, sign='negative')

Multi-group MOFA (group in create_mofa) partitions samples so factor activity can differ across groups while weights stay shared - it asks "do the same axes operate within each group?" It is NOT batch correction and putting the phenotype in as a group does not make MOFA supervised. For samples with a spatial or temporal coordinate, MEFISTO (a GP-prior mode of MOFA2, via mefisto_options + set_covariates) learns factors that vary smoothly along that covariate.

Per-Method Failure Modes

Factor narrated as a mechanism

Trigger: "factor 1 represents immune activation" from the factor's existence. Mechanism: a factor is a direction of covariation the ARD prior kept; it has no intrinsic meaning. Symptom: a biological story with no enrichment, no covariate correlation, no replication. Fix: report the R^2 footprint, the signed-weight enrichment, and the known-covariate correlation; call it hypothesis-generating until validated.

Batch factor mistaken for biology

Trigger: interpreting a top factor without a technical-covariate check. Mechanism: MOFA captures the largest variance, and unregressed batch is often largest. Symptom: the top factor tracks run date / plate better than phenotype. Fix: regress known batch out upstream (before HVG selection); always correlate_factors_with_covariates and exclude technical factors.