$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
Following the successful execution of the workflow, several tables and figures are generated as indicated by Figure 2. Figures are placed in the /figures folder (Figure 6, Figure 7, Figure 8, Supplementary Figure 1, Supplementary Figure 2, Supplementary Figure 3, Supplementary Figure 4), and tables will be placed in the specified /results folder.
In case the workflow execution is not successful, this might be mainly due to: technical errors caused, for example, by insufficient memory (especially in the first step where a large single-cell dataset is loaded), incorrectly formatted data (e.g., non-matching sample_id columns between datasets) or incorrect specifications in the configuration files (e.g., excluding to many features). In this case, usually, an error message within the Jupyter-notebook script will occur during execution and no plots and data will be generated. It is recommended to use the default configuration files as generated during the script execution and only modify specific parameters as described in the protocol.
A successful execution is indicated by the generation of the resulting plots and tables, and each step will reveal additional information about the data and the main patterns of variance inherent in it. Still, not necessarily each execution will produce biologically useful and interpretable results. Often, the data is characterized by large technical effects and different distributions, which need to be accounted for in the 'Data Pre-Processing and Harmonization' step or in the 'MOFA9 Model' (which also allows specifying different distributions for the input data types) to be able to extract the variation of the data that reflects the underlying biological processes.
Within the presented workflow, different multi-omic datasets may be used as input. Currently, the workflow accepts the popular .h5ad file format for single-cell data and a very general .csv file format for all the other datasets as input (Figure 3). It is common for different omic datasets to have very different file formats. To not limit the execution of the workflow to specific file formats, .csv is used as a very general format. Therefore, all kinds of different omics datasets can be used as input for the workflow but need to be converted to the corresponding .csv format as indicated in Figure 3 first before usage within this workflow. This might be prepared either using a spreadsheet or omic specific software. To preprocess the different omics datasets, several options are available within the workflow to apply different pre-processing and normalization steps (e.g., library size adjustment, log transformation, sample quantile normalization) on the different input datasets by configuring the 02_Pre_Processing_Configs.csv and 02_Pre_Processing_Configs_SC.csv file (Figure 2). Nevertheless, the available options here are mainly based on the specific input data available in the data set presented here (scRNA-seq, cytokine assay, proteomics, prime-seq). In case other omics/data types are used, it might be necessary to apply additional omic-specific normalization steps according to existing best practices. In this case the data can be handed over to the workflow in an already pre-processed form and will be integrated along with the other datasets without applying further pre-processing steps. In many cases, applying the Feature Wise Quantile Normalization step is useful in aligning the distribution of all data types to a normal distribution and making downstream analysis between the different input features more comparable and compatible with the model specification of Gaussian noise.
During the workflow execution several plots and outputs are generated that support the process of data integration and subsequent biological downstream interpretation. For scRNA-seq data, the plot in FIG01_Amount_of_Cells_Overview (Figure 6) indicates which cell types might include too few cells per sample and cell type to reliably measure a gene expression signal as for subsequent analyses, the mean value across all cells of a cell type per sample is used as an expression estimate (psedobulk- approach). In this use-case, we exclude cell types that have less than three cells in the majority of samples.
The variance decomposition plot FIG03_Overview_Variance_Decomposition (Figure 7, Supplementary Figure 1) can indicate how well the different data sources integrate and how much of the variance in the different data sources is shared and unique to each data source. For example, testing different pre-processing strategies on the dataset used here shows for example that removing the Feature Wise Quantile normalization step from the pre-processing leads to latent factors that are more focused on specific data views and reduces the integration of the proteomic data with the other data sources. This can be seen in the reduced amount of explained variance (Supplementary Figure 1B). Running the MOFA model without any filtering of features or without normalization leads to less shared variance between the different views captured by the latent factors (Supplementary Figure 1C). This indicates that latent factors predominantly reflect data type specific technical effects. Besides that, the MOFA9 model itself might also return warnings in case of poorly pre-processed data. An example of such a warning is shown in Supplementary Figure 1 for the alternative pre-processing configurations MI_v2 and MI_v3 (the specific example configuration files are stored in the cloned GitHub repository in the config_examples folder).
Additionally, after running the MOFA model, the results can be evaluated in several downstream analyses by associating the factor with known biological meta-information about the samples, as well as technical and other confounding covariates (04_Downstream_Factor_Analysis) to identify the likely cause for the variation captured by the factors. For example, if one of the factors of the MOFA models associates strongly with one of the technical covariates (such as batch information) this might indicate that this factor rather captures technical variation within the data instead of biological variation.
To narrow down the biological interpretation in the downstream analysis part, a couple of findings based on the input dataset (a more refined interpretation can be found in the original publication11) are outlined here. In the first step, we could observe that with the applied pre-processing strategy, we find several factors that capture variance across multiple cell types but also other omics data types (Figure 7A). For example, Factor 2 captures variance in the clinical input features and in several cell types of the scRNA-seq dataset. Associating the first three factors with relevant clinical covariates like 'CRP' and 'CK' (Figure 7B) and investigating the differences in factor values for the different patient subgroups: 'Control (including CCS and non-CCS) vs. 'ACS' measured at the different time points (TP1-TP4) (Figure 7C), we also find that Factor2 associates significantly with the 'CK' value and Factor3 with the 'CRP' value. At the same time, 'ACS' samples at TP1 and TP2 (which reflect the acute phase of the immune response to myocardial infarction (MI)) show an increase in factor values compared to 'Control' and later timepoint samples (TP3/TP4). CK is a known marker of myocardial damage and is typically characterized by increased values at TP1/TP2, similar to the pattern captured by Factor2.
To generate insights into the biological processes shaping Factor2 we evaluate the top-ranking features of the factor by looking at the feature weights table generated by the model (03_Weight_Data.csv). Analyzing the top 1% of features with the highest absolute weights on the factor, we find mainly CD4.TCM and CD14. Mono-derived features are overrepresented compared to their overall number of input features (Figure 8A), indicating these cell types to be highly relevant in the inflammatory process after MI (NOTE: in case no feature-wise quantile normalization was applied in the pre-processing different distributions of the features might also affect this result and evaluation should be done separately by data type). Analyzing the top-ranking features of the CD4.TCM cell type on the factor, we find several interesting genes like EIF3E18 required for robust T-cell activation and HMGB119, which promotes expansion and activation of T-cells (Figure 8B). Next, we run the pathway enrichment analysis using immune pathways from the REACTOME20 database as a pathway set (Prepared_Pathway_Data.csv). We find enrichment for several 'Interleukin' pathways, including 'Interleukin-6' signaling. The expression levels of several genes in different cell types of the scRNA-seq data and the 'IL6' cytokine values measured by the Cytokine assay contributed to this result (Figure 8C). Identification of these shared patterns across data types highlights the added value of an integrated analysis. Overall, this approach can also identify several other factors that reflect disease state or associate treatment outcome and the underlying multicellular immune programs as is described more in detail in the corresponding publication11.
To further stress the advantage of integrated analyses across multiple omics, the same workflow was also run only including the proteomics input data (Supplementary Figure 4). Analyzing the resulting factors, we find similarly to the integrated analysis a factor (Factor1) that strongly correlates with the 'CRP' value. This pattern describes the main source of variation within the proteomics data and is also aligned with some of the variation in the other datasets as captured by 'Factor3' in the integrated analysis (Figure 7C). However, a similar pattern as indicated by Factor2 that captures the time course of inflammation in the integrated analysis cannot be identified solely based on proteomics data.
The introduced workflow and the MOFA9 model itself are highly customizable with many adjustable parameters. Therefore, it is important to visualize and systematically compare the results produced by different configurations. To facilitate this task, the final output that can be generated by the workflow is a comparison of different named runs of the pipeline with different parameters in the pre-processing and model estimation. For example, the MOFA model can be estimated with different numbers of latent factors (Supplementary Figure 2A) or views with a lower number of features can be weighted (Supplementary Figure 3A). Configuring and running the last script of the workflow '07_Compare_Models' produces several plots to assess the similarity between different pipeline runs. FIG07_Variance_Model_Comparison (Supplementary Figure 2B, Supplementary Figure 3B) shows a comparison of the total explained variance for each view for different runs. The correlation of the factor values and feature factor weights between the different runs can indicate how much results change when modifying a certain parameter (Supplementary Figure 2C, Supplementary Figure 3C). Here, modifying the number of factors only causes minor changes in the estimated factor values and feature weights (Supplementary Figure 2C). Modifying the weighting of the data view results in much higher explained variance in the views with a lower number of features, e.g., the 'clinical' view (Supplementary Figure 3B). Nevertheless, the relevant features within the first three factors are still highly correlated to those inferred with the unweighted version (Supplementary Figure 3C).
With the generated model output .csv files in the results folder (e.g., the estimated factor and feature weights), further individual downstream analyses can be conducted. All code and necessary configuration files (including documentation) are available on GitHub at https://github.com/heiniglab/mofa_workflow. The singularity image that was created to enable an easy installation of the required conda packages for the analysis can be downloaded from https://doi.org/10.5281/zenodo.10815146. A small example dataset that can be used to perform an initial test of the pipeline can be downloaded from the same zenodo record as well.

