CNAmet (Louhimo and Hautaniemi, 2011) is a biomarker-discovery correlation-based method. It consists of two main steps; first, weights are calculated for each gene connecting it to DNA methylation and copy number variation (CNV). Second, each gene’s weights are combined and tested (using a corrected p-value via permutation) to calculate a global score for each gene. These scores help identify whether a gene is hypomethylated (and upregulated) or hypermethylated (and downregulated). The main hypothesis is that amplified copy number and hypomethylation result in gene upregulation. iPAC (in-trans Process Associated and cis-Correlated) (Aure et al., 2013) is an unsupervised, integrative method based on mRNA, and CNV aims to identify the cis-regulated genes. It also uses a sequence of statistical tests to narrow down the list of cancer driver genes. In summary, it takes the matrix of all the genes as the input, first filters the genes based on aberration frequency (>10%), then filters the remaining genes based on in-cis correlation (>0.6), and finally checks the in-trans functionality for the remaining ones that make the final set of gens. MEMo (Mutual Exclusivity Modules) (Ciriello et al., 2012) is a module-discovery method to find a set of genes that exhibit the same genetic alternation among patients. First, it gives a score to each gene and makes a binary-event-matrix based on these scores where its elements are either “1” indicating that the gene is significantly altered or “0” otherwise. Subsequently, it builds a network from genes involved in the same molecular pathway (using curated and nun-curated sources of biological information/interactions). The final step collects the genomic events within this network that show a significant level of mutual exclusivity via a permutation test.

Illustrative Case-studies: CNAmet is applied to a cohort of glioblastoma multiforme (GBM) patients from TCGA to find synergistically regulated genes by DNA methylation CNV. It could identify this synergistic effect on well-known oncogenes, including MDM2, EGFR, and PDGFRA. CNAmet also showed that GBM patients with hypomethylated (upregulated) EGFR had a better prognosis than patients with amplified EGFR. iPAC is applied to a cohort of breast carcinoma patients. It identified a list of significant genes, including ERBB2, MAP3K7, MDM4, FGFR1, CCND1, and FADD, which are well-known cancer-associated genes. It also included some less appreciated genes such as ATAD2, TPD52, and PPM1D, which were reported as cancer genes in previous independent studies (Choschzick et al., 2010; Lambros et al., 2010). iPAC could also identify several novel genes such as MTL5 that can affect multiple proteins/enzymes via its negative correlation with the MT (metallothionein) family of proteins and metal-binding ability.

CCA-based methods (Dolédec and Chessel, 1994) can be applied for module identification, feature selection, and classification in high-dimensional multi-Omics data. Due to the high dimensionality of Omics data, standard CCA cannot be employed directly. Therefore, there have been several extensions of CCA for more than two datasets ( B ≥ 3 ) and/or high dimensional data in an unsupervised setting. Sparse MCCA (Sparse Multiple Canonical Correlation Analysis) (Witten and Tibshirani, 2009) is a sample-subtyping method applicable for more than two data types. It aims to find sparse components (linear combination of the features) by maximizing the following objective function ( X ( b ) s are standardized beforehand): a r g m a x w ( b ) ∑ b w ′ ( b ) X ′ ( b ) X ( b ) w ( b ) , s . t . ‖ w ( b ) ‖ 2 ≤ 1 , P b ( w ( b ) ) ≤ c b , where P b is a convex penalty function, such as lasso or fused lasso, and c b is a tuning parameter for data type b ( b = 1 , … , B ). Two main disadvantages of CCA-based models are (1) they are not capable of handling non-linearity in data, which most of the time is the case in real-world data, and (2) they cannot fully take into account the structural information between and within Omics data (such as gene-gene interaction and PPIs).

CIA-based methods (Dolédec and Chessel, 1994; Dray et al., 2003) is another approach to find the low-dimensional components in two-table data settings where X = [ X n × p 1 ( 1 ) | X n × p 2 ( 2 ) ] . This method was first introduced in ecology to link species abundance with environmental features. Orthonormal directions ( u and v ) are computed by maximizing the covariance between the data tables: a r g m a x u , v     u ′ X ( 1 ) ′ X ( 2 ) v , s . t . ‖ u ‖ = ‖ v ‖ = 1.

