Commonly shared features (genes/probes) and samples are chosen across all the profiles from the multi-omics cancer dataset. Specifically, probes (features) from DNA methylation arrays containing any missing values are discarded. The individual profile is normalized using the zero-mean normalization for each feature (Bandyopadhyay et al., 2013), as described in the following formula: x i k ′ = x i k − μ σ . Here, μ is the mean across the data for the feature i prior to normalization, and σ denotes standard deviation. x _ik and x i k ′ signify the value of the i-th feature at k-th patient (sample) prior and after normalization, respectively. To determine statistically significant features, the empirical Bayes statistical test is applied using the package “Linear Models for Microarray and RNA-Seq Data” (Smyth, 2004; Bandyopadhyay et al., 2013), which works better on the dataset with a small sample size. The moderated t-statistic (Ritchie et al., 2015) is elaborated as follows: t ̃ p r = 1 1 m 1 + 1 m 2 β ̂ p r s ̃ p r , (1) where m ₁ and m ₂ are the number of patients (cases) and that of the normal samples (controls), respectively. Here, β ̂ p r signifies the contrast estimator for the feature pr, whereas s ̃ p r 2 refers to the posterior sample variance for pr. The statistic to compute the contrast estimator for the probe pr is formulated as follows: β ̂ p r | σ p r 2 ∼ N ( β p r , σ p r 2 ) . Here, N represents the normal distribution. The statistic to estimate the posterior sample variance for pr is formulated as follows: s ̃ p r 2 = d 0 s 0 2 + d p r s p r 2 d 0 + d p r , (2) where d ₀ ( < ∞ ) signifies the prior degrees of freedom, and s 0 2 denotes the variance. In addition, d _pr ( > 0 ) symbolizes the experimental degrees of freedom of pr, and s p r 2 denotes the sample variance of pr. The significance of the level of the p-value is then determined from s ̃ p r 2 with the help of the cumulative distribution function (cdf). If the p-value of the feature is less than the standard cutoff of 0.05, the feature is defined as statistically significant. The filtered differentially expressed features are then ordered according to the p-values. Notably, if any gene corresponds to more than one probe (feature), the probe with the lowest p-value will be selected to represent the gene, and the rest of the probes for the gene are simply ignored. We apply the same approach to each layer of the molecular profile, and then, we perform the combination of the significant non-redundant features (genes/probes/copy number variation, etc.) from all layers (let, UF).

Let o _i and o _j denote two object types, namely, gene expression and DNA methylation, in all resulted features UF. The number of genes is N, while each gene is denoted by n _i , where i = 1, 2,…, N. There are M number of DNA methylation samples, while each sample is termed as m _j , where j = 1, 2, …, M. In addition, there is a P set consisting of p types of profiles from the multi-omics datasets. The input to this implemented variant of the 3-FPNMF model is R , which is a relational block matrix shown as follows: R = ∗ R 12 … R 1 p R 21 ∗ … R 2 p ⋮ ⋮ ⋱ ⋮ R p 1 R p 2 … ∗ . (3) Here, ∗ denotes that similar object relationships are not considered in this approach. R _ij denotes the relationship between o _i th and o _j th object types. The respective correlation of the xth object of type o _i (e.g., gene) and the yth object of type o _j (e.g., sample) is represented as R o i o j ( x , y ) . In this implementation, we have experimented with six object types, as described later.

For each object type from each profile, there is a constraint in the input constraint block diagonal matrix, as shown in the following expression: τ P = D i a g τ 1 , τ 2 , … , τ p . (4)

The relational block matrix R is tri-factorized into matrix factors G and S (Žitnik and Zupan, 2014), which is shown as follows: G = D i a g G n 1 × m 1 1 , G n 2 × m 2 2 , … , G n p × m p p , (5) S = * S 12 r 1 × r 2 … S 1 p r 1 × r p S 21 r 2 × r 1 * … S 2 p r 2 × r p ⋮ ⋮ ⋱ ⋮ S p 1 r p × r 1 S p 2 r p × r 2 … * . (6) Here, r denotes rank factorization to the object type o _p inferenced by the 3-FPNMF model. The factor S denotes the block relation between object types o _i and o _j . The factor G o i reconstructs relations specifically to the object type o _i .