Figure 7: MOFA output analysis. After the running of the MOFA model (03_Run_MOFA.ipynb) and the downstream analysis of factor values (04_Downstream_Factor_Analysis.ipynb) several plots are generated: (A) FIG03_Overview_Variance_Decomposition: returns a visualization of the explained variance of the estimated MOFA factors within the different views. Heatmap (left): shows the percentage of total variance of a view captured by a factor for each view. Barplot (right): shows the total percentage of variance that is captured by all factors for each view. (B) FIG04_Factor_Association_Numerical_Features: shows the Pearson correlation of the factor values with chosen numeric sample covariates, here: clinical variables (CRP, CK). (C) FIG04_Factor_Association_Categorical_Features: shows the difference in factor values for categorical sample covariates as a boxplot. Here, the factor values of Factors1-3 for each timepoint of ACS and Control patients are compared. Please click here to view a larger version of this figure.

Figure 8: MOFA feature analysis. After the running of the downstream analyses (04_Downstream_Factor_Analysis.ipynb, 05_Downstream_Investigate_Features.ipynb) several plots are generated. All plots here visualize MOFA Factor 2: (A) FIG04_Top_Feature_Overview_per_Factor: Heatmap (left) shows for each view the percentage of variance that is captured by the selected factor. Barplots (right) indicate the relevance of the features of the different views for the factor. On the left the total amount of features of a specific view within the top 1% highest-ranking features across views on the factor is given. On the right, the percentage is given, dividing the total number among the top 1% by the total number of features of that view. (B) FIG05_Heatmap_Feature_Overview: Heatmap (left) shows for the highest ranking 1% of features of the CD4.TCM cell-type the normalized expression values of each sample comparing the 'Control' group patients (CCS and non-CCS) to the different time points for 'ACS' patients. The barplot (right) shows the weight of the features. The direction of the sign of the weight is indicated before on the left before the cell-type names: '+' positive factor weight; '-' negative factor weight. (C) FIG06_Pathway_and_Genes: shows the weight of the highest 25% ranking genes for the factor that belong to enriched Interleukin pathways. In the heatmap on the top, they are averaged across views, and in the heatmap on the bottom shown per view. Please click here to view a larger version of this figure.
Supplementary Figure 1: Data harmonization effects. The figure shows FIG03_Overview_Variance_Decomposition for several different data-preprocessing configurations: visualization of the explained variance of the estimated MOFA factors within the different views. Heatmap (left): shows for each view the percentage of total variance of a view that is captured by a factor. Barplot (right): shows for each view the total percentage of variance that is captured by all factors. (A) The configuration ('MI_v1') based on which biological downstream results have been analyzed in previous figures (parameters set as in default configuration files in the cloned repository). (B) The same pre-processing configuration as in 'MI_v1' with the modification that no feature wise quantile normalization is applied (parameters set as in exemplary configuration files in 'config_examples' folder of repository). A screenshot of the MOFA model output warning for this configuration is added to the plot below. (C) The resulting variance decomposition when no pre-processing steps are applied, and all data is used as input without any pre-processing or filtering of features (parameters set as in exemplary configuration file in 'config_examples' folder of repository). A screenshot of the MOFA model output warning for this configuration is added to the plot below. Please click here to download this File.
Supplementary Figure 2: MOFA configuration - Factor amount effect. The resulting figures generated by the '07_Compare_Models.ipynb' script using several different configurations to run the MOFA model. (A) '03_MOFA_configs.csv': Example of the different configurations used to run the '03_Run_MOFA.ipynb' script specifying several different amount of factors (10,15,20,25). '07_Comparison_configs.csv': Example of how to specify the configuration input file for the execution of the script '07_Compare_Models.ipynb'. (B) 'FIG07_Variance_Model_Comparison' showing the total explained variance for each view (y-axis) for the different models across all factors specified in the model. (C) 'FIG07_Factor_Correlations' showing the correlation of the factor sample values between the different configurations. Please click here to download this File.
Supplementary Figure 3: MOFA configuration - Effect of weighting views. The resulting figures generated by the '07_Compare_Models.ipynb' script using several different configurations to run the MOFA model. (A) '03_MOFA_configs.csv': Example of the different configurations used to run the '03_Run_MOFA.ipynb' script specifying the 'weighting_of_views' parameter to be either 'TRUE' (MI_v1_MOFA_weighted) or 'FALSE' (MI_v1_MOFA). '07_Comparison_configs.csv': Example of how to specify the configuration input file for the execution of the script '07_Compare_Models.ipynb'. (B) 'FIG07_Variance_Model_Comparison' showing the total explained variance for each view (y-axis) for the different models across all factors specified in the model. (C) 'FIG07_Feature_Correlations' showing the correlation of the feature factor weights between the different configurations. Please click here to download this File.
Supplementary Figure 4: Multi-omic integration effect - using only proteomic data. The resulting patterns captured by the latent factors when only using proteomics data as input. (A) FIG04_Factor_Association_Numerical_Features: Pearson correlation of the factor values with clinical variables (CRP, CK). (B) FIG04_Factor_Association_Categorical_Features: Boxplot comparison of the factor values of each timepoint of ACS and Control patients. Please click here to download this File.
Supplementary File 1: Supplementary_File_
Running_Pipeline_with_Exemplary_Data. Descriptions on how to run the pipeline on the example data and the expected outputs are given in an additionally provided supplementary file. Please click here to download this File.
Supplementary Video File 1: Screen capture video of the protocol. Please click here to download this File.