CIA can be considered as a variation of CCA (Sankaran and Holmes, 2019); the only difference is that in CIA the norm constraint ( ‖ u ‖ = ‖ v ‖ = 1 ) is directly applied on the orthonormal directions ( u and v ). MCIA (Multiple Co-Inertia Analysis) (Meng et al., 2014) is the extension of the CIA for the analysis of more than two data tables. sMCIA (sparse Multiple Co-Inertia Analysis) (Min and Long, 2020) is a sparse version of MCIA via imposing a sparsity constraint on the transformed direction vectors.

Illustrative Case-studies: Sparse MCCA is applied on a diffuse large B-cell lymphoma dataset to assess the relationships (such as co-amplification and codeletion) between copy number changes in genome regions on separate chromosomes (Meng et al., 2014). The results showed a complex relationship between CAN in different chromosomes. MCIA is applied to gene and protein expression for NCI-60 cancer cell lines from different tissues. Results showed that different cell lines were differentiated based on their tissue of origin. That is, cell line-specific features can help improve prediction and biomarker identification. In the 2nd application (Meng et al., 2014), MCIA is applied to a cohort of ovarian cancer patients, including mRNA expression data obtained from microarray and NGS. It identified four known subtypes of ovarian cancer (proliferative, immunoreactive, mesenchymal, and differentiated) along with its first two directions (components). Moreover, gene expression analysis in each component showed the capability of MCIA to detect disease subtype-specific markers.

MOFA (Multi-Omics Factor Analysis) (Argelaguet et al., 2018) is an unsupervised multi-Omics integration method that aims to detect the sources of variation (both technical and biological) in datasets via latent factors. It first decomposes each Omics data ( X ( b ) ) as follows: X ( b ) = Z W ( b ) + ε ( b ) , b = 1 , … , B where Z indicates the factor-matrix which is common for all data types and W ( b ) is a matrix of weights for datatype b , and ε ( b ) is the residual (or noise) for datatype b . It then—following the Bayesian framework—assigns a prior distribution for Z , W ( b ) , and parameters of the noise term. MOFA then applies a two-step regularization on the weight matrices to deal with the high dimensionality of multi-Omics data. it first identifies which factor is more active in which datatype (Omics type) and then applies a feature-wise sparsity to find a smaller set of features with active weights. These latent factors can serve as an input for further downstream analysis, including sample classification and missing data imputation. The most important advantages of MOFA are its interpretability, the ability to visualize samples in the factor space, and the capability of handling missing data and data with different distributions.

Illustrative Case-studies: MOFA is applied to a cohort of patients with chronic lymphocytic leukemia (CLL) to integrate mRNA expression, DNA methylation, somatic mutation, and drug response. It identified two important (already-known) markers, including the IGHV gene (immunoglobulin heavy-chain variable) and trisomy of chromosome 12. However, MOFA could find a more comprehensive and complex sub-structure for IGHV and connect it with multiple Omics, including changes in mRNA expression (LPL, PLD1, ADAM29), DNA methylation (cg17479716, cg19358877, cg26615224), and with drugs (tamatinib, dasatinib, AZD7762) that target kinases in the B-cell receptor pathway (Argelaguet et al., 2018). These changes in mRNA expression and DNA methylation were previously connected to IGHV in different independent studies (Plesingerova et al., 2017). Interestingly, IGHV and trisomy of chromosome 12 explained only <20% of the variation in CLL patients, indicating the presence of other factors and sources of heterogeneity. Therefore, they could find the oxidative stress pathway (with HSP family of proteins as the top-weighted genes) as one of the critical drivers which was previously underappreciated in the context of CLL. The results (factors) of MOFA are then used in a Cox-PH regression model and could predict the time to the next treatment with a reasonably high prediction accuracy (C-index∼75%). In the second and the third applications, MOFA is used to analyze Ustekinumab (UST) drug-response (Verstockt et al., 2019) and mESCs (mouse embryonic stem cells) multi-omics data (Argelaguet et al., 2018) to identify predictive factors (a combination of different Omics data).