Thus, each relation matrix R o i o j obtains matrix factorization as G o i S o i o j G o j T . In a simplified way, this relational block 3-FPNMF model is shown as follows: * G o 1 S o 1 o 2 G o 2 T … G o 1 S o 1 o p G o p T G o 2 S o 2 o 1 G o 1 T * … G o 2 S o 2 o p G o p T ⋮ ⋮ ⋱ ⋮ G o p S o p o 1 G o 1 T G o p S o p o 2 G o 2 T … * . (7)

The objective function of this tri-factor penalized matrix decomposition (PMD) model is to minimize the distance between the input block relational matrix R and its 3-FPNMF system adhering to the constraint matrix τ ^P , which is shown as follows: min G ≥ 0 j R : G , S = ∑ R o i o j ∈ R ‖ R o i o j − G o i S o i o j G o j T ‖ 2 + ∑ p = 1 P t r G T τ p G . (8)

Here, ‖.‖ denotes the Frobenius norm, and tr (.) denotes the trace. Our sparse implementation for this 3-FPNMF model reduces the missing relational matrix problem with zero values. Our model is more suitable for real-life heterogeneous datasets with missing values, which differs from those of Žitnik and Zupan (2014) in its non-negative sparse implementation. Our proposed 3FPNMF − MKL model is shown briefly in Figure 1, while a detailed flowchart is represented in Supplementary Figure S1.

Next, we introduce the multiple Kernel Learning (MKL) algorithm (Xu et al., 2013) with the hinge loss soft margin, in which the classifier and the kernel combination coefficients are optimized by solving the hinge loss soft margin MKL problem.

After using the 3-FPNMF model in the first phase, the approximate sparse relation matrix R ̂ o i o j for target object type pairs o _i and o _j is reconstructed as R ̂ o i o j = G o i S o i o j G o j T . (9) Then, to develop kernel fusion, the resulting kernel matrices are generated using the “Kernel Trick”: K ( o i , o j ) = R ̂ o i o j . R ̂ o i o j T . The kernels are further normalized and smoothed using 2-dimensional linear filters.

Given p base-kernels K = { K 1 , K 2 , … , K p } developed from the reconstructed relational block matrix R ̂ = { R ̂ o i o j | i = 1 , … , p ; j = 1 , … , p } , kernel slack variables for the kernel K p ∈ K are defined as the difference between the target margin θ and the SVM dual objective function D S V M K p , α = max α ∈ R N ∑ n = 1 N α n − 1 2 ∑ m = 1 N ∑ n = 1 N α n α m y n y m K p x n , x m subject to ∑ n = 1 N y n α n = 0 , α n ≥ 0 , ∀ n . Then, the slack variable is ζ _p = θ − DSVM(K _p , α), and the hinge loss is shown as follows: z p = ℓ ζ p = max 0 , ζ p . (10) Therefore, the objective function for this hinge loss soft margin MKL algorithm becomes min θ , α ∈ D o m α , ζ p − θ + π ∑ p = 1 P ζ p . (11) subject to DSVM(K _p , α) ≥θ − ζ _p , ζ _p ≥ 0, p = 1, …, P.

The objective of the aforementioned hinge loss soft margin MKL is to maximize the margin θ while considering the “errors” from the given P-based kernels. The parameter π balances the contribution of the loss term represented by slack variables ζ _p and the margin θ. π should be in the range {π|π ≥ 1/P}. Otherwise, there is no solution to the proposed problem. Our proposed framework for gene signature detection from heterogeneous data sources using the 3FPNMF − MKL model is depicted in Figure 2.

