It is critical to learn effective feature representation in clustering tasks. Compared with traditional methods, deep learning methods can extract more advanced data feature representations and are widely used in various fields. In order to extract the high-level feature representation of omics data, we use the Stacked Autoencoder (SAE) model with the strongest generalization performance to learn the original omics data. The training process of SAE model can be found in Figure 2 (See SAE Model).

Suppose there are l layers in the SAE. In the encoder stage, when SAE is used to learn omics data X v , the learning of the l -th layer is written as Z v l Z v l = ϕ W e v l Z v l − 1 + b v l e (3) In Equation 3 ϕ is the activation function of the full connection layer. Here we use the LeakyRELU activation function. W v l e and b v l e are the weight matrix and bias of the l -th layer in the encoder, respectively. When the encoder starts learning, the feature representation is initialized as: Z v 0 = X v .

In the decoder stage, the input data is reconstructed through multiple fully connected layers, which can be written as Z v l = ϕ W v l d Z v l − 1 + b v l d (4) In Equation 4 W v l d and b v l d are the parameters of l -th layer in the decoder.

The final output Z v l is the output X ^ v of SAE: X ^ v = Z v l . We hope that X ^ v can reconstruct the original omics data X v as much as possible, and then use the following loss function in Equation 5 for SAE model training L r e s = 1 2 N ∑ v V X ^ v − X v F 2 (5)

SAE can learn the advanced feature representation of omics data, but it does not consider the structural information among omics data samples. We introduce Graph Convolutional Network (GCN) to learn the structural representation of each omics data. The training process of GCN model can be found in Figure 2 (See GCN Model).

For omics data X v , GCN learns the structural representation G v l of the l -th layer through the following convolution operations G v l = C A ^ v G v l − 1 W v l − 1 (6) where W v l − 1 is the weight matrix of l − 1 -th layer. A ^ v = A v + I , where I is an identity matrix. According to Equation 6, GCN can learn the representation G v l of the next layer through G v l − 1 , W v l − 1 and the adjacency matrix A ^ v .

The information fusion learning process of single omics data by combining SAE and GCN model can be found in Figure 2. Considering both Z v l − 1 and G v l − 1 , we can obtain a joint representation H v l − 1 with more effective information through the following formula H v l − 1 = 1 − ε G v l − 1 + ε Z v l − 1 (7) where ε is the balance parameter used to balance the relationship between the two representations Z v l − 1 and G v l − 1 . For simplicity, we set it to 0.5. Through Equation 7, we have realized the connection between SAE and GCN network. And H v l − 1 contains both feature representation information and structure representation information.

Next, we need to learn the l -th layer representation G v l of GCN. At this time, H v l − 1 is taken as the input of GCN. Then we have G v l = C A ^ v H v l − 1 W v l − 1 (8)

In the traditional GCN model, after the multi-layer graph convolution operation is adopted, the characteristics of different nodes tend to be homogenized, that is, the characteristics of all nodes within the same connected component are almost the same. This is the so-called over-smoothing phenomenon. The representation information learned by the SAE in each layer is very different, and in Equation 8, the joint representation H v l − 1 contains both the feature information learned and the structured information learned. Therefore, the existence of Equation 7 can alleviate the over-smoothing problem of GCN.

It is worth noting that the input data matrix G v 1 of the first layer can be calculated by using omics data X v . G v 1 can be defined by Equation 9 G v 1 = C A ^ v X v W v 1 (9)

The final output of GCN is determined according to Equation 10 G v l = C A ^ v H v l − 1 W v l − 1 (10)