The input data in the kernel-based methods is the kernel matrix ( k ( x i ( b ) , x j ( b ) ) ,   b = 1 , … , B ), also called inter-patients similarities, instead of the original data ( X = [ X n × p 1 ( 1 ) | … | X n × p B ( B ) ] ). Therefore, the multi-Omics data integration problem is converted to kernel integration in the sample space ( R n ) rather than the multi-Omics (feature) space ( R n × ( p 1 + … + p B ) ). As a result, the optimization problems in the kernel-based methods are called dimension-free, i.e., it does not depend on the total number of the features ( p 1 + … + p B ).

SNF (Similarity Network Fusion) (Wang et al., 2014) is a popular method for multi-Omics data integration and subtype analysis. It first builds a sample-by-sample similarity matrix (or network, where nodes are samples and edges are similarities between samples) for each dataset separately and then fuses them to a global (weighted) sample similarity network. The second step (network-fusion) uses a nonlinear message-passing theory-based method (Pearl, 2014) to fuse the similarity matrices. SNF may lead to false fusion since it does not distinguish between different data types. Another drawback of SNF is that it uses Euclidean distance to calculate the similarity matrices between the samples that often is incapable of capturing the intrinsic similarities between data points. To address this issue, DSSF (Deep Subspace Similarity Fusion) (Yang et al., 2018) employs an auto-encoder to improve the discriminative similarity between samples. AFN (Affinity Network Fusion) (Ma and Zhang, 2018) is also an extension of SNF that enables the consideration of patients' pairwise distances. To handle the unmatched samples (different sample sizes in different Omics-types), NEMO (NEighborhood based Multi-Omics clustering) (Rappoport and Shamir, 2019) is introduced that enables the computation of global kernel matrix without performing any imputation on the missing observation. INF (Integrative Network Fusion) (Chierici et al., 2020) is another extension that utilizes SNF within a predictive framework including RF (Breiman, 2001) (Random Forest) and LSVM (Cortes and Vapnik, 1995) (Linear Support Vector Machine). See Supplementary Appendix Section SA1.1 for more details and case studies.

Illustrative Case-studies: SNF is applied to different multi-omics data integration studies. In the 1^st application (2018) (Chiu et al., 2018), it is applied on a cohort of triple-negative breast cancers (TNBC) patients from TCGA (including CNV, miRNA, and mRNA expressions) to identify the different sub-groups of cancer patients. Results revealed a new TNBC classification scheme with three different clusters of patients. One of the clusters, interestingly, was enriched in the “non-basal” subtype (by PAM50), whereas PAM50 obtained the most common “basal-like” subtype. This nan-basal cluster showed more aggressive clinical characteristics and distinctive oncogenic features (including 38% basal-like2 and 50% luminal androgen receptor subtypes).

Standard factorization methods commonly use singular value decomposition (SVD), such as PCA. However, for some data types, such as genotypes, the original matrices are non-negative. SVD-based factorizations contain negative entries, making it difficult to interpret their results in some applications. In contrast, NMF (Nonnegative Matrix Factorization) (Lee and Seung, 2001) restricts the entries in matrix factors to be non-negative.