In biological datasets such as TCGA, the clinical data are made available. This includes patient sample groups, biological subtypes, drug treatment, and survival/prognosis information. In our current study, we obtain accuracies for different combinations of class labels using the non-matrix factorization technique for the case where there were more than two class labels or subtypes. Among them, the combination of class labels, which produces the highest area under curve (AUC), is chosen for the next step (i.e., module detection). Say, q is the specific combination of class labels, which produces the highest AUC. Find q = {∃i, ∃j}|{∃a, ∃b, ∃k} such that A U C q = arg max ∀ i , j A U C c l i c l j ′ , ∀ a , b , k A U C c l a c l b k ′ , (12) where cl denotes the left part of the group combination, cl′ signifies the right part of any sample group combination, and i ∈ {1, 2, …, (m − 1)}, j ∈ {(i + 1), (i + 2), …, m}, a ∈ {1, 2, …, m}, b ∈ {1, 2, …, m} & b ≠ a, k ∈ {, 2, …, m}, and k ≠ a and k ≠ b.

After selecting the right class-label combination, we extracted the sub-gene expression data consisting of only the selected class labels and then used them for gene module detection and signature identification.

In our procedure, we first evaluated the power of the soft thresholding, which was then applied to evaluate the adjacency matrix using Pearson’s correlation. The topological overlap matrix (TOM) similarity score (Ravasz et al., 2002) was computed from the employed adjacency matrix. The TOM score between two nodes (say, i and j) symbolized as TOM(i, j) is defined as follows: T O M i , j = ∑ v ≠ i , j X i , v X j , v + X i , j min ∑ v ≠ i X i , v , ∑ v ≠ j X j , v − X i , j + 1 , if i ≠ j , 1 , if i = j , (13) where X denotes the corresponding adjacency matrix containing Boolean entries. The entry of 1 indicates that the corresponding two nodes share the same connection (i.e., direct connection), while the entry of 0 signifies that no direct connection exists between them.

After obtaining the TOM score, we computed the distance/dissimilarity value between genes (i and j) denoted by dissTOM(i, j),which is shown as follows: dissTOM(i, j) = 1 − TOM(i, j). We conducted average linkage clustering on the multi-omics dissimilarity matrix dissTOM via considering all potential pairs of genes/features. Finally, the dynamic tree cut technique (Langfelder et al., 2008) was applied on the clustering dendrogram to determine the gene modules. In order to evaluate the quality of the aforementioned clustering, we calculated different cluster validity index measures, viz., cluster coefficient, heterogeneity, Dunn Index, maximum adjacency ratio, centralization, silhouette width, and scaled connectivity.

After finding the gene modules, we estimated Pearson’s correlation coefficient (PCC) between each gene pair within the resulted modules. For each module, the mean of the correlations for each gene pair within that particular module was obtained. The module with the maximum mean correlation coefficient was elected as a gene signature. Notably, genes with the elected gene signature are differentially expressed between case and control samples. In order to validate the classification performance of the employed gene signature, we utilized the Prediction Analysis of Microarrays (PAM) classifier with 10-fold cross-validation (CV) on the expression sub-data to classify the underlying class labels. The entire procedure was then repeated ten times. Moreover, we calculated the average scores of several classification performance metrics such as sensitivity, specificity, precision, accuracy, and AUC, individually.

We carried out KEGG pathway and Gene Ontology (GO) analyses using the Enrichr database (Chen et al., 2013). Notably, GO terms can be categorized into three kinds, viz., biological process (BP), cellular component (CC), and molecular function (MF). Those significant pathways/GO terms with an adjusted p-value less than 0.05 were identified. Meanwhile, literature research studies were also performed to identify disease-related pathways/GO terms.

