IntroductionUsing multiple “omic” measurement technologies on the same large set of biological specimens or samples is a fundamental design principle for studying biological variation1,2,3,4. With thousands or millions of features measured, the manifold hypothesis posits that the data concentrate along an (unknown) low-dimensional manifold inside the high-dimensional data space. Biological variation then corresponds to different positions on that manifold. Empirically, it has been found that useful coordinates on that manifold—what we might call interesting axes of biological variation—can be constructed by linear summary statistics. For instance, in principal component analysis (PCA) or its extension to factor analysis, each principal axis, or factor, is a linear combination of features, and this (weighted) averaging has the effect of reducing noise compared to looking at variation only by measuring individual features. Since typically a relatively small number of top principal axes or factors is considered, such approaches may also be viewed as dimension reduction techniques.To take into account the multimodal structure of data, there are extensions of these methods, including canonical correlation analysis (CCA)5, group factor analysis (GFA)6, multi-omics factor analysis (MOFA)7, iCLuster8, and MuVI9.Applied to data, these methods return factors that, typically, fall into two categories: unsurprising and surprising. Unsurprising factors include those that align with known covariates (e.g., in the case of tissue biopsies, disease diagnosis, age or sex of the patient; in the case of environmental samples, geographical sampling location), or to technical artefacts10 (batch effects) associated with known experimental metadata. Surprising factors have no such easily identifiable explanation and are candidates for scientific novelty. The above-mentioned methods are unsupervised, in the sense that they do not use such already known information, and users are required to make their surprisingness considerations post-hoc11, in a process that is informal and error-prone.Here, we present Semi-supervised Omics Factor Analysis (SOFA), a probabilistic factor model that jointly models multi-modal omics data and sample-level information (hereafter: guiding variables). SOFA approximates the data with a low-rank matrix decomposition and partitions the latent factors into guided factors, each associated with a guiding variable, and unguided factors. In this way, already expected biological and technical drivers of variation can be regressed out from the latent space, making the interpretation of the remaining “surprising” factors simpler and potentially more fruitful.Our model takes inspiration from prior work in the unimodal setting. Principal component regression (PCR) uses PCA factors as covariates in a regression model to predict guiding variables12. However, PCR does not guarantee that its covariates are relevant, since for their construction, it only considers their covariance structure, but not the guiding variables. Partial least squares (PLS) simultaneously decomposes a high-dimensional set of predictors and guiding variables, aiming to find factors that capture the maximum covariance between the two13,14. SPEAR15 and DIABLO16 enable the incorporation of guiding variables of interest in the multimodal setting. Like PLS, these methods primarily uncover “unsurprising” factors that align with the guiding variables, rather than explicitly identifying “surprising” or unexpected factors.We present four applications of SOFA that highlight its broad utility: to the pan-gynecologic cohort of The Cancer Genome Atlas (TCGA)17, to the cancer dependency map18,19 (DepMap), to single-cell multi-omic data of microglial cells20, and to a collection of data sets of acute and chronic human heart failure21,22,23,24,25,26. We provide the SOFA method as an open-source Python package along with comprehensive tutorials.ResultsThe SOFA modelThe input to SOFA is a set of numeric data matrices containing multi-omic measurements from the same or mostly overlapping samples, and one or more covariates for each sample. The covariates can be continuous, binary or categorical and contain information that is known or suspected to be important for driving biological or technical variation between the samples. We also call these covariates sample-level guiding variables, and the different types of omics data views. SOFA models each omics view with a Gaussian likelihood and therefore expects omics data appropriately normalized and transformed. For the guiding variables users can choose between Bernoulli (binary), Gaussian (continuous) and Categorical (categorical) likelihoods. From this input, SOFA extracts a lower-dimensional representation comprising a shared factor matrix and modality-specific loading matrices (Fig. 1a). SOFA partitions the latent factors into guided factors each associated with a guiding variable, and unguided factors (Fig. 1b). The loading weights of each factor represent, for each view, the molecular features that it depends on, and can usually be interpreted in terms of biological processes, gene expression programs, modes of regulation, and the like (Methods, Fig. 1c). The factor matrix positions each sample in the latent space spanned by the factors, and each sample’s coordinates along the various factors can be interpreted as a degree of activity of a factor in the sample, similarly as in PCA or MOFA.Fig. 1: Overview of the SOFA model and downstream analyses.Full size imagea: The SOFA model. Input data modalities X are decomposed into shared latent factors Z and modality-specific feature loadings W. A hierarchical horseshoe prior shrinks the loading weights on the view, factor and feature level. This leads to factors shared between all views, factors active in one view and factors shared between subsets of views. b: Guided factors are linked to the guide variables. c: Factors not linked to a guide variable capture potentially novel sources of variation. Factors and loadings can be used for various sample- and feature-level downstream analyses.SOFA is a hierarchical Bayesian factor model, and we ensured scalability of its inference to large multi-omic data sets by employing stochastic variational inference (Methods). Our inference method for SOFA employs a sparsity-inducing prior27, which pulls the loading weights of some factors for some modalities to zero if the data are not pulling them elsewhere strongly enough. As a result, we can distinguish between variation that is commonly shared across modalities and variation that is present in only one or a few modalities.SOFA identifies cancer type-independent immune infiltration vs proliferation axis in the TCGA pan-gynecologic data setWe applied SOFA to a pan-gynecologic cancer multi-omic data set of the Cancer Genome Atlas (TCGA) project, which profiled the transcriptome (mRNA), proteome, methylome, and miRNA of 2599 samples from five different cancers17,28,29,30,31. Additionally, the study includes data about mutations, metadata, and clinical endpoints, progression-free interval (PFI) and overall survival (OS).We performed a SOFA rank-12 decomposition, using five guided factors with binary labels for each cancer type and seven unguided factors. Guided factors showed a clear separation based on the binary cancer type labels (Fig. 2a). We then verified that the unguided factors did not capture cancer type-related variance by computing Adjusted Rand Index (ARI) scores to assess clustering alignment with cancer labels (Figure S2e, Methods). This analysis showed that the set of unguided factors had the lowest alignment with cancer type labels, suggesting that the guided factors mostly explained the cancer-associated variability. Of note, ARI scores for the unguided SOFA decompositions at both rank-7 and rank-12 exhibited lower cluster alignment in comparison to the guided factorization, indicating that factor guidance helped disentangle cancer-associated variability (Figure. S2a, c, e).Fig. 2: SOFA factors are significantly associated with survival and capture immune infiltration vs proliferation axis.Full size imagea: Distributions of guided factors colored by their guiding variables (orange). b: Top: Hazard ratios from univariate Cox regression for each factor on overall survival (n = 2582). Hazard ratios > 1 indicate that positive factor values are associated with higher risk. Error bars indicate the 95% confidence interval. Hazard ratios are colored if adjusted p values are lower than 0.01 (log-likelihood ratio test). P values were adjusted for multiple testing using the Benjamini-Hochberg (BH) method. Adjusted p values were: 1.31e-20 (Factor 1), 0.35 (Factor 2), 4.5e-47 (Factor 3), 5.66e-9 (Factor 4), 1.15e-14 (Factor 5), 0.35 (Factor 6), 0.49 (Factor 7), 0.055 (Factor 8), 0.001 (Factor 9), 0.0006 (Factor 10), 0.58 (Factor 11), 0.031 (Factor 12). Bottom: Variance explained by each factor for each modality. c: Kaplan-Meier curves showing progression-free survival (left) and overall survival (right) for samples with factor values 95th percentiles of Factor 10. P values were computed using a two-sided log-rank test (n = 2582). d: Features with the highest absolute loading weights for factor 10 in transcriptomics (left) and proteomics loadings (right). e: Gene set overrepresentation analysis (one-sided hypergeometric test) for the transcriptomics features with the highest absolute loadings of factor 10. The -log10 adjusted p values are multiplied by the sign of the loadings to reflect the direction of the association with the factor. Cancer types shown: brca: breast invasive carcinoma; cesc: Cervical squamous cell carcinoma and endocervical adenocarcinoma; ov: Ovarian serous cystadenocarcinoma; ucec: Uterine Corpus Endometrial Carcinoma; ucs: Uterine Carcinosarcoma. Source data are provided as a Source Data file.To evaluate the importance of each factor, we calculated the fraction of variance explained by each factor across the input modalities (Fig. 2b). For example, Factor 3 explained the highest fraction of variance in the miRNA data, consistent with previous research suggesting miRNA expression is critical in ovarian cancers and may serve as a biomarker32,33. Factor 1, representative for breast cancers, showed breast cancer-associated proteins AR34, PR35 and ER-alpha36 among the largest proteomics loading weights (Figure S2f). Factor 3 (ovarian serous cystadenocarcinoma) featured BCAM37, MALAT138and WT139 among the highest transcriptomics loading weights (Figure S2g).We then explored the SOFA factors for potential biological relevance by testing them for association with disease outcome, using univariate Cox regression of each factor on overall survival (OS, Fig. 2b). Notably, Factor 3, linked to ovarian serous cystadenocarcinoma, was significantly linked to poor survival (Hazard Ratio (HR) = 2.44 95% CI: 2.17–2.75). Two out of seven unguided factors showed significant associations with OS, suggesting that most of these factors do not strongly influence survival outcomes.Among the unguided factors predictive of survival, we further investigated Factor 10, which had the highest HR among the unguided Factors (HR = 1.56 95% CI: 1.25–1.94, Fig. 2b). Stratifying samples by Factor 10 values, followed by non-parametric survival curve estimation, revealed substantial differences in median survival times between groups (OS: Median survival time 1701 vs. 2888 days, PFI: 1396 vs 3669 days, Fig. 2c). High values of Factor 10 were associated with increased expression of proliferation-related genes and oncogenic pathways, while low values were enriched for immune-related signatures, suggesting an antagonism between tumor proliferation and immune activity. Negative loadings of Factor 10 featured NTN4, CD59 and CFB, and proteins including PDCD4, indicative of complement system activation and lymphocyte infiltration40,41,42,43 (Fig. 2d). In contrast, positive loadings included proliferation markers TOP2A and TPX2, as well as Cyclin B1, FoxM144,45,46,47, suggesting a signature of increased proliferation. This was further corroborated by performing a gene set overrepresentation analysis on the top transcriptomics loadings. Negative loadings were associated with immune system related gene sets, while those with positive loadings were enriched in cell cycle and proliferation related terms (Fig. 2e). Finally, we note that Factor 10 was not restricted to a single cancer type (Figure S2b, d), indicating that it captured biological signatures common to all cancers. Taken together, we found that the guided SOFA factorization delineated cancer type-specific variation from broader biological processes, such cell proliferation or immune cell infiltration, thereby facilitating their identification and subsequent downstream analyses.SOFA identifies EGFR dependency axis in pan-cancer dependency map, while accounting for known driver mutations and cell growthNext, we analyzed a pan-cancer cell line data set of the Cancer Dependency Map (DepMap)19,48 together with proteomics data from the ProCan-DepMapSanger study18. The DepMap project aims to identify cancer vulnerabilities and drug targets across a diverse range of cancer types. The data set includes multi-omic data, encompassing transcriptomics49, proteomics18, and methylation50, as well as drug response profiles for 627 drugs50, and CRISPR-Cas9 gene essentiality scores48 for 17,485 genes for 949 cancer cell lines across 26 different tissues.To identify CRISPR-Cas9 gene essentiality scores associated with multi-omic features, we performed a rank-20 SOFA decomposition to the DepMap data, while accounting for potential/known drivers of variation such as growth rate, microsatellite instability (MSI) status, BRAF, TP53, and PIK3CA mutation and hematopoietic lineage. Although CRISPR-Cas9 essentiality scores could technically be included as an additional view in the SOFA decomposition, we chose to exclude them to enable independent testing of factor associations with these scores.A t-SNE of all the inferred factors showed that they captured the differences between cancer types (Fig. 3a). Further, we confirmed successful factor guidance by asserting that SOFA extracted high or low factor values in a variable-dependent manner (Fig. 3c). We then focused on the guided factors and observed that Factor 6, consistent with previous studies identifying hematopoietic lineage as a major driver of variation in cell lines18, explained the largest fraction of variance (34%) across all data modalities (Fig. 3b). Factor 1, guided by growth rate, mostly explained variance in drug response (Fig. 3b), which is expected as the natural logarithm of the half maximal inhibitory concentration (lnIC50) drug response is derived from the growth rate of the cell line. We then subjected the loadings of guided factors to gene set overrepresentation analysis. Expectedly, Factor 3, guided with labels for BRAF mutation, a common mutation in skin cancer51, was associated with skin cancer related gene sets (Figure S3c). The loadings of Factor 6, guided with hematopoietic lineage, displayed gene sets related to the immune system (Figure S3d). Statistical testing revealed that Factors 3 and 4 were significantly associated with CRISPR-Cas9 essentiality scores of BRAF and TP53, the mutations used to guide these factors (Figure. S3e, f, Methods).Fig. 3: SOFA factors capture lineage subtypes and EGFR dependency axis in DepMap.Full size imagea: tSNE of all SOFA factors, colored by cancer types. b: Variance explained by each factor for each modality. c: Distributions of guided factors colored by their guiding variables (orange). d: Distribution of Factor 10 cell lines from the respiratory system are colored (left). Distribution of Factor 10 for cell lines from non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC) (right). e: Highest absolute transcriptomics (left) and drug response (right) loadings of Factor 10. f: Heatmap showing CRISPR essentiality scores (ES) of top 10 genes significantly associated with Factor 10, tested using a linear model (two-sided t-test of the slope coefficient; n = 778). The cell lines are ordered based on their values of Factor 10. g: Distribution of Factor 12, cell lines from the circulatory system are colored (left). Distribution of Factor 12 for B-cell, T-cell and myeloid cell lines (right). h: Volcano plot of associations between Factor 10 and mutations, tested using a two-sided two-sample t-test (n = 778). Effect sizes are the differences in means between mutated and unmutated samples. BH-adjusted p values