Getting started with MOFA-FLEX#

MOFA-FLEX is a framework for factor analysis of multimodal data, with a focus on single-cell omics, representing an observed data matrix \(\mat{Y}\) as a sum matrices \(\mat{X}_i\), each of which is a product of a low-rank factor matrix \(\mat{Z}_i\) and a low-rank weight matrix \(\mat{W}_i\): \(\mat{Y} \approx \sum_{i=1}^T\mat{Z}_i \mat{W}_i\). In the simplest and most common case, there is only one additive term, i.e. \(T = 1\). In this case, MOFA-FLEX is essentially a synthesis of MOFA[AAB+20, AVA+18], MEFISTO[VBA+22], MuVI[QB23], nonnegative matrix factorization, and NSF[TE23]. It natively supports MuData and AnnData objects, integrating into the wider scverse ecosystem.

For this notebook, we will use the pbmc3k dataset by 10x Genomics, conveniently packaged in the mudatasets package. The dataset contains 10 000 single cells profiled with the 10x multiome assay, capturing gene expression (scRNAseq) and chromatin accessibility (scATACseq).

import mudata as md
import mudatasets as mds
import muon as mu
import mofaflex as mfl
import scanpy as sc
from plotnine import *

Importing the dtw module. When using in academic works please cite:
  T. Giorgino. Computing and Visualizing Dynamic Time Warping Alignments in R: The dtw Package.
  J. Stat. Soft., doi:10.18637/jss.v031.i07.

theme_set(theme_bw())
md.set_options(display_style="html", display_html_expand=0);

import warnings

warnings.filterwarnings("ignore", category=FutureWarning)

mdata = mds.load("pbmc3k_multiome")

■ File filtered_feature_bc_matrix.h5 from pbmc3k_multiome has been found at /home/kats/mudatasets/pbmc3k_multiome/filtered_feature_bc_matrix.h5
■ Checksum is validated (md5) for filtered_feature_bc_matrix.h5
■ Loading filtered_feature_bc_matrix.h5...
Added `interval` annotation for features from /home/kats/mudatasets/pbmc3k_multiome/filtered_feature_bc_matrix.h5

mdata.var_names_make_unique()
mdata

MuData object 2711 obs × 134920 var in 2 modalities

Metadata

.obs0 elements

No metadata

Embeddings & mappings

.obsm2 elements

rna	bool	numpy.ndarray	1 columns
atac	bool	numpy.ndarray	1 columns

Distances

.obsp0 elements

No distances

rna

2711 × 36601

AnnData object 2711 obs × 36601 var

Matrix

float32 scipy.sparse._csr.csr_matrix

Layers

.layers0 elements

No layers

Metadata

.obs0 elements

No metadata

Embeddings

.obsm0 elements

No embeddings

Distances

.obsp0 elements

No distances

Miscellaneous

.uns0 elements

No miscellaneous

atac

2711 × 98319

AnnData object 2711 obs × 98319 var

Matrix

float32 scipy.sparse._csr.csr_matrix

Layers

.layers0 elements

No layers

Metadata

.obs0 elements

No metadata

Embeddings

.obsm0 elements

No embeddings

Distances

.obsp0 elements

No distances

Miscellaneous

.uns0 elements

No miscellaneous

Preprocessing#

We first perform the basic preprocessing steps outlined in the muon tutorial to remove undetected genes and poor-quality cells.

rna = mdata["rna"]
rna.var["mt"] = rna.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    rna, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
mu.pp.filter_var(rna, "n_cells_by_counts", lambda x: x >= 3)
mu.pp.filter_obs(rna, "n_genes_by_counts", lambda x: (x >= 200) & (x < 5000))

mu.pp.filter_obs(rna, "total_counts", lambda x: x < 15000)
mu.pp.filter_obs(rna, "pct_counts_mt", lambda x: x < 20)

atac = mdata.mod["atac"]
sc.pp.calculate_qc_metrics(atac, percent_top=None, log1p=False, inplace=True)
mu.pp.filter_var(atac, "n_cells_by_counts", lambda x: x >= 10)
mu.pp.filter_obs(atac, "n_genes_by_counts", lambda x: (x >= 2000) & (x <= 15000))
mu.pp.filter_obs(atac, "total_counts", lambda x: (x >= 4000) & (x <= 40000))