To demonstrate the effectiveness of SMMSN, we applied it to clustering tasks on three labeled real-world cancer multi-omics datasets. These datasets include the ROSMAP dataset for Alzheimer’s disease (AD) patients and normal control (NC) classification, the Low Grade Glioma (LGG) dataset for Grade 2 and Grade 3 classification in low-grade glioma, and the Pan Kidney Cohort (KIPAN) dataset for the classification of three kidney cancer types: Chromophobe Renal Cell Carcinoma (KICH), Clear Renal Cell Carcinoma (KIRC), and Papillary Renal Cell Carcinoma (KIRP) (Wang et al., 2021). The ROSMAP dataset is composed of ROS and MAP, both of which are longitudinal clinical-pathologic cohort studies of AD from Rush University (Bennett et al., 2012; De Jager et al., 2018). It is available through the AMP-AD Knowledge Portal (https://adknowledgeportal.synapse.org/) (Hodes and Buckholtz, 2016). The omics data for LGG and KIPAN were obtained from TCGA via Broad GDAC Firehose (https://gdac.broadinstitute.org/). For each dataset, we used three types of omics data (i.e., mRNA expression data, DNA methylation data, and miRNA expression data) for clustering to provide comprehensive and complementary information about the diseases. Only samples with matched omics data were included for each data type. Below are detailed descriptions of the datasets.

• ROSMAP: 55889 genes for mRNA expression, 23788 genes for DNA methylation, 309 genes for miRNA expression, 351 patients (NC: 169 patients, AD: 182 patients).

To further validate the efficacy of SMMSN for cancer subtype discovery, it is used to process multiple omics data sourced from TCGA, as preprocessed by Wang et al. (2014). Our study encompassed five distinct cancer types: Breast Invasive Carcinoma (BIC), Glioblastoma Multiforme (GBM), Lung Squamous Cell Carcinoma (LSCC), Kidney Renal Clear Cell Carcinoma (KRCCC), and Colon Adenocarcinoma (COAD). For each cancer type, we analyzed three types of omics data obtained from different platforms: mRNA expression, DNA methylation, and miRNA expression. Detailed descriptions of these multi-omics datasets for the five cancer types are provided below.

• GBM: 12,042 genes for mRNA expression, 1,305 genes for DNA methylation, 534 genes for miRNA expression, 213 patients.

For labeled cancer multi-omics data, we used the Accuracy (ACC) for evaluation to validate the clustering results. ACC quantifies the consistency between the clustering results and the true labels, and its calculation formula is defined as Equation 19: ACC = ∑ i = 1 n δ t i , map y i n (19) where t i is the true label, y i is the label assigned by the clustering methods, δ ⋅ is the indicator function, which equals 1 when t i = map y i , and 0 otherwise. map ⋅ represents an optimal mapping function that best matches the clustering labels to the true labels. By calculating ACC, we can intuitively evaluate the clustering performance.

For unlabeled datasets, this study performs survival analysis on cancer subtypes identified through clustering to assess survival disparities among sample groups derived from the proposed algorithm. In statistical analysis, hypothesis testing, such as the Cox Log-rank Test (CLT) (Hosmer et al., 2000), is employed to quantify differences in survival curves. CLT is a non-parametric method commonly used to evaluate whether variations in survival between subtypes are significant. A lower p-value from this test suggests stronger evidence against the null hypothesis, indicating substantial differences in survival outcomes that are unlikely to be due to chance alone. Additionally, the Kaplan-Meier estimation method (Hosmer et al., 2000) is utilized to derive survival functions and construct Kaplan-Meier survival curves. These curves plot the survival rate on the y-axis against time from the start of observation to the last recorded time point on the x-axis. They visually illustrate how the event (e.g., survival or recurrence) unfolds over time for different cancer subtypes, providing insight into their respective prognostic outcomes.

For comparison purposes, we included five established traditional multi-view clustering algorithms known for their efficacy in cancer subtype prediction: PFA (Shi et al., 2017), SNF (Wang et al., 2014), ANF (Ma and Zhang, 2017), and MVCMO (Yu et al., 2019). Two deep learning-based cancer subtypes discovering methods, Subtype-Former (Yang et al., 2022) and Subtype-DCC (Zhao et al., 2023), are also taken as the competing methods. These algorithms are widely recognized in the field for their ability to integrate diverse data sources and identify meaningful subtypes within cancer datasets.

The deep learning algorithms involved in this study were implemented using the popular deep learning framework PyTorch 3.9, and all experiments were conducted on an NVIDIA GeForce RTX 4080 GPU with 32 GB RAM, Core I7-12700K. To evaluate the performance of the models and comparison methods, each experiment was run five times, and the average accuracy score along with the standard deviation was reported to ensure the robustness and comparability of the results. The parameter settings for the deep learning models are as follows.

• For the SMMSN algorithm, the network output dimension was set to 100, and the adjustable adjacency matrix parameter k was defined as 40. The Adam optimizer was used during training, with an initial learning rate of 1 × 10⁻⁴ and a decay factor of 1 × 10⁻¹⁵. The model was trained for 500 epochs.

For the benchmark machine learning algorithms, they were implemented by MATLAB 2022a software, and their parameters were set strictly according to the guidelines provided by the authors. Each experiment was run five times, and the average accuracy score along with the standard deviation was reported to ensure the robustness and comparability of the results The specific parameters are set as follows.

• For the PFA algorithm, the local sample-spectrum for each biological data type was captured using the Algorithm_1 function from the PFA package. Next, the global sample-spectrum was captured using the Algorithm_4 function, with the hyperparameter lambda set to 1.

For the labeled cancer multi-omics data, the number of clusters corresponds to the number of cancer subtypes in the data itself. The number of clusters for the three datasets, KIPAN, ROSMAP, and LGG, is set to 3, 2, and 2, respectively. For the unlabeled data, we follow the commonly accepted number of cancer subtypes as reported in the majority of studies, such as in references (Wang et al., 2014; Ma and Zhang, 2017; Yu et al., 2019). The number of clusters for the five datasets, GBM, BIC, KRCCC, LSCC, and COAD, is set to 3, 5, 3, 4, and 3, respectively.

Table 2 shows the p-values from survival analysis between SMMSN and competing methods across five datasets. This comparison evaluates the statistical significance of survival differences among cancer subtypes identified by each algorithm. Across all five cancer types, SMMSN consistently yielded the lowest p-values compared to other algorithms. Figure 4 presents Kaplan-Meier survival curves generated by SMMSN for each cancer type, depicting the survival trends of respective subtypes. Each curve in Figure 4 illustrates the survival times of distinct cancer subtypes, with sample counts annotated for clarity. These results demonstrate SMMSN’s ability to discern significantly distinct cancer subtypes across various cancer types.

To further validate the effectiveness of each module in SMMSN, we conducted ablation studies. The SMMSN algorithm mainly consists of three key components: the Feature Representation Fusion (FRF) module based on SAE, the Structural Information Fusion (SIF) module based on GCN, and the Dual Self-supervised (DSS) module. The results of the ablation study are shown in Table 3. It is important to note that when we use only the FRF module or the SIF module, only a single self-supervised learning operation is required, which is denoted as SS module in Table 3. In other words, when both the FRF and SIF modules are used simultaneously in SMMSN, we apply the dual self-supervised module for model learning. From the ablation results shown in Table 3, we can observe different performance outcomes for three different module combinations across five unlabeled multi-omics datasets (GBM, BIC, KRCCC, LSCC, COAD). The analysis can be broken down as follows.

• When only the GCN module and single self-supervised module are used, the results are relatively poor across all datasets, particularly on the KRCCC and LSCC datasets, with p-values of 9.60E-2 and 2.26E-2, respectively. This suggests that while the GCN module can capture structural features, its performance is limited without the feature representation fusion from the SAE module.

GBM stands as the most prevalent and deadly primary brain tumor in adults, categorized within the glioma group. Numerous studies have extensively explored GBM at the molecular level, identifying distinct cancer subtypes with corresponding clinical implications. For instance, Verhaak et al. (2010) classified GBM based on mRNA expression into four subtypes: Mesenchymal, Classical, Neural, and Proneural. Another study (Noushmehr et al., 2010) differentiated GBM into G-CIMP and non-G-CIMP subtypes based on CpG Island Methylator Phenotype (CIMP).

Using GBM data, we analyzed the distribution of clustering results obtained by SMMSN across the subtypes identified in the aforementioned studies, summarized in Table 4. It is worth noting that the cancer subtypes in references (Verhaak et al., 2010) and (Noushmehr et al., 2010) are classifications derived from different research methods and standards, but they are not considered gold standards for cancer subtypes. Instead, they serve as reference classifications used to help understand and validate the biological differences between the three subtypes identified by the SMMSN algorithm. Table 4 highlights that a majority of patients in subtype 1 align with the Proneural subtype. Subtype 2 shows a closer association with Classical and Proneural subtypes. Subtype 3 predominantly corresponds to the Mesenchymal subtype. It shows that the three subtypes identified have certain differences. Notably, all patients in subtypes 2 and 3 belong to the non-G-CIMP category, while a portion of patients in subtype 1 are classified under G-CIMP. This indicates The difference between the identified subtype 1 and subtype 2–3 (subtype 2 and subtype 3) was obvious, and this conclusion was also verified in Figure 4A that subtype 1 had a longer survival time than subtype 2–3.

Subsequently, we further compared long-term survival subtype 1 with short-term survival subtype 2–3 and looked for their differences in gene mutations. Figure 5 show the difference in Copy Number Variation (CNV) abundance between long-lived subtype 1 and short-lived subtype 2–3. In this figure, each point represents a gene, and its axis is the number of patients carrying a variant of that gene in the two survival differential subtypes. The most abundant mutated genes in subtypes 2–3 were EGFR, SEC61G and RP11-745C15.2. EGFR mutations and amplifications are very common in GBM, especially the EGFRvIII variant, which drives rapid tumor cell proliferation, increased invasiveness, and resistance to treatment. EGFR overexpression is closely associated with the progression of more malignant subtypes, which generally indicate poorer survival outcomes (Hu et al., 2022). As a key component of the SEC61 translocation complex in the endoplasmic reticulum, SEC61G is involved in regulating protein transport and processing. Abnormalities in SEC61G may affect proteins involved in cell proliferation and stress responses, thus promoting tumor growth and progression (Zeng et al., 2023). RP11-745C15.2 represents a class of long non-coding RNAs (lncRNAs), whose role in cancer is becoming increasingly recognized. RP11-745C15.2 may regulate key oncogenes like EGFR or its downstream signaling pathways, enhancing malignant cell behavior. In GBM, several lncRNAs, including RP11 family members, are believed to be involved in tumor progression by regulating gene expression, influencing cell growth, and contributing to the aggressive nature of the tumor (Zhang et al., 2020). To sum up, EGFR mutations drive malignant cell proliferation, while SEC61G and RP11-745C15.2, through their roles in protein transport and gene regulation, further promote tumor cell growth, survival, and invasiveness. Their combined action leads to greater gene variation in subtypes 2 and 3, resulting in higher aggressiveness and worse prognosis. This genetic association suggests that these gene alterations are key drivers of tumor progression in the short-term survival subtypes, helping to distinguish subtype 1 (long-term survival) from subtypes 2 and 3 (short-lived).

Further analysis of the cancer subtypes identified by SMMSN involved accessing clinical data for all GBM patients from the cBio Cancer Genomics Portal database. Figure 6 illustrates boxplots depicting the distribution of survival time and age among these subtypes, demonstrating discernible differences. In Figure 6A, subtype 1 exhibits significantly longer survival compared to subtype 2 and subtype 3, supported by p-values from two-sided Welch’s t-tests: 1.45E-4 and 5.69E-3, respectively. Figure 6B reveals that patients in subtype 1 are younger than those in subtype 2 and subtype 3, with corresponding t-test p-values of 2.37E-7 and 5.21E-5, respectively. Moreover, we conducted an Analysis of Variance (ANOVA) test across the three subtypes, confirming significant differences in both survival time (p = 8.24E-7) and age distribution (p = 2.48E-9). Similarly, the Kruskal–Wallis test also indicated statistically significant distinctions in age (p = 1.97E-7) and survival time (p = 1.84E-04) among the subtypes. These consistent findings underscore the biological relevance and statistical significance of the identified subtypes with respect to both age demographics and survival outcomes.

Figure 7 presents Kaplan-Meier survival curves depicting the response of patients to the drug Temozolomide (TMZ). Patients are stratified into two groups: those treated with TMZ and those not treated with TMZ. The p-values for subtype 1, subtype 2, and subtype 3 are 0.65, 4.12E-5, and 5.42E-2, respectively. These results indicate that TMZ treatment has minimal impact on the survival outcomes of patients in subtype 1, while it significantly affects the survival of patients in subtypes 2 and 3.

Differential gene expression and GO enrichment analyses were conducted on GBM data to assess differences among the three subtypes identified by SMMSN. Initially, significant differentially expressed genes across the subtypes were identified using the ANOVA method. Figure 8 displays a heatmap illustrating the top 1,000 differentially expressed genes in mRNA expression data, with panels A, B, and C representing subtype 1, subtype 2, and subtype 3, respectively. The heatmap reveals distinct clusters among the differentially expressed genes, indicating subtype-specific expression patterns.

Further analysis involved functional enrichment of these differentially expressed genes. Figure 9 presents the results of GO enrichment analysis categorizing the genes into four distinct groups (X1, X2, X3, and X4). Each group is associated with specific GO biological processes, as indicated by the number of enriched genes listed below. Notably, genes related to the regulation of mRNA metabolism/chromosome organization were downregulated in subtype 3 and upregulated in subtype 2. The genes related to the regulation of immune effector process/leukocyte-mediated immunity/lymphocyte-mediated immunity and hetero-cell bonding were downregulated in subtype 2. Genes associated with DNA recombination, nuclear transport/export, and RNA splicing were downregulated in subtype 3. In conclusion, GBM subtypes identified by SMMSN have obvious differences in clinical indicators and molecular levels.

BIC refers to a malignancy in which cancer cells have penetrated the basement membrane of breast ducts or lobular acinus and invaded the stroma. Similar to the GBM case analysis procedure described above, we first made a comparison with previous BIC subtype study. PAM50 is a molecular subtype of BIC based on quantitative detection of the expression levels of 50 functional genes in breast tumor tissues, including Luminal A, Luminal B, HER2-enriched and Basal-like subtypes (Parker et al., 2009). Table 5 depicts the comparison results of the distribution of BIC subtypes identified by SMMSN in PAM50 subtypes.

From Table 5, it can be observed that the five cancer subtypes identified by the SMMSN algorithm show different distribution patterns in the PAM50 subtypes (Luminal A, Luminal B, Basal-like, HER2-enriched) from the study in (Parker et al., 2009). SMMSN subtypes 1, 2, 3, and 4 are mainly concentrated in Luminal A, with subtype 2 almost entirely composed of Luminal A patients, indicating a high level of consistency between these subtypes and the Luminal A subtype. SMMSN subtype 5, on the other hand, is predominantly composed of Basal-like patients, suggesting a strong correspondence with the Basal-like subtype. In contrast, Luminal B and HER2-enriched patients are more dispersed across multiple SMMSN subtypes, especially in subtypes 1 and 4, revealing a certain degree of discrepancy between the subtypes identified by SMMSN and the PAM50 subtypes. These observations reflect a strong alignment between SMMSN subtypes and PAM50 in certain subtypes, while in others, cross-subtype distribution patterns are apparent.

To further validate the biological differences among the cancer subtypes 1, 2, 3, 4, and 5 identified by SMMSN, especially the first four subtypes, differential expression analysis of BIC expression data was performed using the Kruskal–Wallis test. Hierarchical clustering method was used to group the top 1,000 differential genes of BIC, and GO: BP analysis was performed based on the grouping results, and their results are presented in Figures 10, 11. Figure 10 presents the heatmap of the top 1,000 genes whose mRNA expression varies significantly across the five subtypes identified by SMMSN in BIC. It can be observed that there exists distinct differences in gene expression among the five subtypes. Figure 11 shows the functional enrichment analysis on the differentially expressed genes identified in Figure 10. It can be found that genes related to Wnt signaling pathway and epidermal development are upregulated in subtype 5 and downregulated in subtype 2 and 4. Genes related to extracellular matrix/structural organization, metabolic processes, tissue remodeling and bone development were downregulated in subtypes 4 and 5. Genes related to prostate/breast development and organic/carboxylic acid catabolic processes were downregulated in subtype 5. The differential expression analysis reveals significant gene expression differences among the five cancer subtypes identified by SMMSN. Functional enrichment analysis shows distinct up- and downregulation patterns across the subtypes, highlighting their biological divergence.