For our experiment, TCGA acute myeloid leukemia (LAML) multi-omics dataset (https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Acute%20Myeloid%20Leukemia%20(LAML)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) contained six heterogeneous profiles such as the gene expression (IlluminaGA) profile, DNA methylation (Illumina Methylation 27k) profile, exon expression (IlluminaGA) profile, miRNA profile, pathway activity (Paradigm IPLs) profile, and copy number (GISTIC2) profile. Initially, the gene expression profile included 179 samples and 20,113 genes. For the methylation profile, there are 194 samples and 27,578 methylation probes. Particularly, for the methylation profile, many genes are profiled with more than one probe. In the exon expression profile, there are a total of 219,296 chromosome ids and 179 samples. Here, many genes are connected with more than one chromosome id. The miRNA profile contains 705 miRNAs and 188 samples. The pathway activity profile has 7,203 genes and 173 samples, while the copy number profile consists of 24,776 genes and 191 samples. There are three categories of samples (i.e., class labels) for the LAML multi-omics dataset: i) favorable, ii) intermediate (also called normal), and iii) poor. Specifically, every profile consists of 161 commonly shared LAML samples. Among them, 31 samples belong to the first category, 96 samples are in the second category, and the rest of the samples (= 34) are in the third category. In addition, there are 1,501 uniquely matched genes among the five profiles [i.e., gene expression, DNA methylation, exon expression, pathway activity, and copy number variation (GISTIC2) profiles].

First, we selected the sub-data, which contain commonly shared samples (i.e., 161) and genes (i.e., 1,501) for each of the five profiles (i.e., gene expression, DNA methylation, exon expression, pathway activity, and copy number variation profiles). Many matched genes are connected with more than one probe (or chromosome id) for each profile. In the case of the miRNA profile, we started to work with the matched samples (n = 161) and all of its miRNAs (n = 705). The empirical Bayes test is performed by limma software on each gene probe or chromosome id for each of the five profiles (i.e., gene expression, DNA methylation, exon expression, pathway activity, and copy number variation profiles) across all the three classes (viz., favorable, intermediate, and poor).

Notably, since there are three classes/groups of samples, here, limma is initially performed between each group pair (i.e., i) favorable vs. intermediate, ii) intermediate vs. poor, and finally iii) favorable vs. poor), then an F-statistics is computed, and finally, the respective p-value is generated from the F-statistics. After the test, for every gene, we only selected the probe or chromosome id with the lowest p-value achieved among all the probes or chromosome ids connected with that gene. As a result, we obtained 728, 272, 1,100, 265, and 904 significant genes for the gene expression, methylation, exon expression, pathway activity, and copy number profiles, respectively. Thereafter, we took the combination of all the significant gene sets, which led to a molecular set of a total of 1,388 genes. Furthermore, the same significance test was applied on each miRNA of the miRNA profile across all the three classes (viz., favorable, intermediate, and poor) as well. We obtained a total of 229 significant miRNAs.