We further normalize the data and determine highly variable genes. MOFA-FLEX will automatically use only the highly variable genes for the analysis.

sc.pp.normalize_total(rna, target_sum=1e4)
sc.pp.log1p(rna)
sc.pp.highly_variable_genes(rna, min_mean=0.02, max_mean=4, min_disp=0.5)

mu.atac.pp.tfidf(atac, scale_factor=1e4)
sc.pp.normalize_per_cell(atac, counts_per_cell_after=1e4)
sc.pp.log1p(atac)
sc.pp.highly_variable_genes(atac, min_mean=0.05, max_mean=1.5, min_disp=0.5)

mdata.update()

mdata

MuData object 2695 obs × 119335 var in 2 modalities

Metadata

.obs7 elements

rna:n_genes_by_counts	float64	2272.00,3254.00,1798.00,1145.00,1495.00,1342.00,25...
rna:total_counts	float32	4747.00,7761.00,3666.00,2162.00,2910.00,2616.00,66...
rna:total_counts_mt	float32	369.00,693.00,409.00,271.00,293.00,213.00,308.00,6...
rna:pct_counts_mt	float32	7.77,8.93,11.16,12.53,10.07,8.14,4.62,12.36,10.79,...
atac:n_genes_by_counts	float64	5415.00,9572.00,9459.00,6634.00,7107.00,7378.00,56...
atac:total_counts	float32	12757.00,23124.00,22912.00,17106.00,18400.00,18723...
atac:n_counts	float32	6129.60,10004.25,9936.28,6779.99,6849.25,7125.32,5...

Embeddings & mappings

.obsm2 elements

rna	bool	numpy.ndarray	1 columns
atac	bool	numpy.ndarray	1 columns

Distances

.obsp0 elements

No distances

rna

2636 × 21256

AnnData object 2636 obs × 21256 var

Matrix

float32 scipy.sparse._csr.csr_matrix

Layers

.layers0 elements

No layers

Metadata

.obs4 elements

n_genes_by_counts	int32	2272,3254,1798,1145,1495,1342,2589,2385,3009,3320,...
total_counts	float32	4747.00,7761.00,3666.00,2162.00,2910.00,2616.00,66...
total_counts_mt	float32	369.00,693.00,409.00,271.00,293.00,213.00,308.00,6...
pct_counts_mt	float32	7.77,8.93,11.16,12.53,10.07,8.14,4.62,12.36,10.79,...

Embeddings

.obsm0 elements

No embeddings

Distances

.obsp0 elements

No distances

Miscellaneous

.uns2 elements

log1p	dict	1 element	base: None
hvg	dict	1 element	flavor: seurat

atac

2450 × 98079

AnnData object 2450 obs × 98079 var

Matrix

float32 scipy.sparse._csr.csr_matrix

Layers

.layers0 elements

No layers

Metadata

.obs3 elements

n_genes_by_counts	int32	5415,9572,9459,6634,7107,7378,5654,5212,13768,1387...
total_counts	float32	12757.00,23124.00,22912.00,17106.00,18400.00,18723...
n_counts	float32	6129.60,10004.25,9936.28,6779.99,6849.25,7125.32,5...

Embeddings

.obsm0 elements

No embeddings

Distances

.obsp0 elements

No distances

Miscellaneous

.uns2 elements

log1p	dict	1 element	base: None
hvg	dict	1 element	flavor: seurat

Fitting a model#

A MOFA-FLEX model is assembled from one or multiple additive terms, each of which takes parameters influencing the model. A basic understanding of the model is required to correctly set the parameters for a given dataset.

We will start with the simplest case of only one additive term. There is currently one one type of additive term implemented: The MofaFlex term, which is a product of two low-rank matrices \(\mat{Z}\mat{W}\), as described above. Additional term types may be added in the future.