iCluster (Shen et al., 2009) simultaneously considers the association between different data types and the covariance structure within each datatype. It employs the principles of two methods, including probabilistic PCA(Tipping and Bishop, 1999) and a (spectral) relaxed version of k -means (Zha et al., 2001). It first uses the Gaussian latent variable model to compute the posterior mean of the (latent variable) components ( T ). It then calculates the class membership by employing the standard k -means algorithm. The integrative model can be written as follows: X ( b ) = W ( b ) ′ T + ε ( b ) , f o r   b = 1 , … , B , T ∼ N ( 0 , I ) , ε ( b ) ∼ N ( 0 , σ ( b ) ) , where W ( b ) is the ( c − 1 ) × p b matrix of coefficient (also called loading matrix), c is the number of the clusters, and T is the ( c − 1 ) × n matrix of latent variable components that are shared between the data tables and explains the correlation between the different data types ( n is the sample size); ε ( b ) indicates the remaining (unexplained) variances for each data type. Therefore, the final data matrix will become as X = ( X ( 1 ) , … , X ( B ) ) ∼ N ( 0 , W W ′ + σ ) . Authors have applied a lasso-based penalty on W in the final likelihood function. Cluster memberships are then calculated by applying k -means clustering on the posterior mean of the latent variable components ( E ( T | X ) ). iNMF (integrative Non-negative Matrix Factorization) (Yang and Michailidis, 2016) is a multi-table extension of NMF to account for heterogeneity between the multiple datasets by providing heterogenous estimations/combinations ( V b T b ) via minimizing the following loss function (using partitioned factorization structure): a r g m i n W , V b , T b ∑ b = 1 B ‖ X − ( W + V b ) T b ‖ F 2 + λ ∑ b = 1 B ‖ V b T b ‖ F 2 , f o r   b = 1 , … , B s . t . W ≥ 0 , T b ≥ 0 , V b ≥ 0 , where V is a homogeneity parameter and enables to account for different degrees of heterogeneity in the multiple datasets (since larger values of V result in smaller heterogeneous components V b T b ). ‖ . ‖ F is the Frobenius norm (Lee and Seung, 1999). The authors also adopted the sparse version of the iNMF by applying L 1 − n o r m to T b . Whereas iCluster-based methods, NMF-based methods do not rely on any model assumptions and allow each sample to fall in more than one class or be excluded from the classification (see Supplementary Appendix Section SA1.2 for more NMF-based and iCluster-based methods and illustrative case-studies).

Illustrative Case-studies: iCluster is applied in different studies, mainly for cancer subtyping. It is recently applied to a cohort of ovarian carcinoma patients (including CNV, DNA methylation, and mRNA expression) to identify prognostic biomarkers (Zheng et al., 2019). The results revealed three distinct clusters of samples and identified UBB (ubiquitin B) and IL18BP (interleukin 18 binding protein) genes as the most prognostic biomarkers. The results suggested that lower expression of these two genes may result in higher methylation and lower CNV. Therefore, evaluating the expression of these two genes can help in early tumor diagnosis. In another study, iCluster is applied to multi-Omics data of adult soft tissue sarcomas (Lazar et al., 2017). iCluster showed that SS-subtype (synovial sarcoma) was the most distinct sarcoma with partial/complete loss of chromosome 3p (45% of cases), high expression of FGFR3 and miR-183, and methylation of the PDE4A promoter. Another cluster identified by iCluster mainly included LMS (Leiomyosarcoma) cases with high expression of MYLK, MYH11, ACTG2, miR-143, miR-145, lower inferred activity of the apoptosis pathway, and higher hormone receptor (ER/PR) levels. It inferred PI3K/AKT pathway activity. The authors also concluded that copy number changes were the most informatics Omics in characterizing these sarcomas (except SS).

In the Bayesian framework of a clustering task, class memberships are calculated using a probability model (such as a Dirichlet Process Mixture (DPM) model (MacEachern and Müller, 2000)) subject to a priori assumption about what the true relationship between the data might be, which is expressed as a probability distribution. This probability is then updated as new observations become available (which is captured by a posterior distribution). This approach enables the use of prior information informing the clusters (sub-samples, sub-disease, or sub-features). The DPM model is one of the most widely used Bayesian nonparametric methods in the multi-Omics (multitype) clustering-based data integration. For a tutorial on DPM models, refer to (Li et al., 2019).