Using the non-matrix factorization technique, we obtained accuracies for different combinations of class labels such as i) Class 1 (favorite) vs. Class 2 (intermediate), ii) Class 1 vs. Class 3 (Poor), iii) Class 1 vs. classes 2 and 3, iv) Class 2 vs. Class 3, v) Class 2 vs. classes 1 and 3, and vi) Class 3 vs. Classes 1 and 2 (as depicted in Table 1). Among them, the second combination, i.e., Class 1 vs. Class 3 produced the highest area under curve (AUC = 0.7713). Hence, we selected the combination for gene signature discovery since other combinations did not produce better AUC scores. After obtaining right combinations of class labels, we first evaluated the power (=1) for soft thresholding (illustrated in Figure 3A), which was then applied to estimate the adjacency matrix through Pearson’s correlation score. Then, the TOM score and distance matrix were computed. To determine gene modules, we applied average linkage clustering and dynamic tree cut methodologies. As a result, we generated a total of 10 gene modules. The numbers of participating differentially expressed genes (DEGs) for these 10 gene modules (represented by black, blue, brown, green, magenta, pink, purple, red, turquoise, and yellow colors) were 50, 99, 90, 74, 23, 25, 22, 51, 214, and 80, respectively. The dendrogram is represented in Figure 3B. The corresponding cluster validity indices in that module detection are illustrated in Table 2. The Average silhouette width plot generated during clustering is represented in Supplementary Figure S2. PCC was calculated between each gene pair within each module. The mean correlation scores of the 10 modules (depicted by blue, green, turquoise, magenta, brown, red, yellow, black, purple, and pink colors) were 0.0268, 0.2562, 0.0321, 0.3914, 0.1143, 0.0215, 0.0570, 0.4029, 0.3455, and 0.1605, respectively. The black module had the highest mean correlation coefficient score (= 0.4029 in Table 3). Thus, it was selected as the gene signature. The resultant gene signature contained 50 DEGs (see Table 3). To verify the classification performance of the resultant signature, we applied the PAM classifier through the 10-fold cross-validation (CV) on all the features and samples of signature data in order to classify the groups (favorite and poor). The entire procedure was then repeated 10 times. In the experiment, the mean sensitivity, mean specificity, mean precision, mean accuracy, and mean AUC were 69.12%, 84.19%, 82.79%, 76.31, and 0.8273, respectively (see Figure 4; Table 4). Based on the gene set enrichment analysis on the 50 genes of the signature using the Enrichr web database, we extracted significant KEGG pathway and Gene Ontology (GO) terms. Among the KEGG pathways, the Rap1 signaling pathway (hsa04015) is the most significant pathway (adjusted p-value = 7.497 × 10⁻⁰⁶) that contains eight genes (viz., EFNA1, GNAO1, TIAM1, CSF1, ITGB3, ITGA2B, THBS1, and MAPK13). Second, the most significant pathway is the PI3K-Akt signaling pathway (hsa04151) with an adjusted p-value of 1.128 × 10⁻⁰⁵, which consists of nine genes (viz., EFNA1, CSF1, ITGB3, ITGA2B, IL2RB, FASLG, TP53, THBS1, and EPOR). The following eight pathways are the cytokine–cytokine receptor interaction (hsa04060) (adj. p-value = 1.437 × 10⁻⁰⁵), inflammatory bowel disease (IBD) (hsa05321) (adj. p-value = 2.1E-05), proteoglycans in cancer (hsa05205) (adj. p-value = 2.1 × 10⁻⁰⁵), hematopoietic cell lineage (hsa04640) (adj. p-value = 6.752 × 10⁻⁰⁵), T-cell receptor signaling pathway (hsa04660) (adj. p-value = 1 × 10⁻⁴), TNF signaling pathway (hsa04668) (adj. p-value = 2 × 10⁻⁴), osteoclast differentiation (hsa04380) (adj. p-value = 3 × 10⁻⁴), and Ras signaling pathway (hsa04014) (adj. p-value = 3 × 10⁻⁴) (also see Table 5). Among the significant GO:BP terms, the positive regulation of cellular metabolic processes (GO:0031325) (adjusted p-value = 8.02947 × 10⁻⁰⁵) was ranked as the most significant, which contains six genes (EDN1, CSF1, CCL5, GATA3, THBS1, and TP53). The second most significant GO term is the regulation of inflammatory responses (GO:0050727) with an adjusted p-value of 8.029 × 10⁻⁰⁵. This term consists of seven genes (CCL5, CCL4, RORA, GATA3, ETS1, BIRC3, and MAPK13) (Table 5). Among the significant GO:CC terms, the platelet alpha granule (GO:0031091) (adjusted p-value = 4 × 10⁻³) contains four genes (viz., ITGB3, ITGA2B, A2M, and THBS1), while among the GO:MF terms, the core promoter binding factor (GO:0001047) (adjusted p-value = 8 × 10⁻⁴) contains five genes (viz., RORA, GATA3, GATA1, TP53, and ARNTL). For details of the top significant pathways and GO terms, see Table 5.