The MofaFlex term type uses prior distributions on the factors and weights. Available priors are listed in the API documentation and can be passed to the factor_prior and weight_prior arguments of MofaFlex. Priors that do not require any mandatory arguments can be passed both as string, e.g. weight_prior="Normal", and as an instance of the respective class, e.g. weight_prior=mfl.priors.Normal(). Some priors, like the InformedHorseshoe, require mandatory arguments. These priors can only be passed as an instance, e.g. weight_prior=mfl.priors.InformedHorseshoe(annotation_varm_key="some_key").

The examples above specified the same prior for all views. Different priors for different views can be specified by passing a dictionary to the corresponding argument. For example,weight_prior={"view_1": "Normal", "view_2": "SpikeSlab", ("view_3", "view_4"): "Horseshoe"} uses the Normal prior for view view_1, the spike and slab prior for view view_2, and the Horseshoe prior for views view_3 and view_4. The same applies to factor priors and groups.

For this introductory tutorial, we will use the default Normal prior on the factors and the spike and slab prior on the weights

model = mfl.terms.MofaFlex(n_factors=15, weight_prior="SpikeSlab")

This model corresponds to a standard MOFA model. It is an instance of the MOFAFLEX class.

Since we have normalized the data in the previous section, we will use a Normal (Gaussian) likelihood. In general, a negative binomial likelihood is more appropriate for unnormalized count data. Similar to how the priors can be specified either as strings or instances of the respective classes, likelihoods for all or individual views can be given as strings or instances of likelihood classes to fit. If a likelihood class has required arguments, it can not be used with the string form.

With default settings, MOFA-FLEX will plot the number of observations in each view and group. This can be turned off by setting the plot_data_overview argument to fit to False. MOFA-FLEX will also automatically save the final model after training which can later be loaded with MOFAFLEX.load. The name of the saved file can be changed using the save_path argument to fit.

model.fit(mdata, likelihoods="Normal", batch_size=1000, seed=42)

../_images/6e9551a9bac397f7c8b9bd9a2d11703571d49e8989a7339cdff606c789663f8a.png

We can plot the loss curve to get an overview of the training process:

mfl.pl.training_curve(model)

../_images/43a6a1876e61fa059e68a44d62bc2c3943979243510e75169c6d2e2e6db1048e.png

Analysing the model#

The trained model can now be analysed to assess its quality and discover relevant biology. Most relevant information is contained in the addtive terms themselves, as can be seen from the API documentation. Because having only one term in the model is so common, this case is handled specially: The MOFAFLEX object forwards requests for any attributes that it doesn’t know about to its only term. We can therefore largely ignore anything related to additive terms and simply pass the MOFAFLEX object to any downstream analysis function.

First, we plot the correlations between factors, which should be as low as possible: We want factors to be uncorrelated, as each factor should capture a different aspect of the data.

mfl.pl.factor_correlation(model)

../_images/6fcd1541e79a1ac4bd5932897d80415b2f89afaf9800328e9699f00747f1efad.png

This is not great, but good enough for this introduction. We can now plot the fraction of variance in the data explained by each factor to determine the most important factors. Note that for non-Normal likelihoods, these values are approximate.

mfl.pl.variance_explained(model)

../_images/7bb169c72328588efe28b98c18ed10ab705392bfc44ac5b4a0a1766cf5730b41.png

We can plot the most important genes per factor to get an idea of what each factor represents. This plot aggregates over all views.

mfl.pl.top_weights(model, figsize=(10, 10))

../_images/725a1b647096f5ed2f3d9affe7c8faafd621fba69225e7c09c63db480a2b2951.png

It appears that the most important features are all from the RNA dataset. We can use a different plotting function to get the top weights individually for each view.

mfl.pl.weights(model, factors=(1, 2), figsize=(10, 10))

../_images/309b195347b0deacef6de40e15f62588a0e738042489f212b7cb273835549206.png

We can also plot factors agains each other. This may be useful to define clusters of cells with similar factor values.

mfl.pl.factors_scatter(model, 1, 2, alpha=0.5)

../_images/2067188aa4bc4c520da204234591b002492e200a78566e420e1d331680318cbf.png

Some priors, like the spike and slab prior we used here, provide additional, prior-specific functionality. It is documented with the respective prior class, but it is only available through the model or term objects. For example, we can retrieve the estimated probabilites that a weight is non-zero from the spike and slab prior:

model.get_sparse_weight_probabilities()["rna"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5	Factor 6	Factor 7	Factor 8	Factor 9	Factor 10	Factor 11	Factor 12	Factor 13	Factor 14	Factor 15
ISG15	0.907351	0.967472	0.930406	0.989911	0.996477	0.987474	0.971301	0.987840	0.997191	0.989061	0.869097	0.992749	0.908470	0.996628	0.997284
C1orf159	0.020836	0.021660	0.035837	0.030507	0.013763	0.016295	0.041932	0.032445	0.337761	0.024845	0.047659	0.150292	0.034279	0.016074	0.162518
AL390719.3	0.014695	0.008472	0.010901	0.011860	0.006212	0.012545	0.012280	0.008580	0.012701	0.009219	0.007854	0.007598	0.012279	0.009544	0.008046
TNFRSF18	0.018560	0.968644	0.020908	0.996613	0.012198	0.062231	0.022805	0.027726	0.059821	0.015497	0.982226	0.977191	0.039064	0.015577	0.012805
TNFRSF4	0.022974	0.980925	0.041207	0.995690	0.018080	0.015872	0.030083	0.015865	0.035173	0.923943	0.996998	0.990913	0.983734	0.024154	0.065391
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
PDZD4	0.019661	0.985020	0.013107	0.988327	0.022756	0.021707	0.007589	0.028342	0.011482	0.013092	0.013994	0.009747	0.015219	0.012465	0.012936
SLC10A3	0.014703	0.009161	0.010335	0.023204	0.020582	0.014217	0.009556	0.010545	0.009974	0.010174	0.009625	0.009961	0.064897	0.011346	0.012347
MTCP1	0.015480	0.017789	0.011080	0.024393	0.013375	0.013496	0.020470	0.016989	0.015864	0.017766	0.019190	0.021110	0.019990	0.025528	0.012766
BRCC3	0.832569	0.036314	0.044227	0.035232	0.019139	0.021913	0.023558	0.021601	0.030561	0.020681	0.064110	0.029122	0.016254	0.021355	0.022014
VBP1	0.034023	0.937453	0.027866	0.099840	0.075255	0.029844	0.935061	0.115212	0.026995	0.037032	0.015042	0.022530	0.045503	0.022754	0.989106

3718 rows × 15 columns

Of course, MOFA-FLEX does not and cannot provide all imaginable analysis functions. It thus provides methods to access the factor and weight values, such that they can be used for manual analysis.

weights = model.get_weights()
factors = model.get_factors()

weights["rna"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5	Factor 6	Factor 7	Factor 8	Factor 9	Factor 10	Factor 11	Factor 12	Factor 13	Factor 14	Factor 15
ISG15	-0.307310	0.147454	-0.168437	-0.548156	0.550673	-0.615919	-0.209126	0.759273	0.692284	0.353094	-0.180655	0.501951	-0.212744	0.691283	-1.386797
C1orf159	-0.000466	0.000449	-0.001714	0.001626	0.000178	0.000157	-0.001347	0.000987	-0.032292	0.000348	0.003784	-0.010668	0.001494	0.000258	0.015819
AL390719.3	-0.000667	0.000051	-0.000198	0.000128	0.000260	0.000356	0.000367	0.000019	-0.000090	-0.000123	0.000030	-0.000028	0.000026	0.000068	0.000026
TNFRSF18	-0.000911	0.157169	-0.000240	-0.537049	-0.000136	0.002771	0.000333	-0.000901	-0.006406	0.000278	-0.215397	-0.229190	-0.003172	-0.000282	0.000249
TNFRSF4	-0.000034	0.299320	-0.001279	-0.515594	-0.000279	0.000474	0.000132	-0.000084	-0.001507	-0.096456	-0.532599	-0.276027	0.099232	-0.001121	-0.002782
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
PDZD4	-0.000487	-0.158583	0.000009	-0.147062	0.000727	-0.000305	0.000097	-0.001088	0.000008	-0.000218	-0.000235	0.000216	0.000193	0.000226	-0.000087
SLC10A3	0.000132	-0.000037	0.000166	-0.001114	-0.000795	0.000501	0.000170	-0.000342	0.000184	-0.000159	-0.000115	0.000168	-0.002455	0.000172	-0.000367
MTCP1	-0.000035	-0.000052	-0.000117	-0.000556	-0.000222	0.000334	0.000722	0.000105	-0.000074	-0.000273	0.000549	-0.000361	0.000083	-0.000626	-0.000029
BRCC3	0.077248	-0.001735	0.001343	-0.001249	-0.000754	-0.000438	-0.000409	0.001033	0.001251	-0.000324	0.002884	-0.000354	-0.000051	-0.000116	0.000231
VBP1	-0.000056	-0.118339	-0.000335	-0.004080	0.006114	0.000078	0.122198	-0.004814	-0.000183	0.000498	-0.000124	-0.000868	-0.002151	0.000853	-0.211334

3718 rows × 15 columns

weights["atac"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5	Factor 6	Factor 7	Factor 8	Factor 9	Factor 10	Factor 11	Factor 12	Factor 13	Factor 14	Factor 15
chr1:267561-268455	0.001055	0.000107	0.000382	0.000647	0.000507	-0.000404	-0.000364	0.003605	-0.000315	-0.000055	0.000534	0.000200	-0.000228	0.000254	0.000564
chr1:629484-630393	-0.029646	0.001389	0.000089	0.000714	-0.173619	0.016173	0.110926	0.000113	-0.043532	0.002735	-0.003180	0.002305	0.000511	-0.001797	0.000167
chr1:778284-779202	-0.416556	-0.002964	-0.000145	-0.000290	0.000106	0.039963	0.019213	-0.000710	0.001283	-0.000204	-0.093867	-0.000176	-0.000167	-0.042418	-0.000236
chr1:844149-845034	-0.000835	0.000119	-0.459350	0.058479	-0.000066	0.006276	0.001320	0.000291	-0.004066	0.002193	-0.001564	0.000009	0.000193	-0.001403	-0.000167
chr1:857873-858632	-0.075457	-0.000059	-0.129848	0.130800	-0.001041	0.012445	0.001186	-0.000691	-0.010941	0.002422	0.000064	-0.000653	0.000755	-0.001552	-0.000193
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
GL000219.1:125017-125889	-0.005248	-0.004037	-0.000345	-0.000036	-0.001405	0.000428	0.002936	-0.000609	-0.000227	0.000182	-0.000779	0.000203	-0.000119	-0.001376	-0.000422
KI270721.1:2089-2980	-0.000860	-0.000382	0.000067	-0.000019	-0.000029	0.000190	0.000234	-0.000398	-0.002710	0.000156	-0.000878	0.000301	0.000319	-0.001503	0.000112
KI270726.1:27153-28037	-0.000689	0.000280	-0.004013	0.000663	-0.000121	0.003187	0.167951	-0.000183	-0.001330	0.000599	-0.000675	-0.000075	0.007371	-0.001064	-0.000152
KI270726.1:41489-42329	-0.001672	-0.000139	-0.002346	0.000202	-0.000907	0.076016	0.087928	-0.001056	-0.005532	0.000444	-0.000063	-0.000532	0.000161	-0.005621	-0.000135
KI270713.1:29578-30400	-0.010091	-0.000255	-0.132934	0.000914	-0.000212	0.001437	0.000864	-0.000690	-0.011513	0.000631	-0.001101	0.000048	0.000121	-0.001666	-0.000606

20985 rows × 15 columns

factors["group_1"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5	Factor 6	Factor 7	Factor 8	Factor 9	Factor 10	Factor 11	Factor 12	Factor 13	Factor 14	Factor 15
AAACAGCCAAATATCC-1	-0.554957	-0.636028	-0.016920	-1.095950	0.261986	0.371856	-0.242503	-0.301752	0.080081	0.032883	0.202017	-0.029437	-0.875365	0.211991	-0.063476
AAACAGCCAGGAACTG-1	0.612035	0.286613	0.207965	0.055287	0.239285	-0.215778	-0.414638	0.159871	0.277654	0.093675	0.227588	0.188399	-0.056610	0.261375	-0.297113
AAACAGCCAGGCTTCG-1	0.354993	0.127425	0.312819	0.082485	0.288181	-0.405416	-0.373531	0.360576	0.078753	-0.483631	0.211231	0.042947	-0.014426	0.179994	0.821503
AAACCAACACCTGCTC-1	0.032957	0.366913	0.427888	-0.051962	-0.880440	0.334719	0.197149	-0.248163	-0.101071	0.163213	0.011708	0.699828	-0.173349	-0.078428	0.242243
AAACCAACAGATTCAT-1	-0.247417	-1.135546	0.035206	-0.585370	0.270555	0.243265	0.073886	-0.238663	0.012999	-0.032337	0.135770	0.174433	-0.261372	0.190698	-0.162544
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
TGTTGTGCATTGTGTG-1	-0.142208	-0.039271	-0.363117	0.392965	0.006957	0.309290	0.337234	0.133889	0.069136	0.069515	0.023291	-0.076985	-0.062623	0.021563	0.079993
TTCCCACAGGCGAATA-1	0.341569	0.231984	0.340370	0.033196	0.199103	-0.229104	-0.359265	0.217055	0.106105	-0.053146	0.219174	-0.062496	-0.057537	0.245146	-0.057080
TTGCTTAGTCCTTTAA-1	0.145313	0.185044	0.316008	-0.024275	0.248887	0.097686	-0.204765	0.193208	0.095867	-0.170530	0.151687	-0.005056	-0.106132	0.142125	0.115739
TTGTCCGGTGCATTAG-1	0.315805	0.180925	0.285874	0.017694	0.162500	-0.222266	-0.414497	0.239941	0.049338	-0.081717	0.254962	0.000952	-0.045032	0.169013	0.122050
TTTGTTGGTAGGTTTG-1	0.133116	0.259846	0.350169	-0.131188	-0.717402	0.056519	0.042166	-0.121947	-0.270413	0.068795	0.134479	0.561410	-0.050457	-0.131921	-0.052049

2695 rows × 15 columns

Working with multiple additive terms#

The general workflow for multiple terms is very similar to the above. We first assemble a model object from the terms. Note that when working with multiple terms, each term must get a unique name as its first argument. These names will be used later to identify each term during downstream analysis.

model = mfl.terms.MofaFlex("spikeslab", n_factors=5, weight_prior=mfl.priors.SpikeSlab()) + mfl.terms.MofaFlex("hs", n_factors=5, weight_prior="Horseshoe")

Using several terms with the same name (including the default, which is used if no name is specified) will not work:

model + mfl.terms.MofaFlex("hs", n_factors=3)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[22], line 1
----> 1 model + mfl.terms.MofaFlex("hs", n_factors=3)

File /data/ilia/famo/famo/src/mofaflex/_core/mofaflex.py:95, in MOFAFLEX.__add__(self, other)
     93     raise ValueError("Cannot add terms to an already trained model.")
     94 if len(intersection := self._terms.keys() & other._terms.keys()) > 0:
---> 95     raise ValueError(
     96         f"Operands must have unique term names, but terms {', '.join(intersection)} were found in both operands."
     97     )
     98 return __class__(**self._terms, **other._terms)

ValueError: Operands must have unique term names, but terms hs were found in both operands.

The model can be trained as before.

model.fit(mdata, likelihoods="Normal", plot_data_overview=False, batch_size=1000, seed=42)

We can still access global properties of the model through the MOFAFLEX object:

mfl.pl.training_curve(model)

../_images/a572e3c22f209a53c557783008f2871002962aed34acb4788e597717bf46bd29.png

However, each term now has its own set of factors and weights. Terms can be accessed through model.terms. For example, to plot factor correlations, individual terms need to be passed to the plotting function.

mfl.pl.factor_correlation(model.terms["spikeslab"])

../_images/8cf1ae4aac8043463a4dae850be5d4572bdd449a403f6a012a3b48d8c2d95eec.png

mfl.pl.factor_correlation(model.terms["hs"])

../_images/28e3e451fa61039906c4ad3b3a4a56323fcaf0b97bbefab2ead65ed05453ef14.png

We can retrieve factors and weights from each term individually:

model.terms["spikeslab"].get_weights()["rna"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5
ISG15	-0.049461	-0.237090	-0.728369	0.024563	1.160542
C1orf159	-0.014093	0.000382	-0.001730	0.000127	-0.004358
AL390719.3	0.000019	0.000224	0.000252	-0.000075	0.000195
TNFRSF18	0.000858	0.484985	0.000087	0.000328	0.000513
TNFRSF4	0.072538	0.707179	-0.000175	0.000211	0.000584
...	...	...	...	...	...
PDZD4	-0.000029	-0.000045	-0.000055	-0.000123	0.000187
SLC10A3	0.000647	-0.000532	0.000093	0.000324	0.000882
MTCP1	-0.000215	0.000174	-0.000006	0.000055	-0.000247
BRCC3	0.000124	0.000354	-0.000259	0.000462	-0.000260
VBP1	0.001148	0.000934	0.000967	0.001477	0.193498

3718 rows × 5 columns

model.terms["hs"].get_weights()["atac"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5
chr1:267561-268455	0.005559	-0.020516	0.000633	-0.024448	-0.010069
chr1:629484-630393	0.028555	0.011072	0.077857	-0.030109	0.029489
chr1:778284-779202	-0.000888	0.040224	0.012272	0.159778	0.042900
chr1:844149-845034	0.206968	0.129121	0.166077	0.032500	0.044504
chr1:857873-858632	0.111969	0.118554	0.110831	-0.000843	0.002059
...	...	...	...	...	...
GL000219.1:125017-125889	0.012977	0.020929	0.000248	0.021503	0.025031
KI270721.1:2089-2980	0.005144	-0.001238	0.014928	-0.006683	0.027719
KI270726.1:27153-28037	0.033519	0.091075	0.053493	0.007624	0.042490
KI270726.1:41489-42329	0.071078	0.061600	0.012454	-0.011959	0.030795
KI270713.1:29578-30400	0.091335	0.034532	0.040051	-0.010140	0.037878

20985 rows × 5 columns

Similarly, any prior-specific methods and properties can be accessed through the term objects:

model.terms["spikeslab"].get_sparse_weight_probabilities()["atac"]

	Factor 1	Factor 2	Factor 3	Factor 4	Factor 5
chr1:267561-268455	0.018063	0.014639	0.019392	0.019494	0.022256
chr1:629484-630393	0.025180	0.013192	0.105004	0.969513	0.020776
chr1:778284-779202	0.015266	0.024299	0.990157	0.012150	0.014349
chr1:844149-845034	0.043671	0.074466	0.040380	0.018401	0.013577
chr1:857873-858632	0.015195	0.018924	0.012104	0.014554	0.015063
...	...	...	...	...	...
GL000219.1:125017-125889	0.022600	0.012115	0.013100	0.018671	0.011112
KI270721.1:2089-2980	0.012974	0.021038	0.016087	0.012599	0.013158
KI270726.1:27153-28037	0.024110	0.025697	0.021489	0.024748	0.016255
KI270726.1:41489-42329	0.016723	0.021952	0.513801	0.013328	0.015378
KI270713.1:29578-30400	0.040418	0.020021	0.014511	0.014389	0.016820

20985 rows × 5 columns

Adding terms to an already trained model will not work:

model + mfl.terms.MofaFlex(n_factors=6)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[30], line 1
----> 1 model + mfl.terms.MofaFlex(n_factors=6)

File /data/ilia/famo/famo/src/mofaflex/_core/mofaflex.py:93, in MOFAFLEX.__add__(self, other)
     91     return NotImplemented
     92 if hasattr(self, "_model") or hasattr(other, "_model"):
---> 93     raise ValueError("Cannot add terms to an already trained model.")
     94 if len(intersection := self._terms.keys() & other._terms.keys()) > 0:
     95     raise ValueError(
     96         f"Operands must have unique term names, but terms {', '.join(intersection)} were found in both operands."
     97     )

ValueError: Cannot add terms to an already trained model.

Getting started with MOFA-FLEX

Contents

Getting started with MOFA-FLEX#

Preprocessing#

Fitting a model#

Analysing the model#

Working with multiple additive terms#