Multi-omics prognostic marker discovery and survival modelling: a case study on multi-cancer survival analysis of women's specific tumours.

Mult-omics data was obtained from The Cancer Genome Atlas (TCGA) ( https://tcga.xenahubs.net ) using the UCSCXenaTools R package. Table 1 summarizes the number of samples available for each omics data type across the studied cancers, focusing on four women-related cancers: Breast Invasive Carcinoma (BRCA), Ovarian Serous Cystadenocarcinoma (OV), Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma (CESC), and Uterine Corpus Endometrial Carcinoma (UCEC). For each cancer type, data were collected on gene expression (GE), copy number variations (CNV), DNA methylation (DM), miRNA expression (ME), and clinical variables, including patient vital status and the number of days to death or last follow-up, which were integrated for survival analysis. To ensure consistency, only samples with complete data across all five categories were included. Age was considered as a potential confounding variable, and using the BRCA dataset, we divided the age range (26–89 years) into three quantile-based groups (Supplementary Fig. 1). A Kruskal–Wallis test revealed no significant differences (p > 0.05) in the expression levels of well-studied biomarkers such as miR-22 25 and miR-150 26 across these groups. Kaplan–Meier survival analysis (Supplementary Fig. 2) stratified by age also showed no significant survival differences, suggesting age does not confound the relationship between these biomarkers and survival. These findings, detailed in the Supplementary Material, support the exclusion of age as a primary variable in this study. While similar analyses can be conducted for other factors, our focus remains on omics-based integration, with demographic and exposure-related variables considered as potential additional data modalities. Table 1 Overview of the TCGA sample counts for each omics data type across various cancer types. TCGA cancer types Gene expression (GE) Copy number variation (CNV) DNA methylation (DM) MiRNA expression (ME) Clinical Common BRCA: breast invasive carcinoma 1218 1080 888 832 1247 611 OV: ovarian serous cystadenocarcinoma 308 579 616 485 630 287 CESC: cervical squamous cell carcinoma and endocervical adenocarcinoma 308 295 312 311 313 289 UCEC: uterine corpus endometrial carcinoma 201 539 478 430 596 167 The “Clinical” column represents samples with available survival data, while “Common” indicates patient samples with complete data across all five data types. Overview of the TCGA sample counts for each omics data type across various cancer types. The “Clinical” column represents samples with available survival data, while “Common” indicates patient samples with complete data across all five data types. Each omics modality, including clinical data, was obtained from the UCSC Xena platform using the TCGA Data Hub via the UCSCXenaTools R package. GE data were generated using the Illumina HiSeq 2000 RNA-seq platform and provided as log2(x + 1) transformed RSEM-normalized counts. These values represent gene-level (not isoform-level) expression of protein-coding genes and exclude non-coding RNAs such as microRNAs (miRNAs). Gene annotations were mapped using the UCSC Xena HUGO probeMap. CNV data were obtained at the gene level, processed through the TCGA FIREHOSE pipeline using the GISTIC2 algorithm. Segmented CNV profiles were mapped to individual genes to produce gene-level copy number estimates. These values were further discretized by GISTIC2 into thresholds of − 2, − 1, 0, 1, and 2, corresponding to homozygous deletion, single-copy deletion, diploid normal copy, low-level amplification, and high-level amplification, respectively 27 . DM data consist of beta values (0–1) from Illumina 450 K/27 K assays, where 0 signifies no methylation and 1 indicates full methylation. miRNA expression (ME) data were quantified by RNA-seq using the Illumina HiSeq platform. For each sample, expression values of all isoforms corresponding to the same mature miRNA strand were summed and log2(RPM + 1) transformed to obtain the final expression values. Omics data were structured in tabular format, with samples as rows and features as columns. Only samples labeled “01” (Primary Solid Tumor) were retained for consistency across cancer types. Clinical survival data were integrated into each omics matrix. For GE features with more than 20% missing values were removed, and the top 10% most variable genes were selected using a 90th percentile variance threshold. ME features with over 20% missing values were excluded, and only miRNAs present in more than 50% of samples (with non-zero expression) and in more than 10% of samples with expression values greater than 1 were retained, following Zhao et al. 1 . No variance threshold was applied to miRNA data. CNV data, processed via the GISTIC2 algorithm, were already discretized into gene-level values ranging from − 2 to + 2, and contained no missing values, requiring no further imputation or scaling. DM data were restricted to 27 k CpG probes to enable cancer comparisons, particularly OV, which had only 27 k-level data available. Features and samples with over 50% missing values were removed. While no global variance filtering was applied to DM data, an additional 90th percentile variance threshold was used to reduce dimensionality for OV and UCEC, where a large number of features remained. Missing values in GE, ME, and DM were imputed using the mean of each feature. Min–max scaling was applied to GE and ME data, but not to CNV or DM. Supplementary Tables 1 to 4 summarize the number of features and samples retained after preprocessing for each omics type. Pre-processed TCGA data were then passed through our feature selection pipeline to identify the final feature set for each omics, which was evaluated using predictive survival models. Before developing the pipeline, we assessed filter methods against a no-selection baseline (Fig.  1 ). Supplementary Table 5 details the machine learning algorithms, filter methods, R packages, and hyperparameters used. Fig. 1 Evaluation of Individual Filter Methods Pipeline. Workflow for evaluating each filter method independently. Cross-validation is employed to assess the performance of each method, and the results are summarized in a heatmap displaying the mean performance and mean feature selected by each filter method when paired with various predictive models. Evaluation of Individual Filter Methods Pipeline. Workflow for evaluating each filter method independently. Cross-validation is employed to assess the performance of each method, and the results are summarized in a heatmap displaying the mean performance and mean feature selected by each filter method when paired with various predictive models. Figure  1 illustrates the evaluation process for each feature selection method. We applied 5 repeats of fivefold cross-validation, splitting the data into five folds, with four for training and one for validation. Each filter method was applied to the training set, selecting features based on thresholded coefficients and importance scores. For the univariate and multivariate Cox methods, features were selected if they were significant (p < 0.05) and had non-zero coefficients. In the univariate Cox method specifically, p-values were further adjusted for multiple testing using the false discovery rate (FDR) correction, and features passing an FDR-adjusted threshold of 0.05 were retained. For the random survival forest method, features with an importance score greater than 0 were selected. The performance of these selected features was assessed using four survival machine learning algorithms, evaluating the concordance index (C-index) on the validation set. The C-index measures model discrimination, ranging from 0.5 (random guessing) to 1 (perfect prediction) 1 . This process was repeated 5 times, resulting in 25 evaluations per feature selection method and survival model. The final performance of each method was determined by averaging the C-index values. We developed two feature selection pipelines, one that utilizes cross-validation and another that uses bootstrapping. Figure  2 A illustrates our cross-validation (CV) pipeline. We begin with a 70:30 training-test split, using the training set for feature selection. The training set is divided into 5 folds, with one-fold held out each time. Each filter method is applied to the remaining folds to identify important features. This process is repeated 5 times with fivefold CV, resulting in 25 iterations. Features are ranked by their occurrence in the selected sets, with a maximum possible occurrence of 100 (25 iterations × 4 filter methods). Features selected at least 50% of the time are included in the final set. This final set is then evaluated on the test set to ensure no information leakage. The second feature selection method, the bootstrapping pipeline (Fig.  2 B), involves bootstrapping 70% of the data and applying each filter method. The common features selected across all methods form a bootstrapped feature set, repeated 100 times. We track the frequency of each feature’s occurrence, selecting those that appear at or above the 70th quantile. This approach identifies highly predictive features while accounting for data variability. Fig. 2 Feature Selection Pipelines. ( A ) CV Feature selection pipeline. Apply a cross-validation process combined with a voting technique to identify optimal features. Selected features are determined based on occurrence across the chosen filter methods, enhancing the robustness of feature selection. ( B ) Bootstrapping Feature selection pipeline. Employs bootstrapping to improve the reliability of selected features by resampling the data across numerous iterations and applying the same voting technique. Feature Selection Pipelines. ( A ) CV Feature selection pipeline. Apply a cross-validation process combined with a voting technique to identify optimal features. Selected features are determined based on occurrence across the chosen filter methods, enhancing the robustness of feature selection. ( B ) Bootstrapping Feature selection pipeline. Employs bootstrapping to improve the reliability of selected features by resampling the data across numerous iterations and applying the same voting technique. To evaluate both pipelines, we test the final selected features using the separate test set. The training set, filtered to include only the selected features, is used to train predictive models, and performance is assessed on the independent test set. This ensures unbiased evaluation and robust assessment of predictive performance. All model training and evaluation were performed using 30 bootstrap iterations. Following feature selection for each omics modality, we integrate multi-omics data using feature-level fusion, refining the signature panel to enhance predictive performance while ensuring clinical feasibility (Fig.  3 ). To achieve this, we employ two refinement strategies. In One-Stage Refinement, selected features from each omics modality are first concatenated, and Recursive Feature Elimination (RFE) is applied to identify the most predictive features. In contrast, Two-Stage Refinement applies RFE at the single-omics level before fusion, then again at the multi-omics level, ensuring that only the most informative features are retained at each stage. This sequential approach enhances model robustness, minimizes overfitting, and improves overall predictive accuracy. To determine the optimal combination of omics modalities for survival prediction, we systematically evaluate all possible integrations, similar to Tong et al. 19 who used autoencoders for feature fusion. We assess model performance using various machine learning methods, including Cox Proportional Hazards, ElasticNet, Random Forest (RF), and GLMBoost. Ultimately, RF is selected for multi-omics integration due to its superior predictive accuracy (Fig.  4 , Table 2 ). Fig. 3 Multi-omics pipelines. ( A ) One-Stage Refinement integration pipeline, where two modalities are combined at the feature level. Recursive Feature Elimination (RFE) is then applied to reduce dimensionality, further refining the feature set to enhance model performance. ( B ) Two-Stage Refinement integration pipeline, where individual modalities first undergo RFE to reduce dimensionality. The reduced feature sets are then fused together, followed by an additional round of RFE on the fused data to obtain the final set of multi-omics features. Fig. 4 Heatmap summarizing feature selection outcomes and model performance using CESC DM data. ( A ) Mean number of features selected by each feature selection method across models based off importance. ( B ) Mean concordance index (C-index) achieved by each model–feature selection combination. The colour gradient (light to dark purple) reflects increasing values—either the number of features ( A ) or C-index ( B ). Error bars are not shown in the figure but corresponding 95% confidence intervals for all mean values are reported in the Supplementary Tables. Abbreviations: RF = Random Forest; No FS = No Feature Selection; Uni/Multivariate = Univariate/Multivariate CoxPH; Var Imp = Variable Importance; Min Depth = Minimal Depth; Max Stat = Variable Hunting; GLMBoost = Boosted CoxPH; ElasticNet = Penalized CoxPH. Table 2 Performance evaluation of feature selection pipelines applied to the BRCA dataset. BRCA Omics None CV BS CoxPH (1) GE 0.461 ± 0.0406 (0.393–0.531) 0.447 ± 0.0908 (0.313–0.634) 0.351 ± 0.066 (0.25–0.465) ME 0.558 ± 0.115 (0.366–0.738) 0.454 ± 0.0936 (0.277–0.615) 0.545 ± 0.0654 (0.413–0.67) DM 0.453 ± 0.0828 (0.327–0.61) 0.434 ± 0.0812 (0.302–0.591) 0.354 ± 0.0528 (0.266–0.477) CNV 0.467 ± 0.0608 (0.393–0.599) 0.554 ± 0.0574 (0.459–0.634) 0.506 ± 0.0446 (0.441–0.577) Ranger (2) GE 0.578 ± 0.0799 (0.456 – 0.705) 0.546 ± 0.0682 (0.443 – 0.674) 0.604 ± 0.0794 (0.436 – 0.704) ME 0.605 ± 0.0745 (0.458 – 0.727) 0.557 ± 0.0766 (0.426 – 0.651) 0.52 ± 0.08 (0.369 – 0.645) DM 0.638 ± 0.0771 (0.49 – 0.774) 0.605 ± 0.0819 (0.44 – 0.722) 0.611 ± 0.0514 (0.528 – 0.698) CNV 0.517 ± 0.0639 (0.418 – 0.632) 0.52 ± 0.056 (0.417 – 0.612) 0.549 ± 0.0616 (0.429 – 0.659) GLMBoost (3) GE 0.377 ± 0.0851 (0.212–0.521) 0.327 ± 0.0759 (0.212–0.465) 0.278 ± 0.0635 (0.183–0.387) ME 0.472 ± 0.0856 (0.337–0.625) 0.414 ± 0.0772 (0.281–0.549) 0.477 ± 0.0757 (0.369–0.642) DM 0.424 ± 0.072 (0.298–0.543) 0.428 ± 0.0923 (0.304–0.62) 0.564 ± 0.0791 (0.459–0.734) CNV 0.513 ± 0.0562 (0.416–0.611) 0.489 ± 0.0637 (0.357–0.592) 0.506 ± 0.0757 (0.391–0.639) ElasticNet (4) GE 0.355 ± 0.0643 (0.252–0.46) 0.304 ± 0.0533 (0.215–0.401) 0.5 ± 0 (0.5–0.5) ME 0.406 ± 0.0915 (0.406 ± 0.0915) 0.406 ± 0.0915 (0.338–0.59) 0.43 ± 0.101 (0.271–0.605) DM 0.455 ± 0.0725 (0.298–0.565) 0.455 ± 0.0697 (0.341–0.581) 0.615 ± 0.0718 (0.486–0.725) CNV 0.5 ± 0 (0.5–0.5) 0.5 ± 0 (0.5–0.5) 0.452 ± 0.0892 (0.295–0.614) Results are shown for baseline models without feature selection (“None”) and models using cross-validation (CV) or bootstrapping (BS) based feature selection. Model abbreviations: Ranger = Random Forest, GLMBoost = Boosted Cox model, ElasticNet = Penalized Cox model. Performance is reported as the mean concordance index (C-index) ± standard deviation, with 95% confidence intervals provided in parentheses. Overall highest performing model values are in bold. Multi-omics pipelines. ( A ) One-Stage Refinement integration pipeline, where two modalities are combined at the feature level. Recursive Feature Elimination (RFE) is then applied to reduce dimensionality, further refining the feature set to enhance model performance. ( B ) Two-Stage Refinement integration pipeline, where individual modalities first undergo RFE to reduce dimensionality. The reduced feature sets are then fused together, followed by an additional round of RFE on the fused data to obtain the final set of multi-omics features. Heatmap summarizing feature selection outcomes and model performance using CESC DM data. ( A ) Mean number of features selected by each feature selection method across models based off importance. ( B ) Mean concordance index (C-index) achieved by each model–feature selection combination. The colour gradient (light to dark purple) reflects increasing values—either the number of features ( A ) or C-index ( B ). Error bars are not shown in the figure but corresponding 95% confidence intervals for all mean values are reported in the Supplementary Tables. Abbreviations: RF = Random Forest; No FS = No Feature Selection; Uni/Multivariate = Univariate/Multivariate CoxPH; Var Imp = Variable Importance; Min Depth = Minimal Depth; Max Stat = Variable Hunting; GLMBoost = Boosted CoxPH; ElasticNet = Penalized CoxPH. Performance evaluation of feature selection pipelines applied to the BRCA dataset. Results are shown for baseline models without feature selection (“None”) and models using cross-validation (CV) or bootstrapping (BS) based feature selection. Model abbreviations: Ranger = Random Forest, GLMBoost = Boosted Cox model, ElasticNet = Penalized Cox model. Performance is reported as the mean concordance index (C-index) ± standard deviation, with 95% confidence intervals provided in parentheses. Overall highest performing model values are in bold. RFE is applied iteratively, ranking features based on RF variable importance scores and progressively eliminating the least informative ones. The final multi-omics signature is determined by maximizing predictive performance, measured by the c-index. By refining feature selection both within individual omics and across the integrated dataset, our approach optimizes the signature panel, reducing its size while preserving predictive power. This ensures a more interpretable and clinically applicable model for survival prediction in cancer. With the optimal multi-omics combinations identified for BRCA, OV, CESC, and UCEC, we now turn to our comparative cancer analysis to uncover shared multi-modal signatures across these cancers. First, we will conduct an overlap analysis to identify common features across the four cancers and explore potential overlaps by analyzing miRNA targets associated with these features. Next, we will examine disease association networks to understand the relationships between these overlapping targets, offering insights into their interconnected roles in cancer biology. To uncover shared biological themes, we will perform Gene Set Enrichment Analysis (GSEA) on the targets from each cancer, focusing on common Gene Ontology (GO) terms and KEGG pathways 28 – 30 . Finally, to ensure the robustness and biological relevance of our findings, we will explore key features using Kaplan–Meier plots and cross-reference them with existing literature. While our feature selection was driven by predictive power in survival analysis, this validation step will confirm that these features have meaningful biological implications. Through this integrated approach, we aim to demonstrate that our results are not only statistically significant but also provide valuable biological insights into the commonalities and differences across the cancers under study.

All analyses were conducted on the University of New South Wales (UNSW) Katana high-performance computing platform, utilizing 16 CPUs and 124 GB of memory per job to efficiently process large-scale omics data. While optimized for HPC, the pipeline can also run on local machines by reducing data size or enabling parallelization for improved runtime efficiency. Before developing our feature selection pipeline, we evaluated the four filter methods used in this study (Fig.  1 ). Due to computational constraints, we limited our analysis to these four methods. The results are presented as heatmaps showing the mean C-index values and the number of features selected across 5 repeats of fivefold cross-validation for each combination of ML algorithms (rows) and filter methods (columns). The Cox Proportional Hazards (CoxPH) model is included as the baseline in our benchmarking, due to its longstanding role as a standard survival modelling approach. Column 1, representing results without feature selection, serves as the baseline for comparison. Figure  4 presents the results from the CESC DM dataset, which reflect trends consistent across other modalities and cancer types (Supplementary Tables 6 to 9), serving as a representative example. We expected the performance of each filter method to at least match the baseline (no feature selection). As shown in Fig.  4 A, feature selection substantially reduced the number of features, while predictive performance (Fig.  4 B) remained stable or slightly declined, indicating effective removal of redundancy. Random Forest generally performed best across methods but dropped to a C-index of 0.39 when paired with the univariate Cox filter, which selected only 2–3 features on average after FDR correction—insufficient for models that depend on multivariate interactions. ElasticNet and GLMBoost also underperformed relative to the baseline CoxPH model without feature selection. This may be due to over-shrinkage from regularization penalties, which can exclude informative features and lead to underfitting. Supporting this, both models retained fewer features than RF and CoxPH (Fig.  4 A). Due to its inconsistent effects and limited utility across models, the univariate filter was excluded from further analyses. In this study, we developed two feature selection pipelines: one based on cross-validation and the other on bootstrapping , both utilizing our four main filter methods. Table 2 compares the performance of these pipelines to a baseline with no feature selection. Overall, both the cross-validation and bootstrapping pipelines demonstrated strengths depending on the context. While cross-validation generally showed robust and consistent improvements, in certain cases, bootstrapping outperformed cross-validation, particularly for specific omics types or cancer datasets. Random Forest consistently delivered the best results across most settings, confirming its effectiveness as a feature selection method. Moving forward, for each omics data type and cancer cohort, we will adopt the feature selection pipeline—either bootstrapping or cross-validation—that yields the superior performance, ensuring an optimized and tailored approach for each scenario. Supplementary Tables 10 to 12 show the results for the other cancers investigated. Based on these observations (Supplementary Tables 13 to 16), we selected the optimal refinement strategy for each cancer type. In three out of four cancers, the Two-Stage refinement strategy yielded superior performance. In contrast, for BRCA, the One-Stage strategy performed only slightly better, making it the preferred choice for that specific case. Table 3 summarizes the performance of each multi-modality combination across all cancer types, demonstrating consistent improvements over their corresponding single-modality Random Forest baselines. The enhanced performance is further supported by the key features driving each model: Supplementary Tables 17–20 list the contributing feature sets from each omics layer. To support the biological validity of our signatures, Supplementary Tables 21–24 highlight well-established biomarkers that overlap with or are associated with our model outputs. Table 3 Performance of each modality combination for all cancers using the best refinement strategy. Cancer ME + GE ME + DM GE + DM GE + CNV ME + CNV BRCA 0.698 0.643 0.682 0.683 0.661 OV 0.608 0.618 0.580 0.585 0.537 CESC 0.737 0.754 0.639 0.600 0.743 UCEC 0.716 0.651 0.707 0.700 0.678 Cancer DM + CNV DM + GE + ME DM + GE + CNV DM + ME + CNV ME + GE + CNV ME + GE + CNV + DM BRCA 0.642 0.672 0.672 0.674 0.688 0.692 OV 0.599 0.581 0.592 0.605 0.608 0.609 CESC 0.643 0.734 0.631 0.750 0.745 0.730 UCEC 0.677 0.730 0.707 0.734 0.726 0.754 Highest C-index value for each cancer and its modality combination are in bold. Performance of each modality combination for all cancers using the best refinement strategy. Highest C-index value for each cancer and its modality combination are in bold. For OV, which exhibited the lowest overall performance among the cancers, the optimal modality combination was ME and DM, achieving a c-index of 0.618. CESC shared the same optimal combination, yielding a higher c-index of 0.754. For BRCA, the best performance was achieved using ME and GE (c-index = 0.698), while UCEC attained its highest c-index (0.754) when integrating all four modalities despite involving fewer features than other cancers. Importantly, ME appeared in the optimal combination for every cancer type, suggesting its strong and consistent contribution to multi-omics integration in survival prediction. After obtaining the final set of multi-modal features for each cancer using our multi-omics pipeline, we conducted a comparative cancer analysis to identify overlapping features among cancers. Table 4 presents the results of our direct overlap analysis, highlighting shared molecular features across the four cancer types. The highest number of overlapping features was observed between BRCA and OV, with all shared features corresponding to miRNA signatures. Interestingly, no single feature was found to be common across all four cancers. However, a small number of miRNAs were shared among three cancer types. Specifically, miR-150 and miR-9 were common to BRCA, CESC, and UCEC, while miR-4662a and miR-1287 were shared by BRCA, OV, and UCEC. These results suggest that miRNA-based features dominate the shared signature space across cancers and may reflect conserved post-transcriptional regulatory mechanisms in these female-specific malignancies. Table 4 Directly overlapping features across cancers. Cancer types Overlapping features BRCA, OV ME_hsa.miR.186.5p, ME_hsa.miR.191.5p, ME_hsa.miR.4662a.5p, ME_hsa.miR.542.3p, ME_hsa.miR.184, ME_hsa.miR.616.3p, ME_hsa.miR.1287.5p, ME_hsa.miR.486.5p BRCA, CESC ME_hsa.miR.150.5p, ME_hsa.miR.628.5p, ME_hsa.miR.378a.3p, ME_hsa.miR.451a, ME_hsa.miR.9.5p BRCA, UCEC ME_hsa.miR.146a.5p, ME_hsa.miR.150.5p, ME_hsa.miR.4662a.5p, ME_hsa.miR.671.3p, ME_hsa.miR.1287.5p, ME_hsa.miR.9.5p OV, CESC ME_hsa.miR.378a.5p OV, UCEC ME_hsa.miR.1224.5p, ME_hsa.miR.363.3p, ME_hsa.miR.1287.5p, ME_hsa.miR.4662a.5p CESC, UCEC ME_hsa.miR.150.5p, ME_hsa.miR.9.5p BRCA, OV, UCEC ME_hsa.miR.4662a.5p, ME_hsa.miR.1287.5p BRCA, CESC, UCEC ME_hsa.miR.150.5p, ME_hsa.miR.9.5p ME = miRNA. Directly overlapping features across cancers. ME = miRNA. The goal of this analysis is to uncover further overlaps among our cancers beyond the direct feature overlap comparison. To investigate functional overlap across cancers, we examined the gene targets of selected miRNA features using the multiMiR R package. We restricted results to experimentally validated interactions by querying the miRecords, miRTarBase, and TarBase databases via the summary = TRUE setting. Only gene targets with supporting experimental evidence were retained. To ensure biological relevance, we further filtered these targets to include only genes expressed in the corresponding cancer transcriptomic datasets. These filtered gene sets were used to assess inter-cancer overlap at the target level (Fig.  5 ). Fig. 5 Upset plot showing direct overlapping miRNA Targets between cancers. Upset plot showing direct overlapping miRNA Targets between cancers. We identified that 143 gene targets were shared in all four cancers. Notably, BRCA contributed the most gene targets, while CESC provided the least. Using the 143 overlapping gene targets, we performed enrichment analysis using DisGeNET, a database of gene-disease associations. Figure  6 A presents an enrichment map that highlights a large, interconnected disease cluster, where genes with shared associations group into functionally related modules. Within this overarching cluster, we observe that several cancers relevant to the gynecologic system and hormone-regulated tissues—such as endometrial neoplasms, endometrial adenocarcinoma, invasive carcinoma of breast, and noninfiltrating intraductal carcinoma—are densely connected. These are highly relevant to BRCA, UCEC, and CESC, which originate from hormonally responsive tissues and are frequently driven by dysregulation in estrogen, progesterone, and other hormone signaling pathways 31 . Other disease terms like endometriosis of the ovary, dyslipidemias, and acute coronary syndrome may reflect systemic inflammatory or metabolic comorbidities often reported in patients with gynecologic cancers, especially those undergoing hormonal therapy or post-menopausal transitions 32 . The cluster involving urogenital disorders (e.g. pyelonephritis, kidney calculi) may point toward anatomical proximity or shared pathways in epithelial tissue remodeling and inflammation—both critical in ovarian and cervical cancer progression 33 . Fig. 6 Disease network using the overlapping targets and miRNA. ( A ) Enrichment map of overlapping miRNA targets disease network. ( B ) Overlapping miRNA disease associations across cancers. Disease network using the overlapping targets and miRNA. ( A ) Enrichment map of overlapping miRNA targets disease network. ( B ) Overlapping miRNA disease associations across cancers. Figure  6 B focuses on the overlapping miRNAs identified in our direct overlap analysis (Table 4 )—those present in at least two of the studied cancers. To assess their disease relevance, we mapped these miRNAs to disease associations using the miR2Disease database. The network shows clear connections between key miRNAs and multiple malignancies. For example, hsa-miR-9-5p, detected in BRCA, UCEC, and CESC, is linked to neuroblastoma, melanoma, and glioma, but has also been implicated in EMT regulation and metastasis in breast and cervical cancers 34 , 35 . Notably, hsa-miR-378a-3p and hsa-miR-378a-5p are connected to gastric, kidney, and colorectal cancers, all of which share similar stromal remodeling and angiogenic programs with gynecological tumors 36 . Additionally, miRNAs such as hsa-miR-146a-5p and hsa-miR-105-5p show associations with prostate and leukemia subtypes, indicating broader systemic roles in inflammation and immune regulation—processes that are also critical in the tumor microenvironment of CESC and UCEC 37 , 38 . These associations lend biological credibility to the miRNA signatures found in multiple cancers, reinforcing their potential as shared or cooperative biomarkers across tissue types. We performed Gene Set Enrichment Analysis (GSEA) 28 – 30 on our miRNA targets to identify shared biological themes and validate our signatures across cancers. Figure  7 A shows enrichment analysis for BRCA miRNA targets, with GO terms such as mesenchymal to epithelial transition, epidermis development, and calcium-regulated neurotransmitter exocytosis—processes linked to cellular plasticity, differentiation, and intercellular communication, all relevant to breast cancer progression. KEGG pathways like ECM-receptor interaction and IL-17 signaling are well-established in tumor invasion, immune modulation, and microenvironment remodelling 39 , 40 . Figure  7 B presents OV enrichment results, highlighting developmental and neurobiological GO terms such as embryonic morphogenesis, axon guidance, and synaptic signaling, reflecting stem-like traits and emerging roles of tumor-nerve interactions in ovarian cancer 41 , 42 . KEGG terms include cell adhesion molecules and cAMP signaling, supporting roles in peritoneal spread and chemoresistance 43 . Figure  7 C displays CESC results, with GO terms including reproductive system development, axon guidance, and neuroinflammatory signaling—linked to HPV-driven transformation and immune infiltration in cervical tumors 44 . KEGG pathways such as IL-17 signaling, cytokine interactions, and complement cascades reflect inflammation and stromal remodelling central to cervical cancer pathogenesis 45 , 46 . Figure  7 D shows UCEC targets enriched for GO terms related to axon extension, thymic T cell differentiation, and embryonic development, suggesting roles in immune modulation and tumor plasticity. KEGG pathways like Wnt signaling—a key driver of endometrial tumorigenesis—and the renin-angiotensin system emphasize hormone-sensitive and angiogenic features of UCEC 47 , 48 . Fig. 7 GSEA enrichment plots of Top 10 GO Terms and KEGG Pathways. ( A ) BRCA, ( B ) OV, ( C ) CSEC and ( D ) UCEC. GSEA enrichment plots of Top 10 GO Terms and KEGG Pathways. ( A ) BRCA, ( B ) OV, ( C ) CSEC and ( D ) UCEC. After reviewing the GSEA results for each cancer’s miRNA targets and identifying notable trends, we analyzed the overlapping GO Terms and KEGG Pathways across the cancers to identify shared biological themes and pathways, further elucidating underlying mechanisms and potential therapeutic targets. Figure  8 A shows an upset plot of GO Term and KEGG Pathway overlaps, revealing the intersection of all four cancers (BRCA, OV, UCEC, and CESC) sharing 73 common GO terms and only 2 KEGG Pathways. Figure  8 B highlights conserved biological processes and signaling pathways derived from miRNA target genes shared across BRCA, OV, CESC, and UCEC. GO terms such as ureteric bud morphogenesis, mesonephric tubule development, and branching morphogenesis suggest activation of developmental programs, consistent with the hallmark of cancer-related reprogramming that recapitulates embryonic organogenesis 49 . Enrichment of terms like neuron projection guidance, axonogenesis, and neurogenesis points to shared neurodevelopmental signatures, in line with evidence linking axon guidance molecules (e.g., SLITs, semaphorins) to cancer migration, perineural invasion, and metastasis 50 , 51 . Terms such as mesenchymal to epithelial transition and nervous system development further support tumor plasticity and reactivation of embryonic networks. On the KEGG side, IL-17 signaling appears consistently across all cancers, reflecting its role in inflammation-driven tumorigenesis and immune evasion 40 , while enrichment of cell adhesion molecules (CAMs) highlights shared mechanisms in EMT, tumor-stroma interaction, and metastatic potential 52 . Together, these conserved pathways underscore the functional convergence of developmental, immune, and cell communication programs in gynecologic and breast cancers, reinforcing the biological relevance of the core miRNA targets. Fig. 8 ( A ) Upset plot showing GSEA overlaps across all cancers for miRNA targets. ( B ) Enrichment dot plot showing Top 10 GSEA overlaps across all cancers. ( A ) Upset plot showing GSEA overlaps across all cancers for miRNA targets. ( B ) Enrichment dot plot showing Top 10 GSEA overlaps across all cancers. To assess the clinical relevance of our shared multi-modal features, we conducted Kaplan–Meier survival analyses and selected miR-150 and miR-9 as representative examples due to their recurrence in BRCA, CESC, and UCEC, and their stronger prognostic signals. miR-150 is a well-characterized tumor suppressor that targets oncogenes such as MYB and AKT2, and its higher expression has been linked to improved survival in triple-negative breast cancer 53 . In CESC, miR-150 regulates cell growth via FOXO4 and has shown similar tumor-suppressive behaviour 53 . Consistently, our analysis (Fig.  9 A) showed that low miR-150 expression correlated with poorer survival in both BRCA and CESC. miR-9, by contrast, is a context-dependent miRNA frequently associated with tumor progression. In UCEC, we observed that high miR-9 expression was linked to worse survival (Fig.  9 B), suggesting a potential oncogenic role. This aligns with previous findings showing that miR-9 promotes metastasis by targeting E-cadherin and FOXO1, and its high expression predicts poor outcomes across several cancers 54 , 55 . These findings underscore the biological and clinical importance of our core miRNA features, particularly miR-150 and miR-9, in driving survival differences across cancer types. Fig. 9 Kaplan–Meier survival plots ( A ) miR-150, ( B ) miR-9. Kaplan–Meier survival plots ( A ) miR-150, ( B ) miR-9.

PRISM was developed with two primary objectives: (1) to improve cancer survival prediction using multi-omics data and robust feature selection strategies, and (2) to identify a concise set of multi-modal biomarkers most predictive of patient outcomes. It offers a comprehensive framework that systematically integrates feature selection and survival modeling to advance prognostic marker discovery. By benchmarking multiple feature selection methods in combination with various survival models, PRISM enables rigorous and reproducible performance assessment. The pipeline incorporates both cross-validation and bootstrapping, along with ensemble-based voting, to ensure robustness across patient cohorts. Crucially, PRISM applies recursive feature elimination to distill high-dimensional omics data into a compact and non-redundant feature set—prioritizing minimal biomarker panels that maintain or exceed the predictive power of full high-throughput datasets. This is central to its design philosophy: enabling clinically feasible and cost-effective biomarker deployment. Through a two-stage refinement approach, PRISM first identifies informative features within each omics layer, then refines across modalities to construct optimal single- and multi-omics signatures. As an end-to-end solution, it automates the workflow from data acquisition and preprocessing to survival modeling and biological validation, supporting broad applicability beyond TCGA and toward translational implementation in real-world clinical settings. We drew inspiration from previous studies, including Zhao et al. 1 , who identified pan-cancer prognostic biomarkers through multi-omics integration, and Tong et al. 19 , who explored deep learning-based feature integration for breast cancer survival. Spooner et al. 3 further informed our approach to machine learning–based survival modeling and feature selection. We compared our results to Tong et al. 19 , who proposed a multi-view learning framework for breast cancer survival prediction by integrating GE, DM, ME, and CNVs using complementary (ConcatAE) and consensus (CrossAE) autoencoders. These models preserved modality-specific features while reducing cross-omics noise. Their best performance for breast cancer overall survival was achieved using DM and ME, yielding a C-index of 0.641 ± 0.031 with ConcatAE and 0.630 ± 0.081 with CrossAE. Wen and Li 56 applied MMOSurv, a meta-learning framework based on deep Cox models, to enable few-shot survival prediction using ME and GE data. By leveraging meta-knowledge across related cancer types, their model achieved C-index values of 0.662 for BRCA, 0.695 for CESC, and 0.713 for UCEC. Cheerla and Gevaert 17 developed a multimodal neural network integrating clinical data, ME and GE, and histopathology images using modality-specific encoders and unsupervised feature fusion. Their model reported C-index scores of 0.73 for BRCA, 0.74 for CESC, 0.59 for OV, and 0.66 for UCEC using their baseline configuration for ME and GE integration. Even when restricting our analysis to the same modality combinations used in previous studies, our results are not only competitive but often superior (Table 3 ). Using the ME and GE integration, PRISM achieved C-index scores of 0.698 for BRCA, 0.608 for OV, 0.737 for CESC, and 0.716 for UCEC—surpassing the reported benchmarks from models such as ConcatAE 19 and MMOSurv 56 . Importantly, our framework is not limited to these fixed modality pairings; by systematically evaluating cross-modality combinations, PRISM identified combinations beyond ME and GE that yielded even higher performance, achieving a high C-index of 0.754 for both CESC and UCEC. Unlike deep learning models that require the full high-throughput omics profile, PRISM derives compact, robust, and interpretable multi-modal signatures. This not only enhances performance but also aligns with the goal of clinical feasibility, offering a cost-effective and scalable solution for real-world prognostic applications. Contrary to conventional expectations, GE was not among the most informative modalities in two of the four cancer types analysed. Instead, ME consistently emerged as the top-performing modality in survival prediction. This observation aligns with findings from Tong et al. 19 and Cheerla and Gevaert 17 , both of whom reported that miRNA data outperformed mRNA-based features in survival analysis when comparing modality contributions. From a biological perspective, this is well supported: miRNAs are small non-coding RNAs that regulate gene expression post-transcriptionally, often targeting multiple mRNAs and exerting broad control over key cellular processes involved in tumor progression. Their dysregulation has been repeatedly linked to oncogenesis, metastasis, and therapy resistance, making them strong candidates for prognostic biomarker development 57 , 58 . As such, the prominence of ME in our models likely reflects its ability to capture regulatory and functional changes that are not always evident at the mRNA expression level. We assessed the biological relevance of our model outputs through disease associations using multiMiR (Fig.  6 ), targeted literature validation (Supplementary Tables 21–24), and gene set enrichment analysis (Fig.  7 ), confirming that our signatures are firmly grounded in biological context. Despite these encouraging results, several areas warrant improvement before expanding this study. First, while we adopted pre-processing strategies similar to Zhao et al. 1 and Tong et al. 19 , our use of low-variance filtering may have excluded biologically significant markers—such as BRCA1 and BRCA2 in the BRCA dataset—due to their low variance in the population. More nuanced filtering approaches will be needed to mitigate this. Another limitation is that feature selection was performed independently for each omics modality prior to integration, potentially overlooking features that lack individual significance but contribute jointly. However, as PRISM is designed to identify compact, clinically feasible biomarker panels, incorporating joint modeling of the entire feature space would compromise the interpretability and cost-effectiveness of our approach. While our framework delivered competitive C-index values comparable to deep learning models using full omics profiles, PRISM offers a more clinically translatable solution, achieving similar performance using a reduced set of low-cost biomarkers. Additionally, our use of simple imputation (e.g., mean imputation) may have limited data quality; more sophisticated methods such as Multiple Imputation by Chained Equations (MICE) could be considered. Finally, we did not incorporate differentially expressed genes (DEGs) or differentially methylated regions (DMRs), which could enhance the biological richness of gene and methylation feature sets. Addressing these limitations will be essential to further improve the robustness, interpretability, and translational potential of our approach. For our models, particularly GLMBoost and Elastic Net, observed performance was somewhat lower than anticipated, which may be partly attributable to over-shrinkage of feature effects. While we used the default boosting parameters for GLMBoost, which controls model complexity primarily through the number of boosting iterations, the Elastic Net model employed regularization with an alpha value of 0.5, balancing L1 and L2 penalties. We acknowledge that more thorough hyperparameter optimization—such as tuning the number of boosting iterations (mstop) for GLMBoost and systematically exploring the alpha and lambda parameters for Elastic Net—could potentially improve model performance. An important consideration for future work is the potential integration of deep learning methods to enhance multi-omics data analysis. However, PRISM aims to identify a minimal yet highly predictive multi-omics signature panel to improve clinical utility and commercial viability. To achieve this, we implemented a rigorous feature selection strategy that combines machine learning-based selection with recursive feature elimination, optimizing panel size while maintaining competitive predictive performance. In contrast, deep learning models, such as autoencoders and transformer-based models, learn features directly from high-dimensional data without explicit feature selection, often resulting in less interpretable features. While prior studies 2 , 10 , 17 – 19 , 56 successfully applied deep neural networks to multi-omics data, they require the inclusion of all omics features, which can be economically unfeasible due to the high costs of high-throughput omics technologies. Notably, PRISM’s performance on TCGA data is comparable to, and in some cases exceeds, deep learning models. While direct comparisons are limited, these results emphasize the competitive advantage of PRISM. Nevertheless, we acknowledge the future potential of deep learning methods for multi-omics analysis as computational and economic barriers evolve.

Cancer’s complex pathophysiology is shaped by diverse genetic, environmental, and molecular factors, leading to considerable variability in patient outcomes even within the same cancer types, which complicates treatment strategies 1 , 2 . High-throughput molecular profiling technologies, such as next-generation sequencing and mass spectrometry, have become fundamental in precision medicine, enabling comprehensive analysis of DNA, RNA, and proteins to discover biomarkers. However, relying on single omics data provides only partial insights into the intricate mechanisms of cancer, potentially missing critical biomarkers and therapeutic opportunities 3 , 4 . The heterogeneity of cancer, reflected in its diverse subtypes and molecular profiles, requires an integrated approach. Combining multiple omics data types is crucial for gaining a holistic understanding of cancer biology and enabling personalized treatment strategies 5 . The vast array of omics data—spanning genomics, transcriptomics, proteomics, metabolomics, and epigenomics—offers a comprehensive view of cancer biology, with immense potential to identify novel biomarkers and improve clinical outcomes 6 , 7 . However, multi-omics data poses challenges like high dimensionality, data imbalance, noise, missing values, and heterogeneity, complicating robust analysis and biomarker discovery 5 . To address these challenges, integrative machine learning methods such as multimodal learning, ensemble strategies, and network-based approaches are increasingly crucial in advancing biomarker discovery and clinical outcome models for disease diagnosis, prognosis, and treatment response monitoring 8 – 12 . In this context, survival analysis is pivotal, helping researchers and clinicians evaluate factors influencing patient outcomes over time. Predicting survival beyond critical milestones, such as five-year survival, informs treatment decisions and enhances our understanding of cancer progression 3 . This growing recognition of the importance of multi-omics survival modeling has spurred numerous studies in the past decade. Yuan et al. 13 conducted a pioneering study by integrating multiple omics data—including somatic copy-number alterations (SCNA), DNA methylation, mRNA, miRNA, and protein expression—to predict patient survival across four cancer types from The Cancer Genome Atlas (TCGA) project. This work established the value of multi-omics integration for survival analysis. Building on this, Zhu et al. 14 systematically evaluated the prognostic value of various omics types—clinical data, mRNA, SCNA, DNA methylation, miRNA—across 14 cancers using a kernel machine learning approach. They found gene and miRNA expression to have the strongest predictive power. Liu et al. 15 focused on multi-omics variable selection combined with Cox-regression model for cancer prognosis prediction. Subsequently, Chai et al. 2 integrated multi-omics data to derive representative features, which were then utilized in a Cox proportional hazards (Cox-PH) model for enhanced survival prediction. Later, Chaudhary et al. 16 employed an autoencoder to integrate DNA methylation, miRNA expression, and RNA-Seq data, selecting survival-associated features using univariate Cox-PH analysis and identifying robust subgroups in hepatocellular carcinoma (HCC). Similarly, Cheerla et al. 17 utilized a deep learning approach with multi-modal representations to enhance cancer prognosis and Cox-PH for outcome prediction. Hao et al. 18 extended deep learning methodologies by proposing a gene- and pathway-based deep neural network (DNN) where the Cox-PH model served as the output layer, embedding multi-omics features for survival prediction. Tong et al. 19 further advanced the field by introducing multi-view learning with autoencoders, integrating multi-omics data through complementary and consensus principles. Their approach demonstrated how DNA methylation and miRNA expression offer both unique and shared prognostic insights, with the integrated model significantly outperforming single-omics methods. While these studies emphasize the use of multi-omics data for cancer survival modeling, they often overlook clinical feasibility, relying on high-throughput multi-omics profiles as input. Such approaches are currently limited in real-world clinical settings due to cost and logistical constraints, highlighting the need for frameworks that identify compact, cost-effective biomarker panels without sacrificing predictive performance. To address these gaps, we developed PRISM (PRognostic marker Identification and Survival Modelling through Multi-omics Integration), an end-to-end framework for survival prediction and prognostic biomarker identification. PRISM systematically analyzes TCGA multi-omics data, including gene expression (GE), DNA methylation (DM), miRNA expression (ME), and copy number variations (CNV). It employs a comprehensive feature selection and survival modeling pipeline, using statistical and machine learning techniques to extract key biomarkers, integrate them via two fusion methods, and evaluate their predictive power. PRISM benchmarks feature selection methods, such as univariate/multivariate Cox filtering and Random Forest importance, alongside survival models including CoxPH, ElasticNet, GLMBoost, and Random Survival Forest 3 . Through cross-validation, bootstrapping, ensemble voting, and recursive feature elimination (RFE), it enhances robustness and minimizes signature panel size without compromising performance. Designed for adaptability beyond TCGA, PRISM offers a complete pipeline from data retrieval to functional analysis and visualization. It improves survival prediction, aiding patient stratification, personalized treatment, and precision medicine. To demonstrate its utility, PRISM was applied to women-related cancers from TCGA—BRCA, OV, CESC, and UCEC—major contributors to female cancer burden 20 , 21 . Despite distinct molecular profiles, these cancers share pathways influencing progression and therapy response 22 , 23 . Comparative cancer analyses further reveal shared oncogenic drivers and therapeutic targets. While applied here to women-related cancers, PRISM is a general-purpose framework extendable to any cancer type, advancing precision oncology and improving patient outcomes 24 .

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1 Supplementary Material 2 Supplementary Material 2