LRAcluster (Low-Rank-Approximation) (Wu et al., 2015) is a low-rank probabilistic method (similar to iClusterPlus) for molecular classification that takes both continuous and categorical data as input. LRAcluster first models each datatype using a probabilistic model and combine them as L ( Θ ) = ∑ b = 1 B L ( Θ ( b ) , X ( b ) ) , for b = 1 , … , B , where Θ ( b ) is the parameter matrix for datatype b . Θ is the overall parameter matrix and is assumed to be a low-rank matrix that leads to the following optimization problem: a r g m i n Θ L ( Θ ) + λ | Θ | ∗ , where λ is the tuning parameter and | Θ | ∗ indicates the nuclear norm of Θ . An iterative, fast LRA of the parameter matrix is then applied to solve this optimization problem. The final clustering task is applied to the low-dimension subspace (using k -means) to find feature subtypes (refer to (Subramanian et al., 2020) for more details and examples). MDI (Multiple Dataset Integration) (Kirk et al., 2012), TMD (Transcriptional Modules Discovery) (Savage et al., 2010), PSDF (Patient-Specific Data Fusion) (Yuan et al., 2011), and BCC (Bayesian Consensus Clustering) (Lock and Dunson, 2013) are four closely related integrative methods that all adopt a DPM. However, MDI and BCC have the same objective (clustering and subtyping), and all can integrate more than two data types (for more details, see Supplementary Appendix Section SA1.3).

Illustrative Case-studies: LRAcluster is applied in a study of hepatocellular carcinoma (HCC) as the major subtype of liver cancer (Wang et al., 2019) to characterize the molecular alternation of the metastatic HCCs. The results identified a list of individualized molecules (including TNC, LAMA2, LAMC3, PDGFRA, CYP2E1, CYP3A4, CYP2C8, CYP1B1, CPS1, TAT, and HPD) significantly expressed between the primary tumor compare and portal vein tumor thrombosis. Therefore, an individualized differential analysis for sequencing data was proposed to automate the process of finding these individualized genes.

COCA (Cluster-Of-Cluster Assignment) (Hoadley et al., 2014) integrates the single-Omics clusters using hierarchical clustering based on pairwise concordance between different Omics platforms (including mRNA, miRNA, DNA methylation, and mutation). PINS (Perturbation clustering for data INtegration and disease Subtyping) (Nguyen et al., 2017) is a disease sub-typing method. It first partitions the samples into k ( k ∈ [ 2 … K ] ) clusters, then builds the patient connectivity matrices based on the pairwise connectivity for each possible cluster (see Supplementary Appendix Section SA1.4 for more information).

Illustrative Case-studies: COCA has recently been applied to a cohort of Ugandan cervical carcinoma patients (both H I V + and H I V − ) that is the first comprehensive profiling (genomic, transcriptomic, epigenomic) of sub-Saharan African patients (Gagliardi et al., 2020). They could identify human papillomavirus (HPV)-clade-specific (clade A7 and A9) patterns of multiple Omics features, including DNA methylation and gene expression. For instance, upregulated genes in clade A7-samples (such as PXDN) are also upregulated in cancers that progress through the epithelial-mesenchymal transition; and DNA methylation is closely regulated through cell differentiation. The clustering result showed the loss of E2 expression in the A7-enriched cluster due to HPV integration in clade A7-samples. However, the A9-enriched cluster showed partial HPV integration supporting the higher expression of episomal HPV genes (due to E2 expression) in these patients. Therefore, the authors hypothesized that clade A9-infected samples might have a more active HPV infection. In another application on glioblastoma cancer patients (Yuan et al., 2020), COCA could identify two novel subtypes, including HX-1 and HX-2 categorized by three CpG regions (∼DUSP1, PHOX2, HOXA7) and 15 gene mutations, including PCDH1, CYP27B1, LPIN3, GPR32, BCL6, OR4Q3, MAGI3, SKIV2L, PCSK5, AKAP12, UBE3B, MAP4, TP53BP1, F5, RHOBTB1.

NBS (Network-Based Stratification) (Hofree et al., 2013) is a sample-stratification method that uses both network propagation algorithm and matrix factorization to construct the final subtypes. Therefore, it can be categorized under either of these categories. NBS integrates genome-scale somatic mutations with a gene-interaction network. It first maps the mutations for each sample onto a gene-interaction network from STRING (https://string-db.org/), Pathway Commons (https://www.pathwaycommons.org/), and HumanNet (https://www.inetbio.org/humannet/download.php), and constructs the patient-by-gene matrix ( M 0 ). Then network propagation is used to smooth the sample-mutation-gene network as follows: M t + 1 = α   M t   A +   ( 1 − α ) M 0 , where A is a normalized adjacency matric of the gene-interaction network, α is a tuning parameter controlling the mutation diffusion. The smoothing (propagation) function runs iteratively till convergence. The result of this step is a network-smoothed profile where the elements indicate the network proximity of each gene to the mutated genes for a specific sample. FUSENET (Žitnik and Zupan, 2015) and DFMF (Žitnik and Zupan, 2014) are flexible about input data and their distributions. The latter does not treat the entire input data as a single matrix and therefore, enables the identification of data-specific factors. Medusa (Zitnik and Zupan, 2016) is a module-discovery method that partly uses the same methodology as DFMF (Data Fusion by Matrix Factorization) (Žitnik and Zupan, 2014) to construct a fused network (see Supplementary Appendix Section SA2.1 for more information). MAE (Multi-view factorization AutoEncoder) (Ma and Zhang, 2019) is a combination of matrix factorization and an autoencoder that enables the simultaneous embedding of both features (Omics) and samples via more complex nonlinear transformations. It first constructs an interaction graph for each datatype. To do so, the interactions among the feature in each datatype are represented as a network ( N ∈ R + p × p ). For instance, for proteome datatype, network N will be protein-protein interaction networks (PPIs) that are publicly available (such as Reactome https://reactome.org/). Note that MAE can also be categorized as a supervised (deep) neural networks method. RAIMC (RBP-AS Target Prediction Based on Inductive Matrix Completion) (Qiu et al., 2021) is based on inductive matrix completion (IMC), where integrated RNA-binding proteins (RBP) similarities were calculated based on RBP-regulating similarity and integrated alternative splicing (AS) event similarities were computed based on AS module-similarity. Then Gaussian interaction profiles (GIP) for RBPs and AS events are computed and combined using the fast kernel learning (FKL). Before completing the association matrix with IMC, a top-kk nearest neighbor model is applied to denoise the integrated similarity matrix. See (Ou-Yang et al., 2022) for a comprehensive review of matrix factorization methods for biomedical link prediction, including, IMCDriver (Zhang et al., 2021) and DisoFun (Wang et al., 2020).

Illustrative Case-studies: NBS is applied for patient-subtype identification and discriminating the somatic mutation profiles in uterine, ovarian, and lung cancer studies obtained from TCGA. The survival result based on the identified subtypes showed that ovarian cancer patients with the most aggressive tumor had a mean survival of 32 months compared to others (∼80 months). The fibroblast growth factor (FGF) signaling pathway was enriched for this sub-network of patients with the worst survival in concordance with previous studies indicating the FGF signaling pathway as a driver of tumor progression resistant to anti-VEGF therapy (Cole et al., 2010). The next subtype of patients with relatively better (higher) survival was mainly enriched in DNA damage–response genes (including ATM, ATR, BRCA1, BRCA2, RAD51, and CHEK2) that have been referred to as BRCAness in previous studies (Konstantinopoulos et al., 2010).

Bayesian networks (BNs) are a combination of (directed acyclic) graph/network theory and probability models. Suppose N = ( V ,   E ) is a network where V is a vector of nodes and E is the set of edges. The structure of N in BNs is a directed acyclic graph (DAG) that defines the factorization of the joint probability of V = { X ( 1 ) , … , X ( B ) } into a set of local probability distributions (one for each X ( b ) ) via Markov property (Korb and Nicholson, 2010) of BNs: P ( X ( 1 ) , … , X ( B ) ) = ∏ b = 1 B P ( X ( b ) | π X ( b ) ) , showing that each node (random variable X ( b ) ) directly depends only on its parents π X ( b ) (Scutari, 2009). The main disadvantage of BNs is their computational complexity since the number of network structures grows exponentially with the number of nodes. However, using the Monte Carlo Markov Chain (MCMC) approach can partially help the situation (Lin and Lane, 2017). PARADIGM (PAthway Recognition Algorithm using Data Integration on Genomic Models) (Vaske et al., 2010) can also be categorized as a BN approach. It uses the prior knowledge of the given pathways to model the nodes (Omics data). CONEXIC (Akavia et al., 2010) is another BN that aims to find cancer-driver mutations by integrating gene expression and CNVs.

Illustrative Case-studies: CONEXIC is applied to gene-CNV paired data from melanoma patients (Lin et al., 2008) to identify a list of cancer driver genes. First, a list of candidates was generated using CNV data, and then the most likely drivers were collected by integrating CNV and mRNA expression. It resulted in several modulators that explain the behavior of 7869 genes. Many of the top modulators were involved in melanoma-related pathways and included known oncogenes and tumor suppressors. CONEXIC could successfully pick known cancer-related genes out of a large region with many underlying genes. For instance, CCNB2 (cell-cycle regulator) was selected from a large, amplified region. Finally, an automated literature-mining method called LitVAn (literature vector analysis) was used to find overrepresented terms in published studies. It resulted in a few well-known activated features in melanoma (such as PI3K, MAPK) and a novel process called “RAB” (Rabs regulate vesicular trafficking).

Network Propagation (NP) ( Cowen et al., 2017) is a stochastic process that tracks each node’s flow and tries to amplify the signals through prior information and pass them to its neighborhoods over time. Suppose N = ( V ,   E ) is a network with an adjacency matrix W . Suppose p 0 ( v ) indicates the starting value of prior (known) information for node v ∈ V . For instance, it can be a vector of 0 and 1, 1 indicating the genes known to be related to the disease, and 0 otherwise. The value of p 0 ( v ) is the amount of information that we want to flow (diffuse) from each node to its neighborhoods. Therefore, the amount of information of node v (also called the state of node v ) at time t ( p t ( v ) ) can be formulated as the sum of the information of its neighbor ( N ( v ) ) at the previous time ( t − 1 ): p t ( v ) = ∑ u ∈ N ( v ) p t − 1 ( u )   w ( u , v ) , (1) where w ( u , v ) indicates the (normalized) weights between nodes u and v and is based on the relationship/interaction between these two nodes. The result of this iterative propagation process (for t times) is the gene-ranks ( p t ( v ) ). Eq. 1 can be re-written with matrix notation as follows: p t = W ∗ p t − 1 , where W ∗ is a transition matrix and calculated from the adjacency matrix W . The random walk with restart—RWR (Tong et al., 2008) is a propagation algorithm that allows a walker (an imaginary particle) to start a walk (flow) from the initial node v 0 ∈ V (with prior probability p 0 ) to node v t − 1 at a discrete-time step t − 1 ( ∼ p t − 1 ). It then walks from node v t − 1 to the next (randomly selected) neighbor v t by following a given transition matrix. Therefore, p t can be written as: p t = α p 0 + ( 1 − α ) W t − 1 ∗ , where α is called restart probability and controls the amount of prior information considered in the network, and W ∗ is a normalized transition matrix. Different algorithms may use different transition matrices. RWRM (Random Walk with Restart on Multigraphs) (Li and Li, 2012) is one of the first extensions of network propagation for integrating multigraph gene networks. It enables multiple edges between two nodes [see Supplementary Appendix Section SA2.2 for extensions of RWR algorithm, including RWR-M (Valdeolivas et al., 2019) and RWR-MH (Valdeolivas et al., 2019)]. TieDIE (Tied Diffusion through Interacting Events) (Paull et al., 2013) accepts a biological graph/pathway (such as PPIs or gene interaction networks) and a set of prior scores for each node indicating the involvement of each node in the network. SNF (Wang et al., 2014) can be considered as both a clustering and network-based method. We have discussed SNF and its extensions in Unsupervised Multi-omics Data Integration Methods.

Based on a benchmarking study (Picart-Armada et al., 2019) for network propagation methods, selecting a prominent network analysis method is not clear-cut. The authors concluded that network propagation methods enable the biomarker discovery, but their efficiency greatly depends on the input biological network and the nodesʼ initial score (see Supplementary Appendix Section SA2.2 for more information).

Illustrative Case-studies: SNF is applied to identify GBM subtypes vis integrating DNA methylation, mRNA, and miRNA expressions (Wang et al., 2014). The results indicated that most edges in the similarity network (patientsʼ similarities) were only detectable when two or more types of Omics information has applied. SNF could successfully distinguish the previously reported IDH subtype (Sturm et al., 2012) consisting of younger patients with an IDH1 mutation. SNF could further identify a subtype of patients who were more responsive to temozolomide, TMZ (a common GBM treatment), whereas another distinct subtype of patients with overexpressed CTSD and less responsive to TMZ [which is consistent with an in vitro study (Sun et al., 2012)]. SNF has recently been applied to TNBC (Chiu et al., 2018) and pancreatic cancers (Sinkala et al., 2020) to identify disease subtypes.

WGCNA (Weighted Gene Co-expression Network Analysis) (Langfelder and Horvath, 2008) is a gene-prioritizing correlation-network-based method that also enables gene module identification. Correlation networks are based on the correlation between a node and an outcome. The significance of a node (such as a gene) is then determined based on either the correlation coefficient or a regression-based p-value. WGCNA can be employed to find gene modules, sub-modules, and marker-prioritization. Lemon-Tree (Bonnet et al., 2015) is a biomarker-discovery method that first processes (normalize) the expression data (mRNA) and finds the co-expressed clusters of genes via a model-based Gibbs sampler (Joshi et al., 2008). It then employs ensemble methods (including spectral edge clustering algorithm) to identify gene modules and regulatory based on the co-expressed genes. Other Omics features (including miRNA, CNV, DNA methylation, and genotype) are added to the model as additional candidates and are combined to calculate the regulatory scores. Lemon-Tree also enables the gene ontology enrichment analysis for the modules. DBN (Deep Belief Network) (Liang et al., 2014) is a sample-classification (deep-learning, DL) method that integrates mRNA, miRNA, and DNA methylation data. DL methods are initially constructed from multi-layered (or deep) artificial neural networks (Bishop, 1995) (ANNs), inspired by actual NNs in the brain. ANN is a parallel system that accepts the input data in its first layer (input layer). It then passes the data into one or more hidden layers that ultimately connect them to an output layer. ANNs used for DL have more hidden layers where each of them helps to refine its previous layer by running a feature construction task. DBN applies a Gaussian restricted Boltzmann machines (Gaussian RBM) (Hinton, 2012) model to obtain the featuresʼ conditional distribution. An RBM consists of a visible layer (a layer of p visible Omics features) and a hidden layer (a layer of g hidden variable). The Gaussian RBM model assumes that the conditional distributions of visible variables (Omics features) given hidden variables follow a Gaussian distribution.

Illustrative Case-studies: Lemon-Tree has been applied to TCGA glioblastoma expression and copy-number data (Bonnet et al., 2015). It resulted in a module network composed of 121 clusters of co-expressed genes and a list of prioritized (high-scored) genes, mostly associated with amplified/deleted regions. Several of these high-scored genes were already reported as cancer genes in glioblastoma (including EGFR, PDGFRA, FGFR3, PIK3CA, MDM4, CDKN2A/B, and PTEN) where all involved in glioblastoma driver pathways, including proliferation, apoptosis, and angiogenesis pathways. Besides the well-known genes, Lemon-Tree could also identify a few novel markers that have rarely or never studied glioblastoma. For instance, INSR was involved in several modules. It stimulates cell proliferation and is aberrantly expressed in cancer cells (Belfiore et al., 2009); therefore, amplification of INSR in glioblastoma may enhance proliferation. PAOX (polyamine oxidase) is another novel marker that might have tumor suppressor activity via amine oxidase activity and their primary involvement in cancer growth inhibition and progression (Guzeloglu-Kayisli et al., 2004). Interestingly, PAOX was biologically relevant based on its prognostic value via a survival analysis.