Machine Learning Model on Multi-Omics Data Enables Risk Stratification and Identifies Molecular Heterogeneity and Therapeutic Targets in Glioblastoma

Machine Learning Model on Multi-Omics Data Enables Risk Stratification and Identifies Molecular Heterogeneity and Therapeutic Targets in Glioblastoma | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Machine Learning Model on Multi-Omics Data Enables Risk Stratification and Identifies Molecular Heterogeneity and Therapeutic Targets in Glioblastoma Zhenyu Zhang, Zilong Wang, Ran Li, Dongling Pei, Jingdian Liu, and 30 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7510857/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 07 Mar, 2026 Read the published version in Molecular Cancer → Version 1 posted 9 You are reading this latest preprint version Abstract Multimodal data integration reveals causal features often missed by single-modality analyses, offering a more comprehensive view of glioblastoma (GBM) complexity. We collected radiomic, pathomic, genomic, transcriptomic, and proteomic data from patients with IDH-wild-type GBM to construct a machine learning–based risk stratification model. While sample sizes varied across modalities, 147 patients with complete data across all five omics layers were used for integrative analysis. This approach identified two clinically distinct subgroups. The low-risk group, linked to favorable outcomes, showed enhanced neurodevelopmental signatures, increased neuronal infiltration, and more oligodendrocytes. In contrast, the high-risk group, associated with poor prognosis, exhibited strong proliferative signals and hyperactive cell cycle pathways. Downstream multi-omics analysis identified PDIA4, EIF3I, and RFT1 as potential prognostic biomarkers and therapeutic targets in high-risk GBM. These findings underscore the utility of multimodal machine learning in refining prognostic models, characterizing tumor heterogeneity, and informing personalized treatment strategies. Multi-Omics Risk Stratification Molecular Heterogeneity Glioblastoma Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 1. Introduction Adult-type diffuse glioma is the most common primary malignant tumor of the central nervous system( 1 ), with IDH mutation serving as a key molecular marker that stratifies gliomas into two prognostically and genetically distinct groups( 2 ). The latest WHO classification recognizes IDH-mutant gliomas and IDH wild-type GBM as separate tumor types( 3 ). Standard GBM treatment includes surgical resection followed by chemoradiotherapy and temozolomide( 4 ), yet the disease's profound intertumoral heterogeneity-spanning the genome, transcriptome, and proteome—poses significant therapeutic challenges( 5 , 6 ). Thus, a comprehensive investigation of IDH wild-type GBM heterogeneity is crucial. High-throughput sequencing and multi-omics analyses have significantly advanced our understanding of the molecular characteristics of GBM( 6 – 8 ). Verhaak et al. identified four transcriptomic subtypes, each linked to distinct genetic alterations( 9 ), while Ceccarelli et al. used DNA methylation and gene expression profiling to further refine glioma subtypes( 10 ). Other studies have classified GBM based on pharmacological response( 11 ), tumor microenvironment (TME) composition( 12 , 13 ), and radiogenomic features, revealing distinct biological and immune landscapes( 14 ). However, most existing classification models rely on single-omics data, limiting their ability to fully capture GBM complexity. With advancements in sequencing and intelligent imaging analysis, multimodal data integration for precision oncology has emerged as a frontier in cancer research( 15 ). Each modality offers unique insights into tumor biology, and their integration provides a more comprehensive understanding of tumor heterogeneity( 16 ). In this context, artificial intelligence has transformed cancer diagnostics, classification, molecular profiling, and treatment selection( 17 ). This study leverages a late-fusion strategy combined with multiple machine learning algorithms to integrate multimodal data—including radiomics, pathology, genomics, transcriptomics, and proteomics—for the risk stratification of glioblastoma (GBM) patients. This stratification framework enables a comprehensive characterization of GBM heterogeneity across bulk multi-omics, single-cell transcriptomics, and spatial transcriptomics, while also identifying potential therapeutic targets for phenotypic validation. 2. Results 2.1 Multimodal Data Summary of GBM in this study All training datasets used in this study were obtained from the First Affiliated Hospital of Zhengzhou University. A total of 1,187 adult patients diagnosed with IDH wild-type GBM were recruited. Among them, 185 patients underwent RNA sequencing (RNA-seq), 166 patients were analyzed by mass spectrometry (MS) for proteomics research, and 167 patients were subjected to whole-exon sequencing (WES). Meanwhile, 936 patients had preoperative multiparametric magnetic resonance imaging (MRI) scans, including T1-weighted imaging (T1WI), contrast-enhanced T1-weighted imaging (CE-T1WI), T2-weighted imaging (T2WI), fluid - attenuated inversion recovery (FLAIR), and apparent diffusion coefficient (ADC). Image segmentation and feature extraction were performed on the MRI data of these 936 patients. Histological whole-slide images (WSIs) were obtained from 1122 patients ( Figure 1A ). Additionally, single-cell transcriptome (scRNA-seq) data were available for 26 patients, and spatial transcriptome (ST) data were accessible for 10 patients. For a comprehensive overview of all FAHZZU multimodal data, please refer to Figure 1B and Data S1 . The external validation datasets in this study were obtained from multiple institutions. WES data from 512 patients were sourced from The Cancer Genome Atlas‌ (TCGA). RNA-seq data included 130 cases from TCGA, 72 from Chinese Glioma Genome Atlas (CGGA) -325, and 107 from CGGA-693. MRI data were collected from 108 patients at Henan Provincial People's Hospital (HPPH) and 28 from TCGA. WSI data included 125 cases from HPPH and 144 from TCGA. MS data were obtained from 31 patients at Samsung Medical Center (SMC) and 91 from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database. 2.2 Development and Validation of the Multimodal Prognostic Model During model development, we implemented a late fusion strategy ( Figure 1A )( 18 ). Initially, the optimal combination of machine learning algorithms was identified at each unimodal level to construct prognostic models for patient risk scoring ( Figure S1-S5 ). These five unimodal models were then integrated using a Cox proportional hazards regression model with 10-fold cross-validation, resulting in a comprehensive multimodal fusion-based risk stratification model. Comparative analysis demonstrated that the multimodal model significantly outperformed single-omics models across multiple evaluation metrics. It achieved superior risk stratification, with a highly significant survival difference between high- and low-risk groups ( Figure 1C ). Decision Curve Analysis (DCA) confirmed greater net clinical benefits ( Figure 1D ), while the Time-dependent Concordance Index (Time-Cindex) and overall Concordance Index (C-index) consistently exceeded those of single-omics models, indicating enhanced predictive accuracy ( Figure 1E ). The model also exhibited lower Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values, reflecting better fitness and parsimony, along with reduced Integrated Absolute Error (IAE) and Integrated Squared Error (ISE), highlighting improved calibration and reliability ( Figure 1F ). Further, the multimodal model showed significant improvements in Integrated Discrimination Improvement (IDI) and Net Reclassification Improvement (NRI), confirming its superior stratification capability ( Figure 1G-H ). A comparative evaluation of four-omics versus five-omics integration models revealed the five-omics model as the superior approach, demonstrating substantial advantages across multiple dimensions ( Figure 1I ). 2.3 Genomic Alteration Profiles of High- and Low-Risk Group To assess the impact of gene mutations on different risk groups, we analyzed key mutational features, including mutation types, genetic variation patterns, single nucleotide substitution profiles, and tumor mutational burden (TMB) ( Figure S6A ). Splice site mutations (P = 0.008) and C>A substitutions (P = 0.041) were significantly enriched in the high-risk group, though their overall proportions were low ( Figure S6B ). TMB showed no significant differences between risk groups and had a weak correlation with risk scores, suggesting a limited role in tumor progression ( Figure S6C-D ). Comparing the top 20 most frequently mutated genes in both groups revealed six shared mutations (TNN, TP53, EGFR, MUC16, NF1, PTEN), with TP53 showing a 9% higher mutation frequency in the high-risk group ( Figure 2A-B ). Functional analysis indicated that high-risk group mutations were enriched in cell cycle regulation, gene expression, and immune activity, while low-risk group mutations were associated with neuronal function, tumor stroma interactions, and substance transport ( Figure 2G ). OncodriveCLUST analysis identified risk-group-specific driver genes: high-risk drivers were linked to cell cycle regulation and immune activity, whereas low-risk drivers were associated with DNA repair and metabolic pathways ( Figure 2D, 2H )( 19 ). Co-mutation analysis further highlighted distinct biological processes, with high-risk group co-mutations enriched in gene expression regulation and disease progression, and low-risk group co-mutations associated with metabolism, neuronal function, and tumor stroma interactions ( Figure 2E, 2I ). Pathway-level mutation analysis using KEGG( 20 ), Hallmark( 21 ), and Reactome( 22 ) databases revealed distinct functional alterations: high-risk group mutations affected cellular stimulation responses, material metabolism, and gene expression, whereas low-risk group mutations impacted immune activity, coagulation, and signal transduction ( Figure 2F, 2J ). Given the impact of copy number variations (CNVs) on gene expression, we examined CNV frequencies ( Figure 3A ). GISTIC 2.0( 23 ) identified 16 chromosomal cytobands with significant CNV differences, categorized into amplifications and deletions specific to each risk group ( Figure 3B-D ). Pathway enrichment analysis of CNV-altered genes revealed that high-risk amplifications were associated with tumor signaling, cell cycle regulation, and metabolic dysregulation, while deletions affected DNA repair and apoptosis. Low-risk deletions were linked to abnormal glycosylation and oncogenic pathways (EPH-Ephrin, Rho GTPase, TGF-β), whereas low-risk amplifications showed no significant pathway enrichment ( Figure 3E-F ). Transcriptome-Proteome Heterogeneity in GBM To investigate GBM molecular heterogeneity at both transcriptomic and proteomic levels, we performed differential expression analyses to identify key altered genes and proteins between risk groups. Pathway enrichment clustering revealed distinct core biological processes: in the high-risk group, pathways related to GPCR signaling, immune responses, extracellular matrix remodeling, RNA splicing, protein synthesis, intracellular transport, and metabolism were enriched ( Figure 4A ). The low-risk group showed enrichment in pathways associated with GPCR signaling, neuronal processes, calcium ion signaling, ion transport, and neuron-tumor communication ( Figure 4B ). Gene Set Enrichment Analysis (GSEA) further validated these findings, highlighting consistent transcriptomic and proteomic enrichment patterns. In the high-risk group, dominant pathways included cell cycle regulation, gene expression, DNA repair, metabolism, and immune responses ( Figure 4C-D ), while the low-risk group enriched in synaptic activity, signal transduction, lipid metabolism, and substance transport ( Figure 4E-F ). Gene Set Variation Analysis (GSVA) confirmed the robustness of these results, further supporting pathway activity consistency ( Figure S7B ). 2.4 Immune Microenvironment Differences Between Risk Groups To explore immune microenvironment disparities between the high- and low-risk groups, we conducted GSEA on immune-related pathways, revealing heightened immune response in the high-risk group ( Figure S8A ). However, ESTIMATE analysis showed no significant differences in stromal score, immune score, or tumor purity between the groups ( Figure S8B ). Immune cell Proportion Score (IPS) analysis indicated higher effector and suppressor cell scores in the low-risk group, suggesting increased activity of these immune cells ( Figure S8C ). Tumor Immune Dysfunction and Exclusion (TIDE) and Cancer Immunity Cycle (CIC) analyses revealed contrasting immune features: the high-risk group exhibited significant Cancer-associated fibroblasts (CAF) infiltration, suggesting a strong immunosuppressive response, while the low-risk group showed greater T-cell dysfunction ( Figure S8D-E ). CIC analysis further highlighted immune differences, with higher macrophage and neutrophil recruitment in the high-risk group at Step 4 and more pronounced T-cell infiltration in the low-risk group at Step 5 ( Figure S8F-G ). Tumor Immune Prediction (TIP) analysis revealed no significant differences in overall immune activity scores between groups ( Figure S8H ). In summary, the high-risk group demonstrated elevated immune pathway activity but strong immunosuppressive CAF features, while the low-risk group exhibited higher effector cell activity but greater T-cell dysfunction. These distinct immune profiles may influence patient prognosis and therapeutic outcomes. 2.5 Integrated Biological Pathway Analysis Across Multi-Omics Levels To investigate the heterogeneity of GBM across different multi-omics levels, we performed differential analysis, Over-Representation Analysis (ORA), and GSVA for each omics dataset. The most significantly altered biological pathways at different omics levels were visualized ( Figure 5A ). The results revealed consistency in the heterogeneity patterns across multi-omics levels. In the high-risk group, pathways related to cell proliferation, cell cycle regulation, immune responses, and oncogenic signaling were highly active, whereas in the low-risk group, pathways associated with neuronal activity and molecular transport exhibited greater enrichment. Furthermore, we compared our risk classification with previously published GBM subtype classifications ( Figure 5B-C )( 8, 9, 11, 12 ). The results demonstrated that our risk stratification differed considerably from prior transcriptomic subtypes of GBM. Notably, each of our risk groups encompassed multiple transcriptional subtypes from different studies, suggesting that the multimodal risk stratification captures a broader spectrum of GBM heterogeneity beyond traditional transcriptomic subtyping. 2.6 Potential Therapeutic Targets Identification and Experimental Validation Based on risk stratification using the multimodal fusion model, we screened highly expressed genes in the high-risk group across multiple omics levels and multi-omics databases as potential therapeutic targets. Figure 6A illustrates the target screening workflow. The results of each screening step are shown in Figure S9 , ultimately identifying 18 target candidates, which were ranked based on their prognostic significance across multiple omics layers and databases. Following a review of previously published literature and preliminary experiments, PDIA4, RFT1, and EIF3I were selected as experimental targets. First, we quantified the expression levels of EIF3I, PDIA4, and RFT1 in different GBM cell lines by RT-qPCR. Cell lines with higher baseline expression levels (A172/U251 for EIF3I, A172/U118 for PDIA4, and A172 for RFT1) were selected for functional validation. Biological effects of the three targets were evaluated using EdU and CCK-8 proliferation assays, wound-healing assays, and Transwell invasion assays ( Figure 6B/D/F/H/J ). The corresponding statistical results are presented in Figure 6C/E/G/I/K . The proliferation assay results demonstrated that the silencing of any of PDIA4, RFT1, or EIF3I significantly inhibited the proliferation of GBM cells. Subsequent wound healing assays further confirmed that the knockdown of PDIA4, RFT1, or EIF3I markedly suppressed the migratory abilities of GBM cells. Additionally, the Transwell assays revealed that the PDIA4, RFT1, or EIF3I deficiency inhibited the invasive capacity of GBM cells. Collectively, these findings suggest that PDIA4, RFT1, and EIF3I play pivotal roles in regulation of GBM cell proliferation, migration, and invasion. 2.7 Cell-Type–Specific Risk Stratification and Therapeutic Target Mapping in Single-Cell Resolution We performed quality control (QC) on the single cell sequencing data and clustered and annotated 12 different cell types ( Figure 7A and Figure 8B ). We then generated a quasi-bulk mRNA profile for each patient and used a multi-algorithm classifier to integrate multimodal risk stratification to divide patients into high-risk and low-risk groups, allowing comparison of the composition of cell subpopulations ( Figure 7B ). The results of subpopulation ratio difference analysis showed that oligodendrocyte subpopulation was significantly enriched in low-risk patients. No significant enrichment of cell type subpopulations was found in the high-risk group ( Figure 7C ). Gene differential expression and pathway enrichment analyses revealed profound biological distinctions between risk groups ( Figure 7D-E ). Across a spectrum of cell type–specific subpopulations in GBM, consistent functional divergence emerged between high- and low-risk groups, underscoring distinct transcriptional programs driving prognostic heterogeneity. In the AC-like, MES-like, and OPC-like subgroups, high-risk cells were characterized by elevated oxidative phosphorylation and enhanced mitochondrial activity, indicative of metabolic reprogramming that supports cellular proliferation. In contrast, low-risk counterparts were enriched for pathways associated with signaling regulation, inflammatory equilibrium, and homeostatic maintenance. Similarly, G2M-like and NPC-like cells in high-risk patients exhibited increased protein synthesis and nutrient stress responses, whereas low-risk groups favored DNA repair, replication fidelity, and genomic stability—features associated with restrained proliferation and improved prognosis. In Oligodendrocyte-like cells, high-risk groups showed upregulated translation and ribosomal function, while low-risk cells preferentially activated heat shock, MAPK, and TNF signaling pathways, reflecting a shift from synthetic to stress-adaptive states. Endothelial and Mural subgroups in the high-risk cohort displayed enhanced ECM remodeling, epithelial–mesenchymal transition (EMT), and interferon-driven inflammatory signaling, whereas low-risk cells maintained vascular integrity, solute transport, and transcriptional balance. Immune lineages mirrored these patterns: high-risk T cells exhibited oxidative metabolism and exhaustion-like phenotypes, while low-risk T cells-maintained cytokine signaling and immune competence. B cells in high-risk tumors displayed interferon-stimulated transcription and translational stress, contrasting with antigen sensing and chromatin-modifying signatures observed in the low-risk group. In the TAM_MG and Mono subpopulations, high-risk states were defined by mitochondrial respiration and immunosuppressive features, while low-risk cells were enriched for innate immune activation via TLR, NOD, NF-κB, and IL-10 signaling pathways, reflecting preserved immunoregulatory capacity. Together, these findings highlight cell-type–specific metabolic and immunological reprogramming as fundamental to GBM risk stratification and point toward risk-specific therapeutic vulnerabilities, particularly in metabolism- and immunity-targeted strategies. To investigate the relationship between single-cell transcriptional states and patient prognosis, we performed Scissor analysis by integrating single-cell RNA-seq profiles with bulk RNA-seq–derived survival information. This approach identified two prognostically relevant cellular compartments: Scissor+ cells, which were positively correlated with poor survival, and Scissor– cells, associated with favorable outcomes ( Figure 8A ). Notably, scissor+ cells were predominantly composed of malignant populations, including AC-like, MES-like, and G2M-like subtypes, suggesting that transcriptional programs within these tumor cell states contribute directly to adverse prognosis. In contrast, scissor– cells were primarily of immune origin, particularly TAM_MG, T cells, and Monocytes, indicating that these immune populations may be associated with better clinical outcomes. These results highlight the prognostic relevance of intra-tumoral cellular heterogeneity and emphasize that tumor-intrinsic transcriptional programs, rather than immune compartment dysfunction, are the primary drivers of poor prognosis in this cohort. We further investigated the single-cell expression patterns of three candidate therapeutic targets—PDIA4, RFT1, and EIF3I—previously identified through bulk multi-omics analysis ( Figure 8C ). Mapping their expression across single-cell subpopulations revealed that all three targets were predominantly expressed in tumor-associated cell states, rather than in immune or stromal compartments ( Figure 8D ). Specifically, PDIA4 and RFT1 were highly enriched in the G2M-like malignant subpopulation, consistent with their association with proliferative programs. In contrast, EIF3I expression was primarily restricted to the MES-like subpopulation, which has been linked to aggressive and therapy-resistant phenotypes. These results suggest that the identified targets are tightly linked to specific tumor cell states and may offer subtype-selective therapeutic opportunities. 2.8 Dissecting Spatial Heterogeneity and Target Localization in GBM via Spatial Transcriptomics To further investigate the spatial organization of intratumoral heterogeneity, we performed spatial transcriptomic deconvolution using the RCTD algorithm, guided by an annotated single-cell reference atlas ( Figure 9A ). This analysis revealed that the majority of spatial spots across tumor slices were predominantly composed of MES-like malignant cells, followed by TAM-MG (tumor-associated microglia/macrophages) and AC-like tumor cells. Other subpopulations, including NPC-like, OPC-like, and immune cells such as T and B cells, were sparsely distributed. This spatial enrichment pattern suggests that MES-like cells represent the dominant malignant state within the sampled (central) tumor regions, consistent with their known invasive and therapy-resistant characteristics. We next performed a comprehensive niche analysis across all spatial spots to identify regionally coherent multicellular microenvironments. Based on compositional similarities, we defined seven distinct spatial niches ( Figure 10A ). Among them (Figure 9B-C) , niche3 and niche4 exhibited highly similar profiles characterized by enrichment of immune-related populations, including B cells, TAM-MG, T cells, and Endothelial cells, indicative of immune-dominant microenvironments with vascular involvement. In contrast, niche1 and niche5 were enriched for G2M-like, NPC-like, OPC-like, and Oligodendrocyte subtypes, reflecting proliferative and neural progenitor–like tumor ecosystems with minimal immune infiltration, possibly representing stem-like or immunologically "cold" tumor niches. niche2 was defined by the high abundance of MES-like cells with moderate immune infiltration, suggesting a transitional state between immune-active and tumor-dominant environments. niche6 represented a highly tumor-centric ecosystem dominated by G2M-like and MES-like subtypes, with negligible contributions from other cell types. Finally, niche7 was characterized by enrichment of AC-like and OPC-like cells, potentially indicating a tumor region associated with glial differentiation programs. We further visualized the spatial distribution of high- and low-risk regions across these niches ( Figure 10B ), which revealed substantial heterogeneity in risk stratification within tumor tissues. Patient-level analyses of niche composition ( Figures 10C–D ) demonstrated inter-individual differences in niche occupancy and subtype prevalence. In addition, spatial overlap between cellular subpopulations and niche types was assessed ( Figure 10E ), confirming niche-specific enrichment patterns of distinct cell types. Lastly, we examined the spatial expression patterns of three candidate therapeutic targets—PDIA4, RFT1, and EIF3I—originally identified from bulk multi-omics modeling. Consistent with scRNA-seq–based findings, all three targets displayed diffuse and broadly distributed expression across tumor-dominant niches in the spatial transcriptomic data, without evident spatial restriction or compartmentalization ( Figures 10F–G ). This pervasive expression pattern underscores their potential as pan-tumoral therapeutic targets within high-risk GBM regions. Collectively, these results demonstrate that GBM tissues harbor spatially segregated multicellular ecosystems with distinct tumor–immune compositions, and that therapeutic targets identified through bulk profiling remain spatially relevant within tumor-dominant ecological contexts. 3. Discussion This study integrates multimodal data from GBM patients to develop a robust risk stratification model, enhancing the understanding of GBM heterogeneity across different risk groups. Key findings of our study include: ( 1 ) the development of a five-modal fusion risk prediction model that outperforms any single-modality or four-modality model, effectively stratifying patients into distinct risk groups with significant survival differences, offering valuable insights for personalized treatment strategies; ( 2 ) a comprehensive multi-omics analysis of GBM across different resolution levels, providing an in-depth characterization of GBM heterogeneity and identifying potential therapeutic targets, thereby presenting promising prospects for precision-targeted therapy. Several studies have explored GBM risk stratification, but most rely on single-omics or clinical data. While efforts like Yan et al( 24 )., Philipp et al( 25 ). and Junseong et al( 26 ). integrated clinical and radiomic features or transcriptomic data, they lacked multimodal integration. Some studies have utilized multimodal data, such as Wang et al( 6 )., who classified GBM subtypes using proteogenomic and metabolomic data, and Ravi et al. ( 27 ), who examined tumor-host interactions through spatially resolved multimodal analysis. However, their clinical prognostic value remains limited. In contrast, our study integrates genomic, transcriptomic, proteomic, radiomic and pathomic data using a late-stage fusion strategy, enhancing prognostic accuracy and providing novel insights into tumor heterogeneity and potential therapeutic targets. Genomic analysis of the high-risk group revealed specific mutations, including TP53, PI3KR1, ARHGEF18, and NEB, with no significant differences in mutation types or tumor mutation burden compared to the low-risk group. Further enrichment analysis highlighted the impact of pathways related to cell cycle regulation and gene expression. Copy number variation (CNV) analysis corroborated these findings, showing amplifications in genes linked to cell cycle regulation and metabolic dysregulation. These genomic alterations were consistent with transcriptomic and proteomic data, reinforcing the value of multi-omics integration for understanding the molecular heterogeneity of GBM. The immune microenvironment in the high-risk group was marked by increased immune pathway activation but a predominantly immunosuppressive state. This was characterized by extensive cancer-associated fibroblast (CAF) infiltration and recruitment of inflammatory cells (macrophages and neutrophils), suggesting immune tolerance mechanisms. These findings align with previous studies on the role of CAFs in immune suppression and poor prognosis( 28 ). Single-cell transcriptomic analysis further revealed the enrichment of AC-like and NPC2-like cells, which are linked to immunosuppression and tumor progression( 29 ). Pathway enrichment analysis highlighted the involvement of EGFR/ERBB2 signaling and P53 pathway abnormalities, contributing to immune evasion and tumor progression. In contrast, the low-risk group exhibited mutations in EGFR, PARP4, and CHRNB3, with enrichment in neuronal activity and molecular transport pathways. CNV analysis indicated a more stable genomic profile with fewer alterations in genes related to cell cycle regulation and gene expression control. The immune microenvironment in the low-risk group displayed severe T-cell dysfunction, but immune cell infiltration, particularly T-cells, was more pronounced than in the high-risk group( 30 ). By combining quasi-bulk transcriptomes with multimodal risk models, we defined patient-specific risk groups and linked them to distinct cellular compositions. Notably, low-risk patients exhibited significant enrichment of oligodendrocyte-like subpopulations, whereas high-risk patients lacked any dominant non-malignant cell type, suggesting that preservation of glial-like differentiation may be associated with more favorable outcomes. This observation echoes prior reports highlighting the protective role of oligodendrocytic features in gliomas and reinforces the relevance of cell-type composition in prognostic evaluation( 31 ). Gene expression and pathway enrichment analyses across malignant and non-malignant subpopulations revealed risk-specific transcriptional programs. High-risk cells exhibited elevated oxidative phosphorylation, translational activity, and nutrient stress pathways, particularly within G2M-like, AC-like, and MES-like tumor compartments. These features reflect metabolic reprogramming and proliferative signaling that may drive tumor aggressiveness. In contrast, low-risk counterparts preferentially activated homeostatic pathways, including stress response, immune signaling, and DNA repair. Immune subpopulations also mirrored this divergence: high-risk T cells and B cells exhibited exhaustion-like features and interferon stress, while low-risk immune cells retained cytokine signaling and immunocompetence. These results support the notion that risk stratification in GBM reflects a convergence of tumor-intrinsic metabolic pressure and immune dysfunction. To further link cell states with patient survival, we applied Scissor analysis, revealing that high-risk–associated (Scissor⁺) cells were predominantly malignant, whereas low-risk–associated (Scissor⁻) cells were mainly immune-derived. This reinforces the idea that tumor-intrinsic programs are the major drivers of poor prognosis, while immune competence may be associated with improved outcomes. Importantly, Scissor helped resolve the functional relevance of transcriptionally defined subpopulations, providing a bridge between bulk survival signals and cellular identities. We also assessed the expression patterns of three candidate therapeutic targets—PDIA4, RFT1, and EIF3I—identified via bulk multi-omics modeling. Single-cell and spatial transcriptomic analysis revealed that all three targets were highly expressed within tumor-dominant cell states, particularly G2M-like and MES-like subpopulations, with broad and diffuse distribution across tumor regions. These findings suggest that the proposed targets are not spatially restricted and may be amenable to pan-tumoral therapeutic strategies in high-risk GBM. Finally, spatial deconvolution and niche analysis highlighted regionally coherent multicellular ecosystems with distinct tumor–immune compositions. Among the seven identified spatial niches, some were immune-dominant (e.g., niche3 and niche4), others were stem-like or proliferation-driven (e.g., niche1 and niche5), and some were highly tumor-centric (e.g., niche6 and niche7). Risk stratification and target gene expression were non-uniformly distributed across these niches, further reinforcing the spatial dimension of prognostic heterogeneity. Despite promising results, this study has limitations: ( 1 ) While external validation was performed for the single-omics prognostic models, the lack of publicly available multi-omics datasets precluded external validation of the multimodal fusion risk stratification model at total level. ( 2 ) During the development and validation of single-omics models, some models demonstrated high prognostic value in datasets derived from Asian populations but failed to exhibit sufficient prognostic significance in cohorts from TCGA, which is comprised of mainly non-Asian populations. ( 3 ) The classifier used for risk mapping in single-cell and spatial transcriptomics analyses was trained on transcriptomic data, which may limit its ability to capture the full multimodal complexity of GBM. Future studies are needed to integrate more comprehensive multimodal datasets such as several external validation sets, and prospective validation sets to further corroborate the robustness of the model. 4. Materials and Methods 4.1 Experimental Design We designed a multi-step integrative framework to investigate prognostic heterogeneity in glioblastoma (GBM) and identify potential therapeutic vulnerabilities. We collected a comprehensive multi-omics dataset for GBM patients, including genomic, transcriptomic, proteomic, radiomic, and histopathologic data, all matched with survival information. For each omics modality, multiple machine learning–based survival algorithms were applied to construct individual prognostic models, generating modality-specific risk scores for each patient. These risk scores were subsequently integrated into a composite risk score using a Cox proportional hazards model–based fusion strategy. Patients were then stratified into high-risk and low-risk groups according to the optimal cutoff of the composite risk score. To elucidate molecular differences between risk groups, we performed systematic heterogeneity analyses across multiple omics layers, including gene/protein expression profiles, mutation burden, pathway activities, and tumor microenvironment composition. Single-cell and spatial transcriptomic data were further incorporated to dissect the contributions of specific cell types and spatial niches to prognostic stratification. Based on the results of differential analyses, we prioritized candidate therapeutic targets that were consistently overregulated in high-risk patients, strongly associated with poor survival, and possessed druggable potential. These targets were then validated for cell-type specificity and spatial distribution at single-cell and spatial resolution, enhancing their translational applicability. 4.2 Data and Sample Collection This retrospective study collected data from patients with IDH-wildtype GBM who underwent surgical resection at FAHZZU between 2015 and 2021, in accordance with the 2021 WHO Classification of CNS Tumors. Inclusion criteria were: ( 1 ) age ≥ 18 years, ( 2 ) primary diffuse glioma, ( 3 ) IDH-wildtype GBM diagnosis based on the 2021 WHO classification, ( 4 ) no prior radiotherapy or chemotherapy, and ( 5 ) complete clinical and follow-up data. Exclusion criteria included: ( 1 ) history of brain surgery or head trauma and ( 2 ) prior radiotherapy or chemotherapy. A total of 1,187 patients with IDH-wildtype GBM were included (Fig. 1 B). Of these, 1,122 had postoperative pathological specimens, from WSIs of H&E-stained sections were obtained. Preoperative MRI data were available for 936 patients, including T1-weighted, contrast-enhanced T1-weighted, T2-weighted, FLAIR, and ADC maps, all of high quality. Fresh surgical tumor specimens from 211 patients were snap-frozen and stored for tissue sequencing. Of these, 185 underwent RNA sequencing (RNA-seq), 167 underwent WES, 166 had proteomic data, 26 had scRNA-seq data, and 10 underwent spatial transcriptomic analysis. Preoperative MRI scans and postoperative slides were also collected from 108 and 125 GBM patients, respectively, at HPPH. The study was approved by the Human Research Ethics Committees of Henan Provincial People's Hospital (Approval No. 2023 − 174) and FAHZZU (Approval Nos. 2019-KY-176 and 2023-KY-1028), with informed consent obtained from all patients for the use of fresh tumor specimens. 4.3 Collection and Processing of GBM Data from Public Databases In this study, we aimed to collect as much sequencing data on GBM (GBM) as possible from public databases to validate our conclusions and enrich the research content. We obtained five GBM datasets from the following sources: The Cancer Genome Atlas (TCGA; https://www.cancer.gov/ccg/research/genome-sequencing/tcga ), the Chinese Glioma Genome Atlas (CGGA; http://www.cgga.org.cn/ ), and the Clinical Proteomic Tumor Analysis Consortium (CPTAC; https://proteomics.cancer.gov/programs/cptac ). These datasets included TCGA-GBM, CGGA325, CGGA693, and CPTAC-GBM. Additionally, we collected an external proteomics dataset, SMC-GBM, containing proteomic profiles of GBM patients from a previously published study( 32 ). 4.4 MRI Scanning and Imaging Feature Extraction Patient MRI images were acquired during routine examination using a 3.0 T MRI scanner (Siemens Magnetom Skyra/Trio TIM; GE Discovery MR750; Philips Ingenia). Sequences included: axial and sagittal T1-weighted imaging (T1WI), axial T2-weighted imaging (T2WI), axial T2-weighted fluid-attenuated inversion recovery (FLAIR) imaging, as well as axial, sagittal, and coronal post-contrast T1-weighted imaging (CE-T1WI) immediately after intravenous injection of a 0.1 mmol/kg dose of gadolinium-based contrast agent. Apparent diffusion coefficient (ADC) maps were obtained from axial diffusion-weighted imaging (DWI). The acquisition parameters for each sequence were as follows: T1WI and CE-T1WI: repetition time (TR) 220–1750 ms; echo time (TE) 2.3–24 ms; echo train length (ETL) 1–12; slice thickness 5 mm; averages/excitations 1; flip angle (FA) 70°-111°; field of view (FOV) 220×192–240×240 mm²; matrix 256×162–320×256 mm². T2WI: TR 1873–5390 ms; TE 70–117 ms; ETL 16–32; slice thickness 5 mm; averages/excitations 1; FA 90°-142°; FOV 220×192–240×240 mm²; matrix 320×238–512×512 mm². FLAIR: TR 4500–8400 ms; TE 85–150 ms; inversion time (TI) 1670–2250 ms; ETL 1–38; slice thickness 5 mm; averages/excitations 1; FA 90°-150°; FOV 220×192–240×240 mm²; matrix 256×179–256×256 mm². DWI: Images were processed by the corresponding post-processing workstation, and ADC images were calculated from DWI acquired at b-values of 0 and 1000 s/mm². Sequence parameters included: TR 2121–6000 ms; TE 77–119 ms; ETL 1–82; slice thickness 5 mm; averages/excitations 1; FA 90°; FOV 220×220–240×240 mm²; matrix 152×114–192×192 mm². ADC maps for all imaging planes were generated on a voxel-by-voxel basis using a single-exponential model. First, the N4ITK algorithm was employed to correct bias field distortions for all sequences. After isotropic voxel resampling to 1×1×1 mm³ through trilinear interpolation, multi-sequence MRI rigid registration for each patient was performed using the axial resampled CE-T1WI as a template, and mutual information as similarity measure. This process was completed using the 3D Slicer software, generating registered images rT1WI, rCE-T1WI, rT2WI, rFLAIR, and rADC. Histogram matching was used for gray-level normalization on rT1WI, rCE-T1WI, rT2WI, and rFLAIR. We set the histogram level to 1024 and the number of matching points to 10 to achieve a finer match while preserving more details. A deputy chief physician in neuroradiology with over 10 years of experience in head MRI diagnosis manually delineated the tumor region of interest (ROI) on the axial plane of rFLAIR, rT2WI, and rCE-T1WI images using ITK-SNAP software, obtaining the tumor volume of interest (VOI)( 33 ). The VOI was defined as the enhanced area, non-enhanced area, and necrotic area of the tumor. The VOI contour was drawn based on FLAIR images, while rT2WI and rCE-T1WI were used for cross-checking the tumor extent and fine-tuning the tumor contour. Z-score normalization was applied within the VOI for all sequences to adjust the ROI intensity to have a mean of 0 and a standard deviation of 1. This radiologist and a deputy chief physician in neurosurgery with over 10 years of work experience randomly selected 100 patients within the group for VOI redrawing using a simple random sampling method. Interclass correlation coefficients (ICC) were used to evaluate intra-rater reliability analysis for the test-retest dataset and inter-rater reliability analysis for the multiple description dataset, retaining features with ICC ≥ 0.75. The obtained VOI was then overlaid with co-registered rT1WI, rCE-T1WI, rT2WI, rFLAIR, and rADC images. PyRadiomics was used to extract three categories of features, including first-order intensity statistics, shape descriptors, and higher-order texture features( 34 ). Five basic matrices were employed to define texture features: the gray-level co-occurrence matrix (GLCM), gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), gray-level dependence matrix (GLDM), and neighborhood gray-tone difference matrix (NGTDM). In this study, imaging features were extracted from three types of images: original images, wavelet images, and Gaussian Laplace images. A PyRadiomics parameters file was provided in Github repository to enhance the reproducibility of feature extraction. Ultimately, 4788 features were extracted from the five MRI sequences, retaining 3015 features with ICC ≥ 0.8. 4.5 Hematoxylin and Eosin (H&E) Histological Slide Scanning and Feature Analysis Pathology slides were scanned at 20x magnification using a digital pathology scanner (KF-PRO-120-HI) to obtain the original WSIs. Subsequently, the original WSI underwent color space conversion, tissue segmentation, patch selection, and feature extraction. Specifically, the WSI at the 5x resolution was converted from RGB to Lab color space, and Otsu's algorithm was then applied to calculate a segmentation threshold for segmenting the tissue from WSI. The obtained tissue image was tiled into many 1024×1024 patches at 20× magnifications, where these patches were adjacent to one another covering the WSI. A Python package Yottixel was used to select the optimal patches for further analysis( 35 ). Finally, CellProfiler (v4.2.5) software was used to extract features from each selected patch( 36 ). 4.6 Whole Exome Sequencing (WES) and Analysis Tumor tissue and adjacent brain tissue DNA were extracted from patients using the QIAamp Fast DNA Tissue Kit (Qiagen). Blood patients were collected in tubes containing EDTA and centrifuged at 1600 xg for 10 minutes at 4°C within 2 hours of collection. Peripheral blood lymphocyte (PBL) pellets were stored at -20°C until further use, and PBL DNA was extracted using the RelaxGene Blood DNA System (Tiangen Biotech Co., Ltd., Beijing, China). DNA quantification was performed using the Qubit 3.0 Fluorometer and Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Inc., Waltham, MA, USA). DNA collected from tissue and PBL patients were fragmented using dsDNA Fragmentase enzyme (New England BioLabs, Inc., Ipswich, MA, USA), followed by size selection of DNA fragments (150–250 bp) using Ampure XP beads (Beckman Coulter, Inc., Brea, CA, USA). The KAPA Library Preparation Kit (Kapa Biosystems, Inc., Wilmington, MA, USA) was employed for the construction of DNA fragment libraries. Cleanup steps were performed using Agencourt AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA). After DNA fragmentation, end repair and 3' A-tailing were conducted, followed by exon capture using the Agilent SureSelect Human All Exon V6 kit. The Qubit 3.0 Fluorometer and Qubit dsDNA HS Assay Kit were utilized to assess the purity and concentration of DNA fragments. Fragment length was measured using the DNA 1000 kit (Agilent Technologies, Inc., Santa Clara, CA, USA) on a 4200 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). DNA libraries with 150 bp end sequences were sequenced using the Illumina Novaseq 6000 system. Raw data were converted to FASTQ files, and adapter and low-quality reads were trimmed using Trimmomatic (v0.39). We achieved a median coverage depth of 112x for tumor specimens and 128x for non-tumor specimens. GATK (v4.2) tools were used to identify single nucleotide variants (SNVs) and insertions or deletions (INDELs). Paired-end WES reads were mapped to the human reference genome (hg38) using BWA-mem (v0.7.17). BAM files were further processed by reordering, sorting, marking duplicates, and adding read groups using Picard (v2.24.2). Base quality score recalibration was performed using the BaseRecalibrator module in GATK, followed by the assessment of cross-sample contamination using the GetPileupSummaries and CalculateContamination modules. Somatic variants were detected by MuTect2 and annotated using ANNOVAR, with patient-matched normal DNA sequencing reads serving as reference. Candidate somatic variants were distinguished based on the following filtering criteria: ① Variants outside of exonic regions and splice sites were excluded; ② Variants with a variant allele fraction (VAF) ≥ 5% and at least 2 supporting variant reads in tumor patients were retained; ③ Variants with a mutation allele frequency (MAF) ≥ 5% in at least one database, including 1000 Genomes, ESP6500, gnomAD, and ExAC, were removed. Normal patients were sequenced using the same scheme, each sample was reduced to 4% and then pooled as reference. To obtain high-quality and reliable somatic variants, we employed stringent downstream filtering criteria: ① Variants outside of exonic regions and splice sites were excluded; ② Variants with a VAF ≥ 5%, at least 5 supporting variant reads in tumor patients, and variants with a VAF in the tumor that was more than five times the VAF in the normal sample were retained; ③ Variants with more than 100 occurrences in COSMIC (v92) were retained; ④ Variants with a MAF ≥ 1% in at least one variant database (1000 Genomes, ESP6500, gnomAD, and ExAC) were removed; ⑤ Variants predicted as benign in at least two of the following tools: MutationAssessor, MutationTaster2, Polyphen2, and SIFT, were removed. Somatic CNVs were inferred by CNVkit (v0.9.9) based on BAM files generated during the somatic mutation detection process, using the default circular binary segmentation algorithm. Segment-level log2 ratios were calculated and transformed as input for the GISTIC 2.0 software to identify significantly amplified or deleted chromosomal regions in the tumors. CNV amplifications and deletions were defined using a ± 0.3 log2 ratio threshold. 4.7 RNA Sequencing (RNA-seq) and Analysis Total RNA was extracted from tissue patients using the TRIzol Reagent Kit (Ambion, Invitrogen, USA). RNA concentration and integrity were assessed using the Qubit RNA Assay Kit, Qubit 2.0 Fluorometer (Life Technologies), and Agilent 2100 Bioanalyzer (Agilent Technologies). Patients with an RNA integrity number greater than 5 were included in the study. Libraries were prepared from patients with high RNA integrity, no contaminants, and sufficient RNA quantity. Poly-T oligonucleotide magnetic beads were used to purify RNA from total RNA. RNA was fragmented in NEBNext First Strand Synthesis Reaction Buffer (5X) using divalent cations at elevated temperatures. cDNA synthesis, end repair, A-tailing, and NEBNext Adaptor ligation were performed using the NEBNext Ultra RNA Library Prep Kit. Library fragments were purified with AMPure XP (Beckman Coulter, Beverly, USA), selecting cDNA fragments of 150–200 bp in length. Library quality was assessed using the Agilent Bioanalyzer 2100. Libraries were sequenced on the Illumina HiSeq X Ten platform, generating 150 bp paired-end reads. Sequencing data were filtered using Trimmomatic software to remove adaptors and low-quality sequences, followed by data quality assessment using FastQC. STAR (v2.7.6a) was used to align sequences to the reference genome (hg38). Gene expression values were calculated using RSEM (v1.3.3) based on the GENCODE (v35) gene annotation file. HTSeq v0.6.0 was used to count the number of reads aligned to each gene, and gene expression levels were quantified as FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) and TPM (Transcripts Per Kilobase of exon model per Million mapped reads). 4.8 Mass Spectrometry Patients were removed from − 80°C storage, and an appropriate amount of tissue was weighed and placed into a liquid nitrogen pre-chilled mortar. Liquid nitrogen was added, and the tissue was thoroughly ground into a powder. Lysis buffer (1% Triton X-100, 1% protease inhibitor, 1% phosphatase inhibitor, 3 µM TSA, 50mM NAM) was added to each sample at four times the volume of the powder, followed by ultrasonic lysis. Patients were centrifuged at 4°C, 12,000 g for 10 minutes to remove cell debris, and the supernatant was transferred to a new centrifuge tube. Protein concentration was determined using a BCA assay kit. Equal amounts of protein from each sample were digested with trypsin, and the volume was adjusted using lysis buffer. One volume of pre-chilled acetone was added, vortexed, and then four volumes of pre-chilled acetone were added, followed by precipitation at -20°C for two hours. Patients were centrifuged at 4,500 g for 5 minutes, and the supernatant was discarded. The pellet was washed twice with pre-chilled acetone. After air-drying the pellet, it resuspended in 200 mM TEAB, and trypsin was added at a 1:50 ratio (protease: protein, w/w) for overnight digestion. Dithiothreitol (DTT) was added to a final concentration of 5 mM, and patients were reduced at 56°C for 30 minutes. Iodoacetamide (IAA) was added to a final concentration of 11 mM, and patients were incubated in the dark at room temperature for 15 minutes. Patients were separated using an Agilent 300Extend C18 column (4.6×250 mm), with detection at 214 nm, a column oven temperature of 35°C, and 95% buffer A for 30 minutes to equilibrate the column. After baseline stabilization, the fractionation gradient method was initiated, and peptide patients were loaded onto the high-performance liquid chromatography (HPLC) fractionation column. Patients were collected at 1-minute intervals, with fractions 11 to 46 combined into 12 groups and vacuum-dried. Peptides were dissolved in mobile phase A and separated using an EASY-nLC 1200 ultra-HPLC system. Mobile phase A consisted of a 0.1% formic acid and 2% acetonitrile aqueous solution, while mobile phase B consisted of a 0.1% formic acid and 90% acetonitrile aqueous solution. The gradient was set as follows: 0–96 min, 6%-25% B; 96–114 min, 25%-35% B; 114–117 min, 35%-80% B; and 117–1200 min, 80% B, with a flow rate maintained at 500 nl/min. Separated peptides were ionized in the NSI ion source and data were collected using the Orbitrap Exploris 480 mass spectrometer. Liquid chromatography (LC) parameters were consistent with those used during library construction. Peptides were separated using the ultra-high-performance liquid chromatography system and analyzed using the Orbitrap Exploris 480 mass spectrometer. Precursor ions and their fragment ions were detected and analyzed using the high-resolution Orbitrap. FAIMS compensation voltage (CV) settings were − 40 V, -55 V, and − 70 V. The primary mass scan range was set at 350–1350 m/z with a resolution of 120,000; the secondary scan resolution was set at 30,000. The secondary data acquisition mode was set to DIA mode, which followed a primary scan with 20 m/z window peptide ions entering the HCD collision cell using 32% collision energy for fragmentation, and subsequent secondary mass analysis. The automatic gain control (AGC) for the secondary spectrum was set at 600%. 4.9 Development and Validation of a Risk Stratification Model Based on Multimodal Data Integration In this study, the multimodal GBM dataset from the FAHZZU was used as the training set for each omics layer, while other GBM datasets served as external validation cohorts. At each omics level, we implemented a comprehensive analytical framework utilizing the R package Mime1, which integrates 10 distinct machine learning algorithms along with their 101 potential combinations( 37 ). The combination that achieved the highest average concordance index (C-index) in the external validation sets was selected as the optimal algorithm for that omics layer, which was then used to calculate patient-specific risk scores. Subsequently, utilizing a late-fusion strategy combined with 10-fold cross-validation, we integrated the risk scores from different omics layers through Cox proportional hazards regression to derive a comprehensive multimodal fusion risk score. Based on the optimal cutoff value, patients were stratified into high- and low-risk groups, thereby establishing the multimodal fusion-based risk stratification model for GBM(Figure 1 A). 4.10 Genomic Alteration Analysis Using the maftools package, WES data was processed( 38 ). To investigate the differences in gene mutation profiles between the two risk groups, we conducted a comprehensive analysis encompassing multiple mutation categories. Initially, the Wilcoxon rank-sum test was applied to assess variations in gene mutation types (frameshift, non-frameshift, point mutations, splice site mutations, and translation initiation mutations), variant classifications (deletions, insertions, multinucleotide polymorphisms, and single nucleotide polymorphisms), and single base substitution types (C > A, C > T, T > A, C > G, T > C, and T > G). Building upon these findings, we calculated the mutation frequencies of all genes across the risk groups and visualized genes with mutation frequencies exceeding 5% using maftools. To elucidate the functional implications of these mutated genes, over-representation analysis (ORA) was performed with the ClusterProfiler package( 39 ), retaining pathways with p-values < 0.05 as statistically significant. Subsequently, to identify group-specific tumor driver genes, we utilized the OncodriveCLUST algorithm, visualizing the results with maftools and performing ORA with ClusterProfiler to determine enriched pathways associated with these driver genes. In parallel, co-occurring and mutually exclusive mutation patterns between the risk groups were explored using Fisher's exact test. The resulting gene interaction networks were visualized via the GGally package, providing insights into distinct mutational landscapes, which were further analyzed through ORA to identify significant pathways (p < 0.05). To comprehensively assess the mutational impact at the pathway level, we calculated pathway mutation burdens based on gene sets from the Reactome, KEGG, and Hallmark databases using maftools. The mutation burden across pathways was visualized with ggplot2, with particular attention to those affecting more than 10% of patients in each group. Finally, pathways identified from the aforementioned analyses were systematically annotated and classified. This integrative approach enabled a detailed evaluation of how gene mutations influence key biological processes, providing mechanistic insights into the differential risk stratification between the two groups. Furthermore, to explore chromosomal-level alterations, we performed a differential analysis of copy number variations (CNVs) between the high- and low-risk groups to identify significantly altered chromosomal segments. First, significantly different CNV regions were identified based on GISTIC 2.0( 23 ) analysis results, with statistical comparisons between groups conducted using the limma package. Following the identification of these regions, functional analysis was performed to elucidate their biological relevance. Genes located within the significantly altered chromosomal segments were extracted and subjected to ORA using ClusterProfiler, taking into account both amplification and deletion statuses. Pathways with p-values < 0.05 were considered significant, enabling the identification of affected biological pathways and their potential functional consequences associated with CNV differences across risk groups. 4.11 Functional Analysis Based on Transcriptomic and Proteomic Profiles To investigate transcriptomic and proteomic differences between the high- and low-risk groups, we performed a comprehensive differential expression analysis to identify significantly altered genes and proteins. First, using the Wilcoxon rank-sum test, we applied the limma package to compare samples from different risk groups, selecting genes and proteins with a false discovery rate (FDR) 1 as significantly differentially expressed. Based on these differentially expressed genes (DEGs) and proteins, we conducted over-representation analysis (ORA) using the ClusterProfiler package to identify enriched pathways with p-values < 0.05. The enriched pathways were subsequently clustered and visualized using the aPEAR package( 40 ), allowing for a clearer understanding of functional groupings. To further explore biological processes associated with expression changes, we performed GSEA. DEGs were ranked according to their log₂FC values, and GSEA was carried out with ClusterProfiler, identifying pathways with FDR < 0.05. The results were clustered and visualized using aPEAR, providing insights into pathway-level alterations between the risk groups. In addition to ORA and GSEA, we employed gene set variation analysis (GSVA) to assess pathway activity differences across patient transcriptomes. By leveraging gene sets from the Reactome, KEGG, and Hallmark databases, GSVA scores were calculated and filtered based on reliable pathways identified through GSEA. These scores were then visualized using the ComplexHeatmap package( 41 ), facilitating the annotation and comparison of relevant pathways between the two risk groups. Through the integration of differential expression analysis, pathway enrichment, and variation analyses, this comprehensive approach provides a multidimensional understanding of the molecular differences underlying risk stratification in GBM patients. 4.12 Tumor Microenvironment (TME) Analysis Based on Transcriptomic Data Based on the stratification of high and low-risk groups, we inferred the tumor immune microenvironment from the transcriptomic data and identified immune activity-related indicators that exhibited significant differences between the groups. First, we performed a differential analysis of immune-related pathways by selecting immune pathways from the GSEA results, categorizing their functions, and visualizing them using the R package ggplot2 to examine the differences in pathway activation between the two groups. Second, we evaluated the immune and stromal cell infiltration in tumor tissues using the ESTIMATE( 42 ) algorithm and calculated the immune and stromal scores. We then conducted a Wilcoxon rank-sum test to analyze the differences in scores between the two risk groups. Third, we assessed the immune phenotypes, including antigen presentation (MHC-MHC molecules), effector cells (EC-effector cells), suppressor cells (SC-suppressor cells), and immune checkpoints (CP-immune checkpoints) from the transcriptomic data and calculated the scores for each immune activity( 43 ). A Wilcoxon rank-sum test was performed to examine the differences in immune scores between the two risk groups. Fourth, we analyzed the Tumor Immune Dysfunction and Exclusion( 44 ) by predicting the cancer patients' responses to immune checkpoint inhibitors (such as PD-1/PD-L1 inhibitors) and other immune-related parameters. Differences in immune-related indices between the two groups were assessed using the Wilcoxon rank-sum test. Finally, we performed a differential analysis of the Cancer Immune Cycle( 45 ) by predicting the anti-cancer immune states at seven immune stages and calculating the corresponding scores. The Wilcoxon rank-sum test was applied to analyze the differences in anti-cancer immune state scores between the two risk groups. 4.13 Heterogeneity Analysis Between Single-Cell and Spatial Transcriptomic Data The R package seurat( 46 ) was used to perform preprocessing such as single-cell data quality control and annotation to determine the cell subtypes contained in the data, and the RCTD algorithm( 47 ) was used to map the single-cell subtypes to the patient's spatial transcriptome data. Based on the multimodal fusion risk stratification model, single-cell and spatial transcriptomic data were categorized into risk groups, and their heterogeneity was comprehensively analyzed. First, the multimodal cohort was divided into a training set and a validation set in a 7:3 ratio. In the training set, logistic regression and ROC analysis were used to select genes with P 0.7, followed by Lasso regression for feature selection and dimensionality reduction. Genes with non-zero coefficients were retained for model construction, and a classifier was built by combining 14 machine learning algorithms (e.g., AdaBoost, decision tree, GBDT, SVM, XGBoost). The classifier was applied to pseudo-bulk data from the single-cell transcriptomic cohort to assign risk groups. Next, the Wilcoxon rank-sum test was employed to analyze differences in cell subpopulation proportions between the two groups, identifying significantly different subpopulations. Differential gene expression analysis was then performed for the cell subpopulations, selecting genes with FDR 2, followed by enrichment analysis using ClusterProfiler to identify significant pathways (P < 0.05) and analyze functional characteristics of the same subpopulations across risk groups. Additionally, we employed the R package Scissor to integrate single-cell transcriptomic data with bulk RNA-seq and survival information, enabling the identification of single-cell subpopulations associated with patient prognosis( 48 ). This approach allowed us to pinpoint individual cells whose transcriptional profiles are positively or negatively correlated with clinical outcomes. 4.14 Identification of Therapeutic Targets Based on Multi-Omics Data Based on multi-modal risk stratification, we identified EIF3I as a potential therapeutic target associated with the high-risk group. Firstly, differential analyses were conducted at the copy number variation, transcriptomic, and proteomic levels according to the multi-modal risk stratification, and the characteristic genes of the high-risk group were selected, followed by identification of the common genes across these levels. These genes were then subjected to correlation analysis between copy number variation-transcriptome and transcriptome-proteome data, which helped identify genes with high correlations across different modalities. Subsequently, these genes were analyzed in the TCGA and CGGA databases to assess differential expression between GBM and low-grade gliomas at the transcriptomic level, selecting those with high expression in GBM. Finally, survival analysis was performed on these genes across multiple databases at various omics levels and based on the Log-rank P-value and univariate Cox P-value (both < 0.05), potential therapeutic targets related to the high-risk group were screened and ranked. 4.15 Acquisition of GBM Single-cell RNA-seq Data Tissue Dissociation and Single-cell Suspension Preparation: The tissue was transferred to a culture dish containing 1× PBS (without RNase and Ca/Mg ions), placed on ice, and washed to remove surface blood clots and fats. The tissue was then minced into 0.5 mm² pieces and washed again with PBS. Dissociation solvent (0.35% collagenase IV, 2 mg/ml papain, 120 Units/ml DNase I) was added, and the mixture was incubated in a 37°C water bath shaker (100 rpm) for 20 minutes. The reaction was terminated by adding PBS containing 10% fetal bovine serum, followed by gentle pipetting 5–10 times. The cell suspension was filtered through a 70 − 30 µm cell strainer, and centrifuged at 300g for 5 minutes at 4°C to collect the cell pellet. The pellet was resuspended in 100 µl of 1× PBS (0.04% BSA) solution, followed by the addition of 1 ml of 1× red blood cell lysis buffer. The cells were incubated at room temperature or on ice for 2–10 minutes. Afterward, cells were centrifuged for 5 minutes, and the pellet was resuspended in 100 µl of Dead Cell Removal MicroBeads reagent and incubated for 15 minutes at room temperature. Dead cells were removed using MS Columns, followed by another 5-minute centrifugation to collect the viable cells, which were then resuspended in 1× PBS (0.04% BSA). After two centrifugation steps at 4°C, the final viable cells were obtained. Cell viability was assessed using trypan blue staining (> 85%) and the cell count was determined, with a concentration of 700–1200 cells/µl. Chromium 10x Genomics Library Construction and Sequencing: According to the protocol of the 10X Genomics Chromium Single-Cell 3' Kit (V3), the single-cell suspension was loaded onto the 10x Chromium chip, with an expected capture of 8,000 single cells. cDNA amplification and library construction were performed according to standard protocols, and the libraries were sequenced on an Illumina NovaSeq 6000 system with paired-end sequencing (150 bp) at a depth of at least 20,000 reads per cell. The operation was performed by LC-Bio Technology (Hangzhou, China). 4.16 Acquisition of GBM Spatial Transcriptomics Data OCT-embedded tissue blocks were sectioned into 10 µm thick slices measuring 6.5 mm x 6.5 mm, which were then adhered to Visium slides. Hematoxylin and eosin (H&E) staining was performed according to the 10x Genomics Visium Fresh Frozen Tissue Processing Protocol. Fluorescent stitching microscopy (Olympus Fluoview 1000) was used to capture H&E images. The tissue was removed, and library construction was carried out according to the 10x Genomics recommended protocol. Each spatial transcriptomics library was processed using the 10x Genomics Space Ranger software (version 2.0.0) and aligned to the reference genome (mm10 for mouse/GRCh38 for human). UMI counts were compiled for each probe site. To distinguish probe sites covering the tissue from background, only barcodes associated with the tissue-covered regions were retained. The filtered UMI count matrix was generated, followed by manual exclusion of probe sites that were not covered by tissue and undetected by Space Ranger. The UMI count matrix was further optimized. 4.17 Experimental verification of biological functions of the targets ( 1 ) Cell culture: The GBM cell lines U251(RRID:CVCL_0021), A172 (RRID:CVCL_0131) and U118(RRID:CVCL_0633) were obtained from Fenghui Biotechnology Co.Ltd. Cells were cultured in Dulbecco's modified Eagle's Medium (DMEM, Sigma) supplemented with 10% fetal bovine serum (FBS, Gibco) at 37.0°C under a humidified 5.0% CO2 atmosphere. All cell lines were tested for mycoplasma every 2 months.All cell lines were contamination free. ( 2 ) Cell transfection: All siRNAs and pLKO.1-puro shRNA plasmids (RRID:Addgene_8453) were synthesized by Tsingke Biotechnology Co.Ltd. The sequences of shRNAs and siRNAs are summarized in Data S2 . Plasmids and siRNAs transfection was performed with jetPRIME reagent (Polyplus) according to the manufacturer’s instructions. Transfection efficiency was determined by reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR). ( 3 ) RT-qPCR: Total RNA was extracted using Trizol (ABclonal), and 1000 ng of RNA was reversed transcribed into cDNA using the RT Master Mix for qPCR II (MCE). And qPCR was performed using the 2X Universal SYBR Green Fast qPCR Mix (ABclonal) and the CFX Opus 384 Real-Time PCR System (BioRad). Gene expression values were normalized to GADPH, which was used as an internal control. The primer sequences used in this study are listed in Data S2 . ( 4 ) Cell viability assay: Cell viability was measured using Cell Counting Kit 8 (CCK-8, Servicebio) according to the manufacturer's instructions. At 0, 24, 48, and 72h, absorbance was measured at 450 nm by using a microplate reader to estimate the cell viability. All experiments were performed in biological triplicate. ( 5 ) EdU staining assay: The EdU assay was performed according to the manufacturer's instructions (BeyoClick™ EdU Cell Proliferation Kit with Alexa Fluor 555; Beyotime, China). Briefly, cells were cultured in 24-well plates, treated with 500µL of medium containing 10 µM EdU at 37.0°C under a 5.0% CO2 atmosphere for 2h, fixed with 4% paraformaldehyde for 20 min, incubated with PBS containing 0.3% Triton X‐100 for 10 min, and treated with Click Reaction Solution for 30 min. The nuclei were counterstained with Hoechst 33342 for 10 min. Three randomly selected regions from each group were visualized using a fluorescence microscope (Leica, Wetzlar, Germany). ( 6 ) Wound‐healing assay: GBM cells from different transfection groups were seeded in 6‐well plates at 1 × 106 cells/well and grown overnight to form monolayers. A sterile tip was used to create a wound line (cell‐free area) across the culture plate surface, and the deplated cells were removed by washing with PBS. Then, cells were cultured in FBS‐free DMEM under a humidified 5.0% CO2 atmosphere at 37.0°C for 24 hours, and images of the wound lines were acquired using a phase‐contrast microscope (Zeiss, Germany). Each assay was performed in triplicate. The migration ability of GBM cell was estimated by measuring the scratch width. ( 7 ) Transwell invasion assay: For the invasion assay, the Transwell upper chambers (Corning) were pre‐coated with Matrigel (Corning) before cell seeding. After serum starvation for 24 hours, GBM cells from different transfection groups were seeded in the upper chambers at 105 cells/200 µL serum‐free DMEM, while 600µL of DMEM with 20% FBS was added to the lower chambers. After 24 hours, all non‐invaded cells on the filter side of the upper chamber were removed by using a cotton swab, and the polycarbonate membrane of the upper chamber was fixed with 4% paraformaldehyde for 30 min, rinsed three times with PBS, and stained with 0.3% crystal violet for 10 min. Invaded cells were counted under a microscope (Zeiss, Germany). All Transwell assays were repeated at least thrice. 4.18 Statistical Analysis All data processing, statistical analyses, and graphing were performed using R software (version 4.2.2). Normality was assessed using the Shapiro-Wilk test. The Wilcoxon rank-sum test was used to compare continuous variables between two groups. The correlation between two continuous variables was evaluated using Spearman or Pearson correlation coefficient. Survival analysis and Kaplan-Meier curve plotting were performed using the survival and survminer packages. The Benjamini-Hochberg method was used to correct the FDR value obtained from multiple comparisons of P values. All statistical tests were two-sided, with P < 0.05 indicating statistical significance. Declarations Acknowledgements Author contributions Conceptualization: YDM, ZYZ, YCJ, ZCL Data curation: ZLW, ZYM, WCD, MKW, WWW, JY, HRL, WYL, YHY, TC, CYM, MMY, JF, SL, BY, YWZ, RKC, DYS, GY Formal analysis: ZYZ, ZLW, ZQL, QCS, YSZ, JXD Project administration: YDM, ZYZ, YCJ, XZL, ZCL, HRZ, DL Writing – original draft: ZYZ, ZLW, RL, DLP, JDL, YNQ Writing – review & editing: ZYZ, ZLW Competing interests Authors declare that they have no competing interests Data Availability Statement The raw WES and RNA-seq data generated in this study have been deposited in the Genome Sequence Archive (GSA) database under accession code HRA006184, accessible at [https://ngdc.cncb.ac.cn/gsa-human/browse/HRA006184]. Additionally, the raw MS-based proteomics data are deposited in the Open Archive for Miscellaneous Data (OMIX) database under accession code OMIX005975, available at [https://ngdc.cncb.ac.cn/omix/release/OMIX005975]. Due to data privacy laws, the raw radiomics and pathomics data are protected and are not available. However, the processed WES, RNA-seq, MS-proteomics, radiomics, pathomics, scRNA-seq and ST data can be accessed via Zenodo at [https://doi.org/10.5281/zenodo.14897960]. For all other data from this study, inquiries can be made to the corresponding authors upon reasonable request. Funding: National Natural Science Foundation of China grant 82273493 (ZYZ) National Natural Science Foundation of China grant U20A2017 (ZCL) National Natural Science Foundation of China grant 82173090 (XZL) National Natural Science Foundation of China grant U1904148 (XZL) Science and Technology Program of Henan Province grant 242102311107 (YCJ) Science and Technology Research and Development Joint Fund of Henan Province grant 242301420014 (ZYZ) References Q. T. Ostrom et al., CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2016-2020. Neuro Oncol 2023;25(12 Suppl 2):iv1-iv99. H. Yan et al., IDH1 and IDH2 mutations in gliomas. N Engl J Med 2009;360(8):765-773. D. N. Louis et al., The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro Oncol 2021;23(8):1231-1251. M. J. van den Bent et al., Primary brain tumours in adults. Lancet 2018;392(10145):432-446. R. Mathur et al., Glioblastoma evolution and heterogeneity from a 3D whole-tumor perspective. Cell 2024;187(2):446-463.e16. L. B. Wang et al., Proteogenomic and metabolomic characterization of human glioblastoma. Cancer Cell 2021;39(4):509-528.e20. K. H. Kim et al., Integrated proteogenomic characterization of glioblastoma evolution. Cancer Cell 2024;42(3):358-377.e8. H. S. Phillips et al., Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell 2006;9(3):157-173. R. G. Verhaak et al., Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell 2010;17(1):98-110. M. Ceccarelli et al., Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma. Cell 2016;164(3):550-563. M. Wu et al., Multi-omics and pharmacological characterization of patient-derived glioma cell lines. Nat Commun 2024;15(1):6740. Q. Wang et al., Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell 2018;33(1):152. K. White et al., Identification, validation and biological characterisation of novel glioblastoma tumour microenvironment subtypes: implications for precision immunotherapy. Ann Oncol 2023;34(3):300-314. F. Guan et al., Biological underpinnings of radiomic magnetic resonance imaging phenotypes for risk stratification in IDH wild-type glioblastoma. J Transl Med 2023;21(1):841. X. Wu et al., Big data and artificial intelligence in cancer research. Trends Cancer 2024;10(2):147-160. N. Jamshidi et al., Illuminating radiogenomic characteristics of glioblastoma multiforme through integration of MR imaging, messenger RNA expression, and DNA copy number variation. Radiology 2014;270(1):1-2. F. Carrillo-Perez et al., Applications of artificial intelligence in dentistry: A comprehensive review. J Esthet Restor Dent 2022;34(1):259-280. J. Lipkova et al., Artificial intelligence for multimodal data integration in oncology. Cancer Cell 2022;40(10):1095-1110. D. Tamborero et al., OncodriveCLUST: exploiting the positional clustering of somatic mutations to identify cancer genes. Bioinformatics 2013;29(18):2238-2244. M. Kanehisa, S. Goto, KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28(1):27-30. A. Liberzon et al., The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1(6):417-425. M. Milacic et al., The Reactome Pathway Knowledgebase 2024. Nucleic Acids Res 2024;52(D1):D672-D678. M. Reich et al., GenePattern 2.0. Nat Genet 2006;38(5):500-501. J. Yan et al., Image-based deep learning identifies glioblastoma risk groups with genomic and transcriptomic heterogeneity: a multi-center study. Eur Radiol 2023;33(2):904-914. P. Karschnia et al., Development and validation of a clinical risk model for postoperative outcome in newly diagnosed glioblastoma: a report of the RANO resect group. Neuro Oncol 2025;27(4):1046-1060. J. Park et al., Transcriptome profiling-based identification of prognostic subtypes and multi-omics signatures of glioblastoma. Sci Rep 2019;9(1):10555. V. M. Ravi et al., Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma. Cancer Cell 2022;40(6):639-655.e13. Galbo PM Jr et al., Molecular Features of Cancer-associated Fibroblast Subtypes and their Implication on Cancer Pathogenesis, Prognosis, and Immunotherapy Resistance. Clin Cancer Res 2021;27(9):2636-2647. C. Neftel et al., An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell 2019;178(4):835-849.e21.. K. Woroniecka et al., T-Cell Exhaustion Signatures Vary with Tumor Type and Are Severe in Glioblastoma. Clin Cancer Res 2018;24(17):4175-4186. F. S. Ng et al., Progenitor-like traits contribute to patient survival and prognosis in oligodendroglial tumors. Clin Cancer Res 2012;18(15):4122-4135. S. Oh et al., Integrated pharmaco-proteogenomics defines two subgroups in isocitrate dehydrogenase wild-type glioblastoma with prognostic and therapeutic opportunities. Nat Commun 2020;11(1):3288. P. A. Yushkevich et al., User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage 2006;31(3):1116-1128. J. J. M. van Griethuysen et al., Computational Radiomics System to Decode the Radiographic Phenotype. Cancer Res 2017;77(21):e104-e107. S. Kalra et al., Yottixel - An Image Search Engine for Large Archives of Histopathology Whole Slide Images. Med Image Anal 2020;65:101757. D. R. Stirling et al., CellProfiler 4: improvements in speed, utility and usability. BMC Bioinformatics 2021;22(1):433. H. Liu et al., Mime: A flexible machine-learning framework to construct and visualize models for clinical characteristics prediction and feature selection. Comput Struct Biotechnol J 2024;23:2798-2810. A. Mayakonda et al., Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28(11):1747-1756. S. Xu et al., Using clusterProfiler to characterize multiomics data. Nat Protoc 2024;19(11):3292-3320. I. Kerseviciute, J. Gordevicius, aPEAR: an R package for autonomous visualization of pathway enrichment networks. Bioinformatics 2023;39(11):btad672. Z. Gu, Complex heatmap visualization. Imeta 2022;1(3):e43. K. Yoshihara et al., Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. P. Charoentong et al., Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18(1):248-262. P. Jiang et al., Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24(10):1550-1558. Mellman I et al., The cancer-immunity cycle: Indication, genotype, and immunotype. Immunity 2023;56(10):2188-2205. Y. Hao et al., Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42(2):293-304. D. M. Cable et al., Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol 2022;40(4):517-526. D. Sun et al., Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol 2022;40(4):527-538. Additional Declarations No competing interests reported. Supplementary Files Supplementarymaterials.docx FigureS1.tif FigureS2.tif FigureS3.tif FigureS4.tif FigureS5.tif FigureS6.tif FigureS7.tif FigureS8.tif FigureS9.tif DataS1.xlsx DataS2.xlsx Cite Share Download PDF Status: Published Journal Publication published 07 Mar, 2026 Read the published version in Molecular Cancer → Version 1 posted Editorial decision: Revision requested 07 Oct, 2025 Reviews received at journal 30 Sep, 2025 Reviews received at journal 21 Sep, 2025 Reviewers agreed at journal 04 Sep, 2025 Reviewers agreed at journal 04 Sep, 2025 Reviewers invited by journal 04 Sep, 2025 Editor assigned by journal 02 Sep, 2025 Submission checks completed at journal 02 Sep, 2025 First submitted to journal 01 Sep, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7510857","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":510820343,"identity":"1b66245b-cd10-4075-b411-3322d068fa06","order_by":0,"name":"Zhenyu Zhang","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Zhenyu","middleName":"","lastName":"Zhang","suffix":""},{"id":510820344,"identity":"e5055151-2109-4426-aff0-c51dd564abf2","order_by":1,"name":"Zilong Wang","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Zilong","middleName":"","lastName":"Wang","suffix":""},{"id":510820345,"identity":"b477fc9a-eb82-42d4-9fc0-971796488d81","order_by":2,"name":"Ran Li","email":"","orcid":"","institution":"Institute of Basic Medical Sciences, Hangzhou City University","correspondingAuthor":false,"prefix":"","firstName":"Ran","middleName":"","lastName":"Li","suffix":""},{"id":510820346,"identity":"1922641a-dce3-4f33-b54f-c25936f40baa","order_by":3,"name":"Dongling Pei","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Dongling","middleName":"","lastName":"Pei","suffix":""},{"id":510820347,"identity":"516ea526-bb89-4ca2-9545-f11167ea353a","order_by":4,"name":"Jingdian Liu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Jingdian","middleName":"","lastName":"Liu","suffix":""},{"id":510820348,"identity":"ff088cb5-9ae6-4b35-a70a-69e863bda943","order_by":5,"name":"Yuning Qiu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Yuning","middleName":"","lastName":"Qiu","suffix":""},{"id":510820349,"identity":"b06a022f-9b28-4384-88fb-4040c147d207","order_by":6,"name":"Zaoqu Liu","email":"","orcid":"","institution":"Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Zaoqu","middleName":"","lastName":"Liu","suffix":""},{"id":510820350,"identity":"e4f0dcb2-f124-482c-b795-e178a19738ca","order_by":7,"name":"Minkai Wang","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Minkai","middleName":"","lastName":"Wang","suffix":""},{"id":510820351,"identity":"bdb833ed-cfba-4351-ba12-8bb33cb797fc","order_by":8,"name":"Zeyu Ma","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Zeyu","middleName":"","lastName":"Ma","suffix":""},{"id":510820352,"identity":"0a6d775d-2996-4abd-bf3e-ec3ff86e9592","order_by":9,"name":"Wenchao Duan","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Wenchao","middleName":"","lastName":"Duan","suffix":""},{"id":510820353,"identity":"72db2e3e-0eb3-409d-90ae-e3acf283a21a","order_by":10,"name":"Weiwei Wang","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Weiwei","middleName":"","lastName":"Wang","suffix":""},{"id":510820354,"identity":"d5e86597-b748-4457-9f0b-37e49cd5ae4e","order_by":11,"name":"Jing Yan","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Jing","middleName":"","lastName":"Yan","suffix":""},{"id":510820355,"identity":"39c05380-4a21-40ec-90b8-9358f3b6c27d","order_by":12,"name":"Yang Guo","email":"","orcid":"","institution":"Henan Provincial People's Hospital","correspondingAuthor":false,"prefix":"","firstName":"Yang","middleName":"","lastName":"Guo","suffix":""},{"id":510820356,"identity":"3d3656d4-0964-4177-bd22-69eb213e7974","order_by":13,"name":"Haoran Liu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Haoran","middleName":"","lastName":"Liu","suffix":""},{"id":510820357,"identity":"442bdc21-3b56-4d3b-8828-d7e330bb559e","order_by":14,"name":"Wenyuan Li","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Wenyuan","middleName":"","lastName":"Li","suffix":""},{"id":510820358,"identity":"f2cea476-27bf-4488-bd7e-8782dc2fbfc7","order_by":15,"name":"Yinhui Yu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Yinhui","middleName":"","lastName":"Yu","suffix":""},{"id":510820359,"identity":"9d55fa21-bbbe-464e-9839-4e28c6ca3aa5","order_by":16,"name":"Te Chen","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Te","middleName":"","lastName":"Chen","suffix":""},{"id":510820360,"identity":"b7d48ee3-b25a-409d-93fb-db824980aaca","order_by":17,"name":"Caoyuan Ma","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Caoyuan","middleName":"","lastName":"Ma","suffix":""},{"id":510820361,"identity":"f385bf27-4078-4740-b360-95ae9a457e8d","order_by":18,"name":"Miaomiao Yu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Miaomiao","middleName":"","lastName":"Yu","suffix":""},{"id":510820362,"identity":"68404117-f52c-43fc-96d3-f8eb599e2684","order_by":19,"name":"Jing Fu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Jing","middleName":"","lastName":"Fu","suffix":""},{"id":510820363,"identity":"3b230cfa-65eb-4887-a314-8f1a0e35a4a6","order_by":20,"name":"Dingyuan Su","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Dingyuan","middleName":"","lastName":"Su","suffix":""},{"id":510820364,"identity":"31d6505a-7018-4620-bebb-811d1b718bd3","order_by":21,"name":"Sen Li","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Sen","middleName":"","lastName":"Li","suffix":""},{"id":510820365,"identity":"e0e289ed-0938-41cb-acfa-1566e0c5bd75","order_by":22,"name":"Haotian Geng","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Haotian","middleName":"","lastName":"Geng","suffix":""},{"id":510820366,"identity":"f625a1bb-07ed-47a4-9444-38a3307ff894","order_by":23,"name":"Bin Yu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Bin","middleName":"","lastName":"Yu","suffix":""},{"id":510820367,"identity":"458fdca6-8083-4c28-88ba-b52eecf62aaf","order_by":24,"name":"Yingwei Zhen","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Yingwei","middleName":"","lastName":"Zhen","suffix":""},{"id":510820368,"identity":"b8904544-4f3a-48b9-8873-0a90ef76447b","order_by":25,"name":"Ruokun Chen","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Ruokun","middleName":"","lastName":"Chen","suffix":""},{"id":510820369,"identity":"dd22a89f-96b1-4c6f-affb-e98568e7640a","order_by":26,"name":"Qiuchang Sun","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Qiuchang","middleName":"","lastName":"Sun","suffix":""},{"id":510820370,"identity":"19546589-2e0f-4422-927b-073812d8670a","order_by":27,"name":"Yuanshen Zhao","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Yuanshen","middleName":"","lastName":"Zhao","suffix":""},{"id":510820371,"identity":"afc05e9b-6f37-4df1-9e53-3b6e105687a2","order_by":28,"name":"Jingxian Duan","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Jingxian","middleName":"","lastName":"Duan","suffix":""},{"id":510820372,"identity":"a437c6f2-d036-4c29-86df-6d82b5e5d2c2","order_by":29,"name":"Hairong Zheng","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Hairong","middleName":"","lastName":"Zheng","suffix":""},{"id":510820373,"identity":"9e4c8e9a-898a-40ee-9979-2698fb774598","order_by":30,"name":"Dong Liang","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Dong","middleName":"","lastName":"Liang","suffix":""},{"id":510820374,"identity":"51b22595-48c7-4170-99de-181c22cb268d","order_by":31,"name":"Xianzhi Liu","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Xianzhi","middleName":"","lastName":"Liu","suffix":""},{"id":510820375,"identity":"78fbf33d-f60c-4343-a549-f7488248087d","order_by":32,"name":"Zhi-Cheng Li","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Zhi-Cheng","middleName":"","lastName":"Li","suffix":""},{"id":510820376,"identity":"7571ad6a-9243-4055-b0bd-57875b07a00a","order_by":33,"name":"Yuchen Ji","email":"","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":false,"prefix":"","firstName":"Yuchen","middleName":"","lastName":"Ji","suffix":""},{"id":510820377,"identity":"4110a63f-cc68-4b67-a1a4-7579ad653fc3","order_by":34,"name":"Dongming Yan","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAtElEQVRIiWNgGAWjYJACCYYKGIN4LWdI1sLYRooWg+NnD966Oa/O3uAA88HbPAx2eYS1nMlLts7dxsZscIAt2ZqHIbmYsJYDOWbSudt42AwO8JhJ8zAcSGwgqOX8G6CWORI8Bgf4vxGp5QbIlgYDCaAtbMRpkbzxxtg651iCgeRhNmPLOQbJhLXwnc8xvJ1TU2fPd7z54Y03FXaEtSgcgLGYwe4kpB4I5AkaOgpGwSgYBaMAACgfNlFahxlGAAAAAElFTkSuQmCC","orcid":"","institution":"The First Affiliated Hospital of Zhengzhou University","correspondingAuthor":true,"prefix":"","firstName":"Dongming","middleName":"","lastName":"Yan","suffix":""}],"badges":[],"createdAt":"2025-09-01 17:53:19","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7510857/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7510857/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1186/s12943-026-02637-2","type":"published","date":"2026-03-07T15:59:58+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":91190290,"identity":"a2b9b697-0218-41b2-bbff-c5102e6e2c52","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1960071,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDevelopment and validation of the multi-omics prognostic model.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eThe detailed workflow of this study. The best combination of algorithms was screened at each of the five single-omics levels to build the best single-omics prognostic model and calculate the patient risk score, and then the five-omics risk scores were fused together based on a late fusion strategy to obtain a multi-omics risk score.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eMulti-omics data queue presentation. Genomics, bulk transcriptomics, proteomics, pathomics, and radiomics data are used to establish prognostic models. Single-cell transcriptomics and spatial transcriptomics data are used for the analytical validation of risk stratification.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eKaplan-Meier survival curves derived from the multi-omics risk score，demonstrating significant survival differences among risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eDecision curves show a significant advantage of the clinical utility and predictive accuracy of multi-omics models over single-omics models.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E) \u003c/strong\u003eTime-dependent Concordance Index Curve of various omics models demonstrated that the multi-omics model exhibited consistently superior predictive performance across different time points which highlights the enhanced prognostic accuracy and robustness of the multi-omics model compared to individual omics models.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(F) \u003c/strong\u003eThe multi-omics model consistently outperforms individual omics models across all metrics (Cindex, AIC, BIC, IAE, and ISE), highlighting its superior predictive accuracy, reduced error, and balanced complexity.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(G) \u003c/strong\u003eThe curve of Net Reclassification Improvement value from the Mult model versus other single-omics models at various points reflecting the Mult model improvement over the other single-omics models in terms of risk categorization.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(H) \u003c/strong\u003eThe curve of Integrated Discrimination Improvemen value obtained from the Mult model versus other single-omics models at various time points reflecting the Mult model improvement over the other single-omics models in terms of distinguishing individuals with different outcomes.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(I) \u003c/strong\u003eThe bar chart illustrates the changes in the C-index of the Mult models missing any one of the omics datasets, as well as the significance changes in Kaplan-Meier survival analysis, represented by -log10 (p-value) from the log-rank test.\u003c/p\u003e","description":"","filename":"Figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/c78dbb469889c40c6b0119ea.png"},{"id":91189296,"identity":"b3e04c6b-115d-430f-bf49-e55adf4d97ce","added_by":"auto","created_at":"2025-09-12 14:27:16","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":638401,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAnalysis of gene mutation levels revealing the heterogeneity of GBM across different risk groups.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A-B) \u003c/strong\u003eWaterfall plot illustrates the distinct mutation landscapes of the high-risk and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eBar chart shows the genes with significantly different mutation frequencies (\u0026gt; 5%) between high- and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eCancer driver genes identified in patients from different risk groups using the oncodriveCLUST algorithm.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E) \u003c/strong\u003eNetwork diagram of co-mutated genes specific to high- and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(F) \u003c/strong\u003eDot plot depicts the impact of mutated genes on pathways from the KEGG, Hallmark, and Reactome databases, as well as on patients in different risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(G-J) \u003c/strong\u003eThe stacked bar chart illustrates the proportional distribution of pathways obtained from ORA analysis of differentially mutated genes, cancer driver genes, and co-mutated genes between the high- and low-risk groups, as well as pathways with a high mutation burden.\u003c/p\u003e","description":"","filename":"Figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/ccb70ecb55bdfa8e3a0a76a7.png"},{"id":91190287,"identity":"e20a4215-de19-4ef9-a684-a2affe022230","added_by":"auto","created_at":"2025-09-12 14:35:16","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":689524,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAnalysis of CNV levels revealing the heterogeneity of GBM across different risk groups.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eThe ridge plot depicts the landscape of significant copy number variations (CNVs) across the entire genome, comparing high-risk and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eThe volcano plot visualizes the differential analysis of CNV scores between high- and low-risk groups across various chromosomal regions.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eThe cytoband waterfall plot highlights significant differences in CNV amplifications and deletions between high- and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eA detailed G-score representation pinpoints cytobands with significant CNV differences between high- and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E) \u003c/strong\u003eThe lollipop plot illustrates genes associated with cytoband regions exhibiting significant CNV differences, along with their CNV frequencies across high- and low-risk groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(F) \u003c/strong\u003eThe bar chart presents biological pathways impacted by significant copy number variations between high- and low-risk groups.\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/44d20c38efcb1803fc46500c.png"},{"id":91190297,"identity":"00cb8bd7-a1e4-443c-929d-5c26bb663522","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1513069,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional enrichment analysis of integrated transcriptomic and proteomic data.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eThe network diagram illustrates the clustering of significantly activated pathways in the transcriptome and proteome of high-risk patients.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eThe network diagram illustrates the clustering of significantly activated pathways in the transcriptome and proteome of low-risk patients.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C-D) \u003c/strong\u003eThe proportion of all pathway types in the high-risk group identified through Gene Set Enrichment Analysis (GSEA) conducted separately on transcriptomic and proteomic data.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E-F) \u003c/strong\u003eThe proportion of all pathway types in the low-risk group identified through Gene Set Enrichment Analysis (GSEA) conducted separately on transcriptomic and proteomic data.\u003c/p\u003e","description":"","filename":"Figure4.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/dfb147e55e6b3c5eb118ea75.png"},{"id":91189304,"identity":"167aea80-ab8e-48b2-a835-6ecbdbb371bf","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":883294,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCross-omics heterogeneity analysis and comparison with transcriptional subtypes.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eThe heatmap illustrates the activation status of characteristic biological pathways associated with each risk group across different omics levels.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B-C) \u003c/strong\u003eThe Sankey diagram illustrates the relationship between the risk groups defined in this study and various transcriptomic subtypes reported in previous research.\u003c/p\u003e","description":"","filename":"Figure5.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/a27055c0c91fcb7e9c5ab7cc.png"},{"id":91189307,"identity":"ab201e1b-031d-411c-8989-b1a1a0d1fd87","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":6287119,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eExperimental verification of biological functions of the targets.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eSchematic representation of the target screening process.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B-C) \u003c/strong\u003eRepresentative experimental images and statistical results show the effect of EIF3I knockdown on proliferation, migration, and invasion ability of A172 cell.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D-E) \u003c/strong\u003eRepresentative experimental images and statistical results show the effect of EIF3I knockdown on proliferation, migration, and invasion ability of U251 cell.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(F-G) \u003c/strong\u003eRepresentative experimental images and statistical results show the effect of PDIA4 knockdown on proliferation, migration, and invasion ability of A172 cell.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(H-I) \u003c/strong\u003eRepresentative experimental images and statistical results show the effect of PDIA4 knockdown on proliferation, migration, and invasion ability of U118 cell.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(J-K) \u003c/strong\u003eRepresentative experimental images and statistical results show effect of RFT1 knockdown on proliferation, migration, and invasion ability of A172 cell.\u003c/p\u003e\n\u003cp\u003eData are presented as mean ±SD from at least three independent experiments. *P \u0026lt; 0.05, **P \u0026lt; 0.01, ***P \u0026lt; 0.001, and ****P \u0026lt; 0.0001.\u003c/p\u003e","description":"","filename":"Figure6.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/61483207dc223d388b822848.png"},{"id":91190295,"identity":"3f7a0822-8a15-4d88-91bb-7bde8511228f","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":777224,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSingle-cell transcriptome analysis based on risk grouping based on multi-omics fusion model.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eUniformManifold Approximation and Projection (UMAP) plot displaying the integrated cellmap, which consists of 12 annotated cell types. Cells are colored by clusters.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eThe cumulative bar chart illustrates the proportion of various cell subtypes across different groups, highlighting compositional differences between the two groups.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eThe volcano plot visualizes the differential analysis results of cell subtype proportions between groups, identifying significant variations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eThe multi-group volcano plots showcase differential gene expression analysis within the same cell subtype across different groups, revealing transcriptional changes.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E) \u003c/strong\u003eThe bar graph shows the functional differences of different cell subtypes between the two risk groups.\u003c/p\u003e","description":"","filename":"Figure7.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/1fe3844930d9b5c4d643f2bf.png"},{"id":91189299,"identity":"4dfb2fe2-d30d-4017-aa3c-ff05d19b0a9e","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":2385804,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePrognosis and target analysis at the single-cell transcriptome level.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eBar plots and a pie chart showing the distribution of Scissor+ and Scissor− cells identified from the integrated single-cell transcriptomic data. Scissor+ cells (65.74%) were associated with poor prognosis, while Scissor− cells (34.26%) were associated with better survival outcomes. The left and right bar plots show the dominant cell types in Scissor+ and Scissor− populations, respectively.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eHeatmap of differentially expressed genes across major cell subpopulations, annotated with functional enrichment results from Gene Ontology (GO_BP), KEGG, and WikiPathways databases. Enrichment terms highlight pathways associated with immune response, cell cycle regulation, and hypoxia.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eUMAP plots visualizing the expression levels of PDIA4, EIF3I, and RPT1 across the entire single-cell population. High expression levels of PDIA4 and EIF3I are enriched in specific niches.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eViolin plots show the expression distribution of PDIA4, EIF3I, and RPT1 across annotated cell types, indicating their cell type-specific expression patterns.\u003c/p\u003e","description":"","filename":"Figure8.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/da03b01fc11b35e0e02ff71f.png"},{"id":91189317,"identity":"875e4a33-d097-4834-a918-1b3000f39747","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":6985379,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSpatial transcriptome-wide distribution of cell types and cell niches.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eSpatial maps showing the distribution of major cell types across 10 GBM tissue sections based on deconvolution of spatial transcriptomic data. Each row represents a distinct cell type, and each column corresponds to a different GBM sample. Heatmap color intensity reflects the predicted proportion of each cell type within individual spatial spots.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eHeatmap showing the distribution of cell subsets in each cell niche.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eSpatial map showing how spatial points are clustered into different cellular microenvironments in an unsupervised clustering process.\u003c/p\u003e","description":"","filename":"Figure9.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/8ba01e8c4684a6c53c8434e4.png"},{"id":91189334,"identity":"cef69f8f-6efc-42a1-aca1-d02cfb912d08","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"png","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":4220181,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDistribution of targets in the cell niche and on idle slices.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eUMAP plot defining seven spatial niches based on similarities in cell-type composition derived from spatial transcriptomic data.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(B) \u003c/strong\u003eHigh- and low-risk regions were mapped onto the spatial niche UMAP, revealing spatial heterogeneity of prognostic risk across niches.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(C) \u003c/strong\u003eTen patients were projected onto the spatial niche UMAP, illustrating inter-individual variability.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(D) \u003c/strong\u003eBar plot summarizing the relative proportions of each niche type across the 10 patients, indicating heterogeneity in niche occupancy.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(E) \u003c/strong\u003eSpatial distributions of major cell types within the defined niches, highlighting niche-specific enrichment patterns.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(F) \u003c/strong\u003eSpatial expression profiles of the candidate therapeutic targets PDIA4, EIF3I, and RFT1 across all samples.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(G) \u003c/strong\u003eNiche-level spatial expression patterns of PDIA4, EIF3I, and RFT1.\u003c/p\u003e","description":"","filename":"Figure10.png","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/f137b81fd0fff34ecb4bd521.png"},{"id":104251900,"identity":"aeb2633f-6e08-453e-8cd0-fb57fe0114aa","added_by":"auto","created_at":"2026-03-09 16:15:58","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":23546416,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/650b434f-ef63-4248-9f7e-fdb8d085e51b.pdf"},{"id":91190303,"identity":"50d9d4a7-358a-4227-b1cc-9776a7bf7e01","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":10626530,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarymaterials.docx","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/ea7e4bee07fd919b492f0e46.docx"},{"id":91189308,"identity":"3a9a21c2-6d72-4f5d-9314-15294ba001f8","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"tif","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":4983290,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS1.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/cfaf5aa90677fd1a5070b1b8.tif"},{"id":91192795,"identity":"ac406cf8-5ec2-4c38-a3a7-ac315c8254de","added_by":"auto","created_at":"2025-09-12 14:43:17","extension":"tif","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":6005602,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS2.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/4b5591a068155162108735ef.tif"},{"id":91192796,"identity":"8dbcaca0-77aa-42ba-8897-8c451130d081","added_by":"auto","created_at":"2025-09-12 14:43:17","extension":"tif","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":4828144,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS3.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/83dc2e6251b85dbc20a1c4e3.tif"},{"id":91192801,"identity":"eacbc55d-a397-441b-a1fa-f23726f4da08","added_by":"auto","created_at":"2025-09-12 14:43:17","extension":"tif","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":4808716,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS4.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/cbabb38988a6805cc878df92.tif"},{"id":91190307,"identity":"f3f8f966-cc36-4527-a3b4-580b48df6648","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"tif","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":5677588,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS5.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/5164fbfcb9bd4d3a80834e89.tif"},{"id":91194273,"identity":"1b50a833-4ec2-4aaa-b2ea-e7acb9663d30","added_by":"auto","created_at":"2025-09-12 14:51:17","extension":"tif","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":5192504,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS6.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/a86628f02078f9e9f8c008be.tif"},{"id":91194274,"identity":"a3e8d3fd-737d-4382-a2bf-ca2fbad492d2","added_by":"auto","created_at":"2025-09-12 14:51:18","extension":"tif","order_by":7,"title":"","display":"","copyAsset":false,"role":"supplement","size":13169380,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS7.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/bf88688bc62fb598fef72a4a.tif"},{"id":91194271,"identity":"6ac8f5cf-1b81-4a26-b6b8-968bb62e488d","added_by":"auto","created_at":"2025-09-12 14:51:17","extension":"tif","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":4066556,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS8.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/d58dcf8627f74fe2b024157d.tif"},{"id":91189315,"identity":"b9d051a2-2764-45d5-8887-ce7308e5d05d","added_by":"auto","created_at":"2025-09-12 14:27:17","extension":"tif","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":4368956,"visible":true,"origin":"","legend":"","description":"","filename":"FigureS9.tif","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/201a468cae8957de738f9bdf.tif"},{"id":91190309,"identity":"f785f942-d5c9-48ef-b74b-43f8cd880ef6","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"xlsx","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":72010,"visible":true,"origin":"","legend":"","description":"","filename":"DataS1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/c3f270bdf8f480cf1323cef7.xlsx"},{"id":91190308,"identity":"688c38cd-c50b-44c1-8621-407bb77fcc9f","added_by":"auto","created_at":"2025-09-12 14:35:17","extension":"xlsx","order_by":11,"title":"","display":"","copyAsset":false,"role":"supplement","size":92462,"visible":true,"origin":"","legend":"","description":"","filename":"DataS2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7510857/v1/98223fb143dd444cc78dc9e5.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Machine Learning Model on Multi-Omics Data Enables Risk Stratification and Identifies Molecular Heterogeneity and Therapeutic Targets in Glioblastoma","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eAdult-type diffuse glioma is the most common primary malignant tumor of the central nervous system(\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e), with IDH mutation serving as a key molecular marker that stratifies gliomas into two prognostically and genetically distinct groups(\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e). The latest WHO classification recognizes IDH-mutant gliomas and IDH wild-type GBM as separate tumor types(\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e). Standard GBM treatment includes surgical resection followed by chemoradiotherapy and temozolomide(\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e), yet the disease's profound intertumoral heterogeneity-spanning the genome, transcriptome, and proteome\u0026mdash;poses significant therapeutic challenges(\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e). Thus, a comprehensive investigation of IDH wild-type GBM heterogeneity is crucial. High-throughput sequencing and multi-omics analyses have significantly advanced our understanding of the molecular characteristics of GBM(\u003cspan additionalcitationids=\"CR7\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e). Verhaak et al. identified four transcriptomic subtypes, each linked to distinct genetic alterations(\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e), while Ceccarelli et al. used DNA methylation and gene expression profiling to further refine glioma subtypes(\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e). Other studies have classified GBM based on pharmacological response(\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e), tumor microenvironment (TME) composition(\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e), and radiogenomic features, revealing distinct biological and immune landscapes(\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e). However, most existing classification models rely on single-omics data, limiting their ability to fully capture GBM complexity. With advancements in sequencing and intelligent imaging analysis, multimodal data integration for precision oncology has emerged as a frontier in cancer research(\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e). Each modality offers unique insights into tumor biology, and their integration provides a more comprehensive understanding of tumor heterogeneity(\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e). In this context, artificial intelligence has transformed cancer diagnostics, classification, molecular profiling, and treatment selection(\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eThis study leverages a late-fusion strategy combined with multiple machine learning algorithms to integrate multimodal data\u0026mdash;including radiomics, pathology, genomics, transcriptomics, and proteomics\u0026mdash;for the risk stratification of glioblastoma (GBM) patients. This stratification framework enables a comprehensive characterization of GBM heterogeneity across bulk multi-omics, single-cell transcriptomics, and spatial transcriptomics, while also identifying potential therapeutic targets for phenotypic validation.\u003c/p\u003e"},{"header":"2. Results","content":"\u003cp\u003e\u003cstrong\u003e2.1 Multimodal Data Summary of GBM in this study\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll training datasets used in this study were obtained from the First Affiliated Hospital of Zhengzhou University. A total of 1,187 adult patients diagnosed with IDH wild-type GBM were recruited. Among them, 185 patients underwent RNA sequencing (RNA-seq), 166 patients were analyzed by mass spectrometry (MS) for proteomics research, and 167 patients were subjected to whole-exon sequencing (WES). Meanwhile, 936 patients had preoperative multiparametric magnetic resonance imaging (MRI) scans, including T1-weighted imaging (T1WI), contrast-enhanced T1-weighted imaging (CE-T1WI), T2-weighted imaging (T2WI), fluid - attenuated inversion recovery (FLAIR), and apparent diffusion coefficient (ADC). Image segmentation and feature extraction were performed on the MRI data of these 936 patients. Histological whole-slide images (WSIs) were obtained from 1122 patients (\u003cstrong\u003eFigure 1A\u003c/strong\u003e). Additionally, single-cell transcriptome (scRNA-seq) data were available for 26 patients, and spatial transcriptome (ST) data were accessible for 10 patients. For a comprehensive overview of all FAHZZU multimodal data, please refer to \u003cstrong\u003eFigure 1B\u003c/strong\u003e and \u003cstrong\u003eData S1\u003c/strong\u003e.\u0026nbsp;The external validation datasets in this study were obtained from multiple institutions. WES data from 512 patients were sourced from The Cancer Genome Atlas\u0026zwnj; (TCGA). RNA-seq data included 130 cases from TCGA, 72 from Chinese Glioma Genome Atlas (CGGA) -325, and 107 from CGGA-693. MRI data were collected from 108 patients at Henan Provincial People\u0026apos;s Hospital (HPPH) and 28 from TCGA. WSI data included 125 cases from HPPH and 144 from TCGA. MS data were obtained from 31 patients at Samsung Medical Center (SMC) and 91 from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.2 Development and Validation of the Multimodal Prognostic Model\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDuring model development, we implemented a late fusion strategy (\u003cstrong\u003eFigure 1A\u003c/strong\u003e)(\u003cem\u003e18\u003c/em\u003e). Initially, the optimal combination of machine learning algorithms was identified at each unimodal level to construct prognostic models for patient risk scoring (\u003cstrong\u003eFigure S1-S5\u003c/strong\u003e). These five unimodal models were then integrated using a Cox proportional hazards regression model with 10-fold cross-validation, resulting in a comprehensive multimodal fusion-based risk stratification model. Comparative analysis demonstrated that the multimodal model significantly outperformed single-omics models across multiple evaluation metrics. It achieved superior risk stratification, with a highly significant survival difference between high- and low-risk groups (\u003cstrong\u003eFigure 1C\u003c/strong\u003e). Decision Curve Analysis (DCA) confirmed greater net clinical benefits (\u003cstrong\u003eFigure 1D\u003c/strong\u003e), while the Time-dependent Concordance Index (Time-Cindex) and overall Concordance Index (C-index) consistently exceeded those of single-omics models, indicating enhanced predictive accuracy (\u003cstrong\u003eFigure 1E\u003c/strong\u003e). The model also exhibited lower Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values, reflecting better fitness and parsimony, along with reduced Integrated Absolute Error (IAE) and Integrated Squared Error (ISE), highlighting improved calibration and reliability (\u003cstrong\u003eFigure 1F\u003c/strong\u003e). Further, the multimodal model showed significant improvements in Integrated Discrimination Improvement (IDI) and Net Reclassification Improvement (NRI), confirming its superior stratification capability (\u003cstrong\u003eFigure 1G-H\u003c/strong\u003e). A comparative evaluation of four-omics versus five-omics integration models revealed the five-omics model as the superior approach, demonstrating substantial advantages across multiple dimensions (\u003cstrong\u003eFigure 1I\u003c/strong\u003e).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.3 Genomic Alteration Profiles of High- and Low-Risk Group\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo assess the impact of gene mutations on different risk groups, we analyzed key mutational features, including mutation types, genetic variation patterns, single nucleotide substitution profiles, and tumor mutational burden (TMB) (\u003cstrong\u003eFigure S6A\u003c/strong\u003e). Splice site mutations (P = 0.008) and C\u0026gt;A substitutions (P = 0.041) were significantly enriched in the high-risk group, though their overall proportions were low (\u003cstrong\u003eFigure S6B\u003c/strong\u003e). TMB showed no significant differences between risk groups and had a weak correlation with risk scores, suggesting a limited role in tumor progression (\u003cstrong\u003eFigure S6C-D\u003c/strong\u003e). Comparing the top 20 most frequently mutated genes in both groups revealed six shared mutations (TNN, TP53, EGFR, MUC16, NF1, PTEN), with TP53 showing a 9% higher mutation frequency in the high-risk group (\u003cstrong\u003eFigure 2A-B\u003c/strong\u003e). Functional analysis indicated that high-risk group mutations were enriched in cell cycle regulation, gene expression, and immune activity, while low-risk group mutations were associated with neuronal function, tumor stroma interactions, and substance transport (\u003cstrong\u003eFigure 2G\u003c/strong\u003e). OncodriveCLUST analysis identified risk-group-specific driver genes: high-risk drivers were linked to cell cycle regulation and immune activity, whereas low-risk drivers were associated with DNA repair and metabolic pathways (\u003cstrong\u003eFigure 2D, 2H\u003c/strong\u003e)(\u003cem\u003e19\u003c/em\u003e). Co-mutation analysis further highlighted distinct biological processes, with high-risk group co-mutations enriched in gene expression regulation and disease progression, and low-risk group co-mutations associated with metabolism, neuronal function, and tumor stroma interactions (\u003cstrong\u003eFigure 2E, 2I\u003c/strong\u003e). Pathway-level mutation analysis using KEGG(\u003cem\u003e20\u003c/em\u003e), Hallmark(\u003cem\u003e21\u003c/em\u003e), and Reactome(\u003cem\u003e22\u003c/em\u003e) databases revealed distinct functional alterations: high-risk group mutations affected cellular stimulation responses, material metabolism, and gene expression, whereas low-risk group mutations impacted immune activity, coagulation, and signal transduction (\u003cstrong\u003eFigure 2F, 2J\u003c/strong\u003e). Given the impact of copy number variations (CNVs) on gene expression, we examined CNV frequencies (\u003cstrong\u003eFigure 3A\u003c/strong\u003e). GISTIC 2.0(\u003cem\u003e23\u003c/em\u003e) identified 16 chromosomal cytobands with significant CNV differences, categorized into amplifications and deletions specific to each risk group (\u003cstrong\u003eFigure 3B-D\u003c/strong\u003e). Pathway enrichment analysis of CNV-altered genes revealed that high-risk amplifications were associated with tumor signaling, cell cycle regulation, and metabolic dysregulation, while deletions affected DNA repair and apoptosis. Low-risk deletions were linked to abnormal glycosylation and oncogenic pathways (EPH-Ephrin, Rho GTPase, TGF-\u0026beta;), whereas low-risk amplifications showed no significant pathway enrichment (\u003cstrong\u003eFigure 3E-F\u003c/strong\u003e).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTranscriptome-Proteome Heterogeneity in GBM\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo investigate GBM molecular heterogeneity at both transcriptomic and proteomic levels, we performed differential expression analyses to identify key altered genes and proteins between risk groups. Pathway enrichment clustering revealed distinct core biological processes: in the high-risk group, pathways related to GPCR signaling, immune responses, extracellular matrix remodeling, RNA splicing, protein synthesis, intracellular transport, and metabolism were enriched (\u003cstrong\u003eFigure 4A\u003c/strong\u003e). The low-risk group showed enrichment in pathways associated with GPCR signaling, neuronal processes, calcium ion signaling, ion transport, and neuron-tumor communication (\u003cstrong\u003eFigure 4B\u003c/strong\u003e). Gene Set Enrichment Analysis (GSEA) further validated these findings, highlighting consistent transcriptomic and proteomic enrichment patterns. In the high-risk group, dominant pathways included cell cycle regulation, gene expression, DNA repair, metabolism, and immune responses (\u003cstrong\u003eFigure 4C-D\u003c/strong\u003e), while the low-risk group enriched in synaptic activity, signal transduction, lipid metabolism, and substance transport (\u003cstrong\u003eFigure 4E-F\u003c/strong\u003e). Gene Set Variation Analysis (GSVA) confirmed the robustness of these results, further supporting pathway activity consistency (\u003cstrong\u003eFigure S7B\u003c/strong\u003e).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.4 Immune Microenvironment Differences Between Risk Groups\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo explore immune microenvironment disparities between the high- and low-risk groups, we conducted GSEA on immune-related pathways, revealing heightened immune response in the high-risk group (\u003cstrong\u003eFigure S8A\u003c/strong\u003e). However, ESTIMATE analysis showed no significant differences in stromal score, immune score, or tumor purity between the groups (\u003cstrong\u003eFigure S8B\u003c/strong\u003e). Immune cell Proportion Score (IPS) analysis indicated higher effector and suppressor cell scores in the low-risk group, suggesting increased activity of these immune cells (\u003cstrong\u003eFigure S8C\u003c/strong\u003e). Tumor Immune Dysfunction and Exclusion (TIDE) and Cancer Immunity Cycle (CIC) analyses revealed contrasting immune features: the high-risk group exhibited significant Cancer-associated fibroblasts (CAF) infiltration, suggesting a strong immunosuppressive response, while the low-risk group showed greater T-cell dysfunction (\u003cstrong\u003eFigure S8D-E\u003c/strong\u003e). CIC analysis further highlighted immune differences, with higher macrophage and neutrophil recruitment in the high-risk group at Step 4 and more pronounced T-cell infiltration in the low-risk group at Step 5 (\u003cstrong\u003eFigure S8F-G\u003c/strong\u003e). Tumor Immune Prediction (TIP) analysis revealed no significant differences in overall immune activity scores between groups (\u003cstrong\u003eFigure S8H\u003c/strong\u003e). In summary, the high-risk group demonstrated elevated immune pathway activity but strong immunosuppressive CAF features, while the low-risk group exhibited higher effector cell activity but greater T-cell dysfunction. These distinct immune profiles may influence patient prognosis and therapeutic outcomes.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.5 Integrated Biological Pathway Analysis Across Multi-Omics Levels\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo investigate the heterogeneity of GBM across different multi-omics levels, we performed differential analysis, Over-Representation Analysis (ORA), and GSVA for each omics dataset. The most significantly altered biological pathways at different omics levels were visualized (\u003cstrong\u003eFigure 5A\u003c/strong\u003e). The results revealed consistency in the heterogeneity patterns across multi-omics levels. In the high-risk group, pathways related to cell proliferation, cell cycle regulation, immune responses, and oncogenic signaling were highly active, whereas in the low-risk group, pathways associated with neuronal activity and molecular transport exhibited greater enrichment. Furthermore, we compared our risk classification with previously published GBM subtype classifications (\u003cstrong\u003eFigure 5B-C\u003c/strong\u003e)(\u003cem\u003e8, 9, 11, 12\u003c/em\u003e). The results demonstrated that our risk stratification differed considerably from prior transcriptomic subtypes of GBM. Notably, each of our risk groups encompassed multiple transcriptional subtypes from different studies, suggesting that the multimodal risk stratification captures a broader spectrum of GBM heterogeneity beyond traditional transcriptomic subtyping.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.6 Potential Therapeutic Targets Identification and Experimental Validation\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBased on risk stratification using the multimodal fusion model, we screened highly expressed genes in the high-risk group across multiple omics levels and multi-omics databases as potential therapeutic targets. \u003cstrong\u003eFigure 6A\u003c/strong\u003e illustrates the target screening workflow. The results of each screening step are shown in \u003cstrong\u003eFigure S9\u003c/strong\u003e, ultimately identifying 18 target candidates, which were ranked based on their prognostic significance across multiple omics layers and databases. Following a review of previously published literature and preliminary experiments, PDIA4, RFT1, and EIF3I were selected as experimental targets.\u003c/p\u003e\n\u003cp\u003eFirst, we quantified the expression levels of EIF3I, PDIA4, and RFT1 in different GBM cell lines by RT-qPCR. Cell lines with higher baseline expression levels (A172/U251 for EIF3I, A172/U118 for PDIA4, and A172 for RFT1) were selected for functional validation. Biological effects of the three targets were evaluated using EdU and CCK-8 proliferation assays, wound-healing assays, and Transwell invasion assays (\u003cstrong\u003eFigure 6B/D/F/H/J\u003c/strong\u003e). The corresponding statistical results are presented in \u003cstrong\u003eFigure 6C/E/G/I/K\u003c/strong\u003e. The proliferation assay results demonstrated that the silencing of any of PDIA4, RFT1, or EIF3I significantly inhibited the proliferation of GBM cells. Subsequent wound healing assays further confirmed that the knockdown of PDIA4, RFT1, or EIF3I markedly suppressed the migratory abilities of GBM cells. Additionally, the Transwell assays revealed that the PDIA4, RFT1, or EIF3I deficiency inhibited the invasive capacity of GBM cells. Collectively, these findings suggest that PDIA4, RFT1, and EIF3I play pivotal roles in regulation of GBM cell proliferation, migration, and invasion.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.7 Cell-Type\u0026ndash;Specific Risk Stratification and Therapeutic Target Mapping in Single-Cell Resolution\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe performed quality control (QC) on the single cell sequencing data and clustered and annotated 12 different cell types (\u003cstrong\u003eFigure 7A\u0026nbsp;\u003c/strong\u003eand\u003cstrong\u003e\u0026nbsp;Figure 8B\u003c/strong\u003e). We then generated a quasi-bulk mRNA profile for each patient and used a multi-algorithm classifier to integrate multimodal risk stratification to divide patients into high-risk and low-risk groups, allowing comparison of the composition of cell subpopulations (\u003cstrong\u003eFigure 7B\u003c/strong\u003e). The results of subpopulation ratio difference analysis showed that oligodendrocyte subpopulation was significantly enriched in low-risk patients. No significant enrichment of cell type subpopulations was found in the high-risk group (\u003cstrong\u003eFigure 7C\u003c/strong\u003e).\u003c/p\u003e\n\u003cp\u003eGene differential expression and pathway enrichment analyses revealed profound biological distinctions between risk groups (\u003cstrong\u003eFigure 7D-E\u003c/strong\u003e). Across a spectrum of cell type\u0026ndash;specific subpopulations in GBM, consistent functional divergence emerged between high- and low-risk groups, underscoring distinct transcriptional programs driving prognostic heterogeneity. In the AC-like, MES-like, and OPC-like subgroups, high-risk cells were characterized by elevated oxidative phosphorylation and enhanced mitochondrial activity, indicative of metabolic reprogramming that supports cellular proliferation. In contrast, low-risk counterparts were enriched for pathways associated with signaling regulation, inflammatory equilibrium, and homeostatic maintenance. Similarly, G2M-like and NPC-like cells in high-risk patients exhibited increased protein synthesis and nutrient stress responses, whereas low-risk groups favored DNA repair, replication fidelity, and genomic stability\u0026mdash;features associated with restrained proliferation and improved prognosis. In Oligodendrocyte-like cells, high-risk groups showed upregulated translation and ribosomal function, while low-risk cells preferentially activated heat shock, MAPK, and TNF signaling pathways, reflecting a shift from synthetic to stress-adaptive states. Endothelial and Mural subgroups in the high-risk cohort displayed enhanced ECM remodeling, epithelial\u0026ndash;mesenchymal transition (EMT), and interferon-driven inflammatory signaling, whereas low-risk cells maintained vascular integrity, solute transport, and transcriptional balance. Immune lineages mirrored these patterns: high-risk T cells exhibited oxidative metabolism and exhaustion-like phenotypes, while low-risk T cells-maintained cytokine signaling and immune competence. B cells in high-risk tumors displayed interferon-stimulated transcription and translational stress, contrasting with antigen sensing and chromatin-modifying signatures observed in the low-risk group. In the TAM_MG and Mono subpopulations, high-risk states were defined by mitochondrial respiration and immunosuppressive features, while low-risk cells were enriched for innate immune activation via TLR, NOD, NF-\u0026kappa;B, and IL-10 signaling pathways, reflecting preserved immunoregulatory capacity. Together, these findings highlight cell-type\u0026ndash;specific metabolic and immunological reprogramming as fundamental to GBM risk stratification and point toward risk-specific therapeutic vulnerabilities, particularly in metabolism- and immunity-targeted strategies.\u003c/p\u003e\n\u003cp\u003eTo investigate the relationship between single-cell transcriptional states and patient prognosis, we performed Scissor analysis by integrating single-cell RNA-seq profiles with bulk RNA-seq\u0026ndash;derived survival information. This approach identified two prognostically relevant cellular compartments: Scissor+ cells, which were positively correlated with poor survival, and Scissor\u0026ndash; cells, associated with favorable outcomes (\u003cstrong\u003eFigure 8A\u003c/strong\u003e). Notably, scissor+ cells were predominantly composed of malignant populations, including AC-like, MES-like, and G2M-like subtypes, suggesting that transcriptional programs within these tumor cell states contribute directly to adverse prognosis. In contrast, scissor\u0026ndash; cells were primarily of immune origin, particularly TAM_MG, T cells, and Monocytes, indicating that these immune populations may be associated with better clinical outcomes. These results highlight the prognostic relevance of intra-tumoral cellular heterogeneity and emphasize that tumor-intrinsic transcriptional programs, rather than immune compartment dysfunction, are the primary drivers of poor prognosis in this cohort.\u003c/p\u003e\n\u003cp\u003eWe further investigated the single-cell expression patterns of three candidate therapeutic targets\u0026mdash;PDIA4, RFT1, and EIF3I\u0026mdash;previously identified through bulk multi-omics analysis (\u003cstrong\u003eFigure 8C\u003c/strong\u003e). Mapping their expression across single-cell subpopulations revealed that all three targets were predominantly expressed in tumor-associated cell states, rather than in immune or stromal compartments (\u003cstrong\u003eFigure 8D\u003c/strong\u003e). Specifically, PDIA4 and RFT1 were highly enriched in the G2M-like malignant subpopulation, consistent with their association with proliferative programs. In contrast, EIF3I expression was primarily restricted to the MES-like subpopulation, which has been linked to aggressive and therapy-resistant phenotypes. These results suggest that the identified targets are tightly linked to specific tumor cell states and may offer subtype-selective therapeutic opportunities.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2.8 Dissecting Spatial Heterogeneity and Target Localization in GBM via Spatial Transcriptomics\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo further investigate the spatial organization of intratumoral heterogeneity, we performed spatial transcriptomic deconvolution using the RCTD algorithm, guided by an annotated single-cell reference atlas (\u003cstrong\u003eFigure 9A\u003c/strong\u003e). This analysis revealed that the majority of spatial spots across tumor slices were predominantly composed of MES-like malignant cells, followed by TAM-MG (tumor-associated microglia/macrophages) and AC-like tumor cells. Other subpopulations, including NPC-like, OPC-like, and immune cells such as T and B cells, were sparsely distributed. This spatial enrichment pattern suggests that MES-like cells represent the dominant malignant state within the sampled (central) tumor regions, consistent with their known invasive and therapy-resistant characteristics. We next performed a comprehensive niche analysis across all spatial spots to identify regionally coherent multicellular microenvironments. Based on compositional similarities, we defined seven distinct spatial niches (\u003cstrong\u003eFigure 10A\u003c/strong\u003e). Among them \u003cstrong\u003e(Figure 9B-C)\u003c/strong\u003e, niche3 and niche4 exhibited highly similar profiles characterized by enrichment of immune-related populations, including B cells, TAM-MG, T cells, and Endothelial cells, indicative of immune-dominant microenvironments with vascular involvement. In contrast, niche1 and niche5 were enriched for G2M-like, NPC-like, OPC-like, and Oligodendrocyte subtypes, reflecting proliferative and neural progenitor\u0026ndash;like tumor ecosystems with minimal immune infiltration, possibly representing stem-like or immunologically \u0026quot;cold\u0026quot; tumor niches. niche2 was defined by the high abundance of MES-like cells with moderate immune infiltration, suggesting a transitional state between immune-active and tumor-dominant environments. niche6 represented a highly tumor-centric ecosystem dominated by G2M-like and MES-like subtypes, with negligible contributions from other cell types. Finally, niche7 was characterized by enrichment of AC-like and OPC-like cells, potentially indicating a tumor region associated with glial differentiation programs. We further visualized the spatial distribution of high- and low-risk regions across these niches (\u003cstrong\u003eFigure 10B\u003c/strong\u003e), which revealed substantial heterogeneity in risk stratification within tumor tissues. Patient-level analyses of niche composition (\u003cstrong\u003eFigures 10C\u0026ndash;D\u003c/strong\u003e) demonstrated inter-individual differences in niche occupancy and subtype prevalence. In addition, spatial overlap between cellular subpopulations and niche types was assessed (\u003cstrong\u003eFigure 10E\u003c/strong\u003e), confirming niche-specific enrichment patterns of distinct cell types. Lastly, we examined the spatial expression patterns of three candidate therapeutic targets\u0026mdash;PDIA4, RFT1, and EIF3I\u0026mdash;originally identified from bulk multi-omics modeling. Consistent with scRNA-seq\u0026ndash;based findings, all three targets displayed diffuse and broadly distributed expression across tumor-dominant niches in the spatial transcriptomic data, without evident spatial restriction or compartmentalization (\u003cstrong\u003eFigures 10F\u0026ndash;G\u003c/strong\u003e). This pervasive expression pattern underscores their potential as pan-tumoral therapeutic targets within high-risk GBM regions. Collectively, these results demonstrate that GBM tissues harbor spatially segregated multicellular ecosystems with distinct tumor\u0026ndash;immune compositions, and that therapeutic targets identified through bulk profiling remain spatially relevant within tumor-dominant ecological contexts.\u003c/p\u003e"},{"header":"3. Discussion","content":"\u003cp\u003eThis study integrates multimodal data from GBM patients to develop a robust risk stratification model, enhancing the understanding of GBM heterogeneity across different risk groups. Key findings of our study include: (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) the development of a five-modal fusion risk prediction model that outperforms any single-modality or four-modality model, effectively stratifying patients into distinct risk groups with significant survival differences, offering valuable insights for personalized treatment strategies; (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) a comprehensive multi-omics analysis of GBM across different resolution levels, providing an in-depth characterization of GBM heterogeneity and identifying potential therapeutic targets, thereby presenting promising prospects for precision-targeted therapy. Several studies have explored GBM risk stratification, but most rely on single-omics or clinical data. While efforts like Yan et al(\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e)., Philipp et al(\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e). and Junseong et al(\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e). integrated clinical and radiomic features or transcriptomic data, they lacked multimodal integration. Some studies have utilized multimodal data, such as Wang et al(\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e)., who classified GBM subtypes using proteogenomic and metabolomic data, and Ravi et al. (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e), who examined tumor-host interactions through spatially resolved multimodal analysis. However, their clinical prognostic value remains limited. In contrast, our study integrates genomic, transcriptomic, proteomic, radiomic and pathomic data using a late-stage fusion strategy, enhancing prognostic accuracy and providing novel insights into tumor heterogeneity and potential therapeutic targets.\u003c/p\u003e\u003cp\u003eGenomic analysis of the high-risk group revealed specific mutations, including TP53, PI3KR1, ARHGEF18, and NEB, with no significant differences in mutation types or tumor mutation burden compared to the low-risk group. Further enrichment analysis highlighted the impact of pathways related to cell cycle regulation and gene expression. Copy number variation (CNV) analysis corroborated these findings, showing amplifications in genes linked to cell cycle regulation and metabolic dysregulation. These genomic alterations were consistent with transcriptomic and proteomic data, reinforcing the value of multi-omics integration for understanding the molecular heterogeneity of GBM. The immune microenvironment in the high-risk group was marked by increased immune pathway activation but a predominantly immunosuppressive state. This was characterized by extensive cancer-associated fibroblast (CAF) infiltration and recruitment of inflammatory cells (macrophages and neutrophils), suggesting immune tolerance mechanisms. These findings align with previous studies on the role of CAFs in immune suppression and poor prognosis(\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e). Single-cell transcriptomic analysis further revealed the enrichment of AC-like and NPC2-like cells, which are linked to immunosuppression and tumor progression(\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e). Pathway enrichment analysis highlighted the involvement of EGFR/ERBB2 signaling and P53 pathway abnormalities, contributing to immune evasion and tumor progression. In contrast, the low-risk group exhibited mutations in EGFR, PARP4, and CHRNB3, with enrichment in neuronal activity and molecular transport pathways. CNV analysis indicated a more stable genomic profile with fewer alterations in genes related to cell cycle regulation and gene expression control. The immune microenvironment in the low-risk group displayed severe T-cell dysfunction, but immune cell infiltration, particularly T-cells, was more pronounced than in the high-risk group(\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eBy combining quasi-bulk transcriptomes with multimodal risk models, we defined patient-specific risk groups and linked them to distinct cellular compositions. Notably, low-risk patients exhibited significant enrichment of oligodendrocyte-like subpopulations, whereas high-risk patients lacked any dominant non-malignant cell type, suggesting that preservation of glial-like differentiation may be associated with more favorable outcomes. This observation echoes prior reports highlighting the protective role of oligodendrocytic features in gliomas and reinforces the relevance of cell-type composition in prognostic evaluation(\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e). Gene expression and pathway enrichment analyses across malignant and non-malignant subpopulations revealed risk-specific transcriptional programs. High-risk cells exhibited elevated oxidative phosphorylation, translational activity, and nutrient stress pathways, particularly within G2M-like, AC-like, and MES-like tumor compartments. These features reflect metabolic reprogramming and proliferative signaling that may drive tumor aggressiveness. In contrast, low-risk counterparts preferentially activated homeostatic pathways, including stress response, immune signaling, and DNA repair. Immune subpopulations also mirrored this divergence: high-risk T cells and B cells exhibited exhaustion-like features and interferon stress, while low-risk immune cells retained cytokine signaling and immunocompetence. These results support the notion that risk stratification in GBM reflects a convergence of tumor-intrinsic metabolic pressure and immune dysfunction. To further link cell states with patient survival, we applied Scissor analysis, revealing that high-risk\u0026ndash;associated (Scissor⁺) cells were predominantly malignant, whereas low-risk\u0026ndash;associated (Scissor⁻) cells were mainly immune-derived. This reinforces the idea that tumor-intrinsic programs are the major drivers of poor prognosis, while immune competence may be associated with improved outcomes. Importantly, Scissor helped resolve the functional relevance of transcriptionally defined subpopulations, providing a bridge between bulk survival signals and cellular identities.\u003c/p\u003e\u003cp\u003eWe also assessed the expression patterns of three candidate therapeutic targets\u0026mdash;PDIA4, RFT1, and EIF3I\u0026mdash;identified via bulk multi-omics modeling. Single-cell and spatial transcriptomic analysis revealed that all three targets were highly expressed within tumor-dominant cell states, particularly G2M-like and MES-like subpopulations, with broad and diffuse distribution across tumor regions. These findings suggest that the proposed targets are not spatially restricted and may be amenable to pan-tumoral therapeutic strategies in high-risk GBM.\u003c/p\u003e\u003cp\u003eFinally, spatial deconvolution and niche analysis highlighted regionally coherent multicellular ecosystems with distinct tumor\u0026ndash;immune compositions. Among the seven identified spatial niches, some were immune-dominant (e.g., niche3 and niche4), others were stem-like or proliferation-driven (e.g., niche1 and niche5), and some were highly tumor-centric (e.g., niche6 and niche7). Risk stratification and target gene expression were non-uniformly distributed across these niches, further reinforcing the spatial dimension of prognostic heterogeneity.\u003c/p\u003e\u003cp\u003eDespite promising results, this study has limitations: (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) While external validation was performed for the single-omics prognostic models, the lack of publicly available multi-omics datasets precluded external validation of the multimodal fusion risk stratification model at total level. (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) During the development and validation of single-omics models, some models demonstrated high prognostic value in datasets derived from Asian populations but failed to exhibit sufficient prognostic significance in cohorts from TCGA, which is comprised of mainly non-Asian populations. (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e) The classifier used for risk mapping in single-cell and spatial transcriptomics analyses was trained on transcriptomic data, which may limit its ability to capture the full multimodal complexity of GBM. Future studies are needed to integrate more comprehensive multimodal datasets such as several external validation sets, and prospective validation sets to further corroborate the robustness of the model.\u003c/p\u003e"},{"header":"4. Materials and Methods","content":"\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e4.1 Experimental Design\u003c/h2\u003e\u003cp\u003eWe designed a multi-step integrative framework to investigate prognostic heterogeneity in glioblastoma (GBM) and identify potential therapeutic vulnerabilities. We collected a comprehensive multi-omics dataset for GBM patients, including genomic, transcriptomic, proteomic, radiomic, and histopathologic data, all matched with survival information. For each omics modality, multiple machine learning\u0026ndash;based survival algorithms were applied to construct individual prognostic models, generating modality-specific risk scores for each patient. These risk scores were subsequently integrated into a composite risk score using a Cox proportional hazards model\u0026ndash;based fusion strategy. Patients were then stratified into high-risk and low-risk groups according to the optimal cutoff of the composite risk score. To elucidate molecular differences between risk groups, we performed systematic heterogeneity analyses across multiple omics layers, including gene/protein expression profiles, mutation burden, pathway activities, and tumor microenvironment composition. Single-cell and spatial transcriptomic data were further incorporated to dissect the contributions of specific cell types and spatial niches to prognostic stratification. Based on the results of differential analyses, we prioritized candidate therapeutic targets that were consistently overregulated in high-risk patients, strongly associated with poor survival, and possessed druggable potential. These targets were then validated for cell-type specificity and spatial distribution at single-cell and spatial resolution, enhancing their translational applicability.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003e4.2 Data and Sample Collection\u003c/h2\u003e\u003cp\u003e This retrospective study collected data from patients with IDH-wildtype GBM who underwent surgical resection at FAHZZU between 2015 and 2021, in accordance with the 2021 WHO Classification of CNS Tumors. Inclusion criteria were: (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) age\u0026thinsp;\u0026ge;\u0026thinsp;18 years, (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) primary diffuse glioma, (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e) IDH-wildtype GBM diagnosis based on the 2021 WHO classification, (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e) no prior radiotherapy or chemotherapy, and (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e) complete clinical and follow-up data. Exclusion criteria included: (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) history of brain surgery or head trauma and (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) prior radiotherapy or chemotherapy.\u003c/p\u003e\u003cp\u003eA total of 1,187 patients with IDH-wildtype GBM were included (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). Of these, 1,122 had postoperative pathological specimens, from WSIs of H\u0026amp;E-stained sections were obtained. Preoperative MRI data were available for 936 patients, including T1-weighted, contrast-enhanced T1-weighted, T2-weighted, FLAIR, and ADC maps, all of high quality. Fresh surgical tumor specimens from 211 patients were snap-frozen and stored for tissue sequencing. Of these, 185 underwent RNA sequencing (RNA-seq), 167 underwent WES, 166 had proteomic data, 26 had scRNA-seq data, and 10 underwent spatial transcriptomic analysis. Preoperative MRI scans and postoperative slides were also collected from 108 and 125 GBM patients, respectively, at HPPH.\u003c/p\u003e\u003cp\u003e The study was approved by the Human Research Ethics Committees of Henan Provincial People's Hospital (Approval No. 2023\u0026thinsp;\u0026minus;\u0026thinsp;174) and FAHZZU (Approval Nos. 2019-KY-176 and 2023-KY-1028), with informed consent obtained from all patients for the use of fresh tumor specimens.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003e4.3 Collection and Processing of GBM Data from Public Databases\u003c/h2\u003e\u003cp\u003eIn this study, we aimed to collect as much sequencing data on GBM (GBM) as possible from public databases to validate our conclusions and enrich the research content. We obtained five GBM datasets from the following sources: The Cancer Genome Atlas (TCGA; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.cancer.gov/ccg/research/genome-sequencing/tcga\u003c/span\u003e\u003cspan address=\"https://www.cancer.gov/ccg/research/genome-sequencing/tcga\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), the Chinese Glioma Genome Atlas (CGGA; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.cgga.org.cn/\u003c/span\u003e\u003cspan address=\"http://www.cgga.org.cn/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), and the Clinical Proteomic Tumor Analysis Consortium (CPTAC; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://proteomics.cancer.gov/programs/cptac\u003c/span\u003e\u003cspan address=\"https://proteomics.cancer.gov/programs/cptac\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). These datasets included TCGA-GBM, CGGA325, CGGA693, and CPTAC-GBM. Additionally, we collected an external proteomics dataset, SMC-GBM, containing proteomic profiles of GBM patients from a previously published study(\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003e4.4 MRI Scanning and Imaging Feature Extraction\u003c/h2\u003e\u003cp\u003ePatient MRI images were acquired during routine examination using a 3.0 T MRI scanner (Siemens Magnetom Skyra/Trio TIM; GE Discovery MR750; Philips Ingenia). Sequences included: axial and sagittal T1-weighted imaging (T1WI), axial T2-weighted imaging (T2WI), axial T2-weighted fluid-attenuated inversion recovery (FLAIR) imaging, as well as axial, sagittal, and coronal post-contrast T1-weighted imaging (CE-T1WI) immediately after intravenous injection of a 0.1 mmol/kg dose of gadolinium-based contrast agent. Apparent diffusion coefficient (ADC) maps were obtained from axial diffusion-weighted imaging (DWI). The acquisition parameters for each sequence were as follows:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eT1WI and CE-T1WI: repetition time (TR) 220\u0026ndash;1750 ms; echo time (TE) 2.3\u0026ndash;24 ms; echo train length (ETL) 1\u0026ndash;12; slice thickness 5 mm; averages/excitations 1; flip angle (FA) 70\u0026deg;-111\u0026deg;; field of view (FOV) 220\u0026times;192\u0026ndash;240\u0026times;240 mm\u0026sup2;; matrix 256\u0026times;162\u0026ndash;320\u0026times;256 mm\u0026sup2;.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eT2WI: TR 1873\u0026ndash;5390 ms; TE 70\u0026ndash;117 ms; ETL 16\u0026ndash;32; slice thickness 5 mm; averages/excitations 1; FA 90\u0026deg;-142\u0026deg;; FOV 220\u0026times;192\u0026ndash;240\u0026times;240 mm\u0026sup2;; matrix 320\u0026times;238\u0026ndash;512\u0026times;512 mm\u0026sup2;.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eFLAIR: TR 4500\u0026ndash;8400 ms; TE 85\u0026ndash;150 ms; inversion time (TI) 1670\u0026ndash;2250 ms; ETL 1\u0026ndash;38; slice thickness 5 mm; averages/excitations 1; FA 90\u0026deg;-150\u0026deg;; FOV 220\u0026times;192\u0026ndash;240\u0026times;240 mm\u0026sup2;; matrix 256\u0026times;179\u0026ndash;256\u0026times;256 mm\u0026sup2;.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eDWI: Images were processed by the corresponding post-processing workstation, and ADC images were calculated from DWI acquired at b-values of 0 and 1000 s/mm\u0026sup2;. Sequence parameters included: TR 2121\u0026ndash;6000 ms; TE 77\u0026ndash;119 ms; ETL 1\u0026ndash;82; slice thickness 5 mm; averages/excitations 1; FA 90\u0026deg;; FOV 220\u0026times;220\u0026ndash;240\u0026times;240 mm\u0026sup2;; matrix 152\u0026times;114\u0026ndash;192\u0026times;192 mm\u0026sup2;. ADC maps for all imaging planes were generated on a voxel-by-voxel basis using a single-exponential model.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eFirst, the N4ITK algorithm was employed to correct bias field distortions for all sequences. After isotropic voxel resampling to 1\u0026times;1\u0026times;1 mm\u0026sup3; through trilinear interpolation, multi-sequence MRI rigid registration for each patient was performed using the axial resampled CE-T1WI as a template, and mutual information as similarity measure. This process was completed using the 3D Slicer software, generating registered images rT1WI, rCE-T1WI, rT2WI, rFLAIR, and rADC. Histogram matching was used for gray-level normalization on rT1WI, rCE-T1WI, rT2WI, and rFLAIR. We set the histogram level to 1024 and the number of matching points to 10 to achieve a finer match while preserving more details. A deputy chief physician in neuroradiology with over 10 years of experience in head MRI diagnosis manually delineated the tumor region of interest (ROI) on the axial plane of rFLAIR, rT2WI, and rCE-T1WI images using ITK-SNAP software, obtaining the tumor volume of interest (VOI)(\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e). The VOI was defined as the enhanced area, non-enhanced area, and necrotic area of the tumor. The VOI contour was drawn based on FLAIR images, while rT2WI and rCE-T1WI were used for cross-checking the tumor extent and fine-tuning the tumor contour. Z-score normalization was applied within the VOI for all sequences to adjust the ROI intensity to have a mean of 0 and a standard deviation of 1. This radiologist and a deputy chief physician in neurosurgery with over 10 years of work experience randomly selected 100 patients within the group for VOI redrawing using a simple random sampling method. Interclass correlation coefficients (ICC) were used to evaluate intra-rater reliability analysis for the test-retest dataset and inter-rater reliability analysis for the multiple description dataset, retaining features with ICC\u0026thinsp;\u0026ge;\u0026thinsp;0.75. The obtained VOI was then overlaid with co-registered rT1WI, rCE-T1WI, rT2WI, rFLAIR, and rADC images.\u003c/p\u003e\u003cp\u003ePyRadiomics was used to extract three categories of features, including first-order intensity statistics, shape descriptors, and higher-order texture features(\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e). Five basic matrices were employed to define texture features: the gray-level co-occurrence matrix (GLCM), gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), gray-level dependence matrix (GLDM), and neighborhood gray-tone difference matrix (NGTDM). In this study, imaging features were extracted from three types of images: original images, wavelet images, and Gaussian Laplace images. A PyRadiomics parameters file was provided in Github repository to enhance the reproducibility of feature extraction. Ultimately, 4788 features were extracted from the five MRI sequences, retaining 3015 features with ICC\u0026thinsp;\u0026ge;\u0026thinsp;0.8.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e\u003ch2\u003e4.5 Hematoxylin and Eosin (H\u0026amp;E) Histological Slide Scanning and Feature Analysis\u003c/h2\u003e\u003cp\u003ePathology slides were scanned at 20x magnification using a digital pathology scanner (KF-PRO-120-HI) to obtain the original WSIs. Subsequently, the original WSI underwent color space conversion, tissue segmentation, patch selection, and feature extraction. Specifically, the WSI at the 5x resolution was converted from RGB to Lab color space, and Otsu's algorithm was then applied to calculate a segmentation threshold for segmenting the tissue from WSI. The obtained tissue image was tiled into many 1024\u0026times;1024 patches at 20\u0026times; magnifications, where these patches were adjacent to one another covering the WSI. A Python package Yottixel was used to select the optimal patches for further analysis(\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e). Finally, CellProfiler (v4.2.5) software was used to extract features from each selected patch(\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e\u003ch2\u003e4.6 Whole Exome Sequencing (WES) and Analysis\u003c/h2\u003e\u003cp\u003eTumor tissue and adjacent brain tissue DNA were extracted from patients using the QIAamp Fast DNA Tissue Kit (Qiagen). Blood patients were collected in tubes containing EDTA and centrifuged at 1600 xg for 10 minutes at 4\u0026deg;C within 2 hours of collection. Peripheral blood lymphocyte (PBL) pellets were stored at -20\u0026deg;C until further use, and PBL DNA was extracted using the RelaxGene Blood DNA System (Tiangen Biotech Co., Ltd., Beijing, China). DNA quantification was performed using the Qubit 3.0 Fluorometer and Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Inc., Waltham, MA, USA). DNA collected from tissue and PBL patients were fragmented using dsDNA Fragmentase enzyme (New England BioLabs, Inc., Ipswich, MA, USA), followed by size selection of DNA fragments (150\u0026ndash;250 bp) using Ampure XP beads (Beckman Coulter, Inc., Brea, CA, USA). The KAPA Library Preparation Kit (Kapa Biosystems, Inc., Wilmington, MA, USA) was employed for the construction of DNA fragment libraries. Cleanup steps were performed using Agencourt AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA). After DNA fragmentation, end repair and 3' A-tailing were conducted, followed by exon capture using the Agilent SureSelect Human All Exon V6 kit. The Qubit 3.0 Fluorometer and Qubit dsDNA HS Assay Kit were utilized to assess the purity and concentration of DNA fragments. Fragment length was measured using the DNA 1000 kit (Agilent Technologies, Inc., Santa Clara, CA, USA) on a 4200 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). DNA libraries with 150 bp end sequences were sequenced using the Illumina Novaseq 6000 system. Raw data were converted to FASTQ files, and adapter and low-quality reads were trimmed using Trimmomatic (v0.39). We achieved a median coverage depth of 112x for tumor specimens and 128x for non-tumor specimens.\u003c/p\u003e\u003cp\u003eGATK (v4.2) tools were used to identify single nucleotide variants (SNVs) and insertions or deletions (INDELs). Paired-end WES reads were mapped to the human reference genome (hg38) using BWA-mem (v0.7.17). BAM files were further processed by reordering, sorting, marking duplicates, and adding read groups using Picard (v2.24.2). Base quality score recalibration was performed using the BaseRecalibrator module in GATK, followed by the assessment of cross-sample contamination using the GetPileupSummaries and CalculateContamination modules. Somatic variants were detected by MuTect2 and annotated using ANNOVAR, with patient-matched normal DNA sequencing reads serving as reference. Candidate somatic variants were distinguished based on the following filtering criteria: ① Variants outside of exonic regions and splice sites were excluded; ② Variants with a variant allele fraction (VAF)\u0026thinsp;\u0026ge;\u0026thinsp;5% and at least 2 supporting variant reads in tumor patients were retained; ③ Variants with a mutation allele frequency (MAF)\u0026thinsp;\u0026ge;\u0026thinsp;5% in at least one database, including 1000 Genomes, ESP6500, gnomAD, and ExAC, were removed. Normal patients were sequenced using the same scheme, each sample was reduced to 4% and then pooled as reference. To obtain high-quality and reliable somatic variants, we employed stringent downstream filtering criteria: ① Variants outside of exonic regions and splice sites were excluded; ② Variants with a VAF\u0026thinsp;\u0026ge;\u0026thinsp;5%, at least 5 supporting variant reads in tumor patients, and variants with a VAF in the tumor that was more than five times the VAF in the normal sample were retained; ③ Variants with more than 100 occurrences in COSMIC (v92) were retained; ④ Variants with a MAF\u0026thinsp;\u0026ge;\u0026thinsp;1% in at least one variant database (1000 Genomes, ESP6500, gnomAD, and ExAC) were removed; ⑤ Variants predicted as benign in at least two of the following tools: MutationAssessor, MutationTaster2, Polyphen2, and SIFT, were removed. Somatic CNVs were inferred by CNVkit (v0.9.9) based on BAM files generated during the somatic mutation detection process, using the default circular binary segmentation algorithm. Segment-level log2 ratios were calculated and transformed as input for the GISTIC 2.0 software to identify significantly amplified or deleted chromosomal regions in the tumors. CNV amplifications and deletions were defined using a\u0026thinsp;\u0026plusmn;\u0026thinsp;0.3 log2 ratio threshold.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e\u003ch2\u003e4.7 RNA Sequencing (RNA-seq) and Analysis\u003c/h2\u003e\u003cp\u003eTotal RNA was extracted from tissue patients using the TRIzol Reagent Kit (Ambion, Invitrogen, USA). RNA concentration and integrity were assessed using the Qubit RNA Assay Kit, Qubit 2.0 Fluorometer (Life Technologies), and Agilent 2100 Bioanalyzer (Agilent Technologies). Patients with an RNA integrity number greater than 5 were included in the study. Libraries were prepared from patients with high RNA integrity, no contaminants, and sufficient RNA quantity. Poly-T oligonucleotide magnetic beads were used to purify RNA from total RNA. RNA was fragmented in NEBNext First Strand Synthesis Reaction Buffer (5X) using divalent cations at elevated temperatures. cDNA synthesis, end repair, A-tailing, and NEBNext Adaptor ligation were performed using the NEBNext Ultra RNA Library Prep Kit. Library fragments were purified with AMPure XP (Beckman Coulter, Beverly, USA), selecting cDNA fragments of 150\u0026ndash;200 bp in length. Library quality was assessed using the Agilent Bioanalyzer 2100. Libraries were sequenced on the Illumina HiSeq X Ten platform, generating 150 bp paired-end reads. Sequencing data were filtered using Trimmomatic software to remove adaptors and low-quality sequences, followed by data quality assessment using FastQC. STAR (v2.7.6a) was used to align sequences to the reference genome (hg38). Gene expression values were calculated using RSEM (v1.3.3) based on the GENCODE (v35) gene annotation file. HTSeq v0.6.0 was used to count the number of reads aligned to each gene, and gene expression levels were quantified as FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) and TPM (Transcripts Per Kilobase of exon model per Million mapped reads).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec20\" class=\"Section2\"\u003e\u003ch2\u003e4.8 Mass Spectrometry\u003c/h2\u003e\u003cp\u003ePatients were removed from \u0026minus;\u0026thinsp;80\u0026deg;C storage, and an appropriate amount of tissue was weighed and placed into a liquid nitrogen pre-chilled mortar. Liquid nitrogen was added, and the tissue was thoroughly ground into a powder. Lysis buffer (1% Triton X-100, 1% protease inhibitor, 1% phosphatase inhibitor, 3 \u0026micro;M TSA, 50mM NAM) was added to each sample at four times the volume of the powder, followed by ultrasonic lysis. Patients were centrifuged at 4\u0026deg;C, 12,000 g for 10 minutes to remove cell debris, and the supernatant was transferred to a new centrifuge tube. Protein concentration was determined using a BCA assay kit.\u003c/p\u003e\u003cp\u003eEqual amounts of protein from each sample were digested with trypsin, and the volume was adjusted using lysis buffer. One volume of pre-chilled acetone was added, vortexed, and then four volumes of pre-chilled acetone were added, followed by precipitation at -20\u0026deg;C for two hours. Patients were centrifuged at 4,500 g for 5 minutes, and the supernatant was discarded. The pellet was washed twice with pre-chilled acetone. After air-drying the pellet, it resuspended in 200 mM TEAB, and trypsin was added at a 1:50 ratio (protease: protein, w/w) for overnight digestion. Dithiothreitol (DTT) was added to a final concentration of 5 mM, and patients were reduced at 56\u0026deg;C for 30 minutes. Iodoacetamide (IAA) was added to a final concentration of 11 mM, and patients were incubated in the dark at room temperature for 15 minutes.\u003c/p\u003e\u003cp\u003ePatients were separated using an Agilent 300Extend C18 column (4.6\u0026times;250 mm), with detection at 214 nm, a column oven temperature of 35\u0026deg;C, and 95% buffer A for 30 minutes to equilibrate the column. After baseline stabilization, the fractionation gradient method was initiated, and peptide patients were loaded onto the high-performance liquid chromatography (HPLC) fractionation column. Patients were collected at 1-minute intervals, with fractions 11 to 46 combined into 12 groups and vacuum-dried. Peptides were dissolved in mobile phase A and separated using an EASY-nLC 1200 ultra-HPLC system. Mobile phase A consisted of a 0.1% formic acid and 2% acetonitrile aqueous solution, while mobile phase B consisted of a 0.1% formic acid and 90% acetonitrile aqueous solution. The gradient was set as follows: 0\u0026ndash;96 min, 6%-25% B; 96\u0026ndash;114 min, 25%-35% B; 114\u0026ndash;117 min, 35%-80% B; and 117\u0026ndash;1200 min, 80% B, with a flow rate maintained at 500 nl/min. Separated peptides were ionized in the NSI ion source and data were collected using the Orbitrap Exploris 480 mass spectrometer.\u003c/p\u003e\u003cp\u003eLiquid chromatography (LC) parameters were consistent with those used during library construction. Peptides were separated using the ultra-high-performance liquid chromatography system and analyzed using the Orbitrap Exploris 480 mass spectrometer. Precursor ions and their fragment ions were detected and analyzed using the high-resolution Orbitrap. FAIMS compensation voltage (CV) settings were \u0026minus;\u0026thinsp;40 V, -55 V, and \u0026minus;\u0026thinsp;70 V. The primary mass scan range was set at 350\u0026ndash;1350 m/z with a resolution of 120,000; the secondary scan resolution was set at 30,000. The secondary data acquisition mode was set to DIA mode, which followed a primary scan with 20 m/z window peptide ions entering the HCD collision cell using 32% collision energy for fragmentation, and subsequent secondary mass analysis. The automatic gain control (AGC) for the secondary spectrum was set at 600%.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec21\" class=\"Section2\"\u003e\u003ch2\u003e4.9 Development and Validation of a Risk Stratification Model Based on Multimodal Data Integration\u003c/h2\u003e\u003cp\u003eIn this study, the multimodal GBM dataset from the FAHZZU was used as the training set for each omics layer, while other GBM datasets served as external validation cohorts. At each omics level, we implemented a comprehensive analytical framework utilizing the R package Mime1, which integrates 10 distinct machine learning algorithms along with their 101 potential combinations(\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e). The combination that achieved the highest average concordance index (C-index) in the external validation sets was selected as the optimal algorithm for that omics layer, which was then used to calculate patient-specific risk scores. Subsequently, utilizing a late-fusion strategy combined with 10-fold cross-validation, we integrated the risk scores from different omics layers through Cox proportional hazards regression to derive a comprehensive multimodal fusion risk score. Based on the optimal cutoff value, patients were stratified into high- and low-risk groups, thereby establishing the multimodal fusion-based risk stratification model for GBM(Figure \u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec22\" class=\"Section2\"\u003e\u003ch2\u003e4.10 Genomic Alteration Analysis\u003c/h2\u003e\u003cp\u003eUsing the maftools package, WES data was processed(\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e). To investigate the differences in gene mutation profiles between the two risk groups, we conducted a comprehensive analysis encompassing multiple mutation categories. Initially, the Wilcoxon rank-sum test was applied to assess variations in gene mutation types (frameshift, non-frameshift, point mutations, splice site mutations, and translation initiation mutations), variant classifications (deletions, insertions, multinucleotide polymorphisms, and single nucleotide polymorphisms), and single base substitution types (C\u0026thinsp;\u0026gt;\u0026thinsp;A, C\u0026thinsp;\u0026gt;\u0026thinsp;T, T\u0026thinsp;\u0026gt;\u0026thinsp;A, C\u0026thinsp;\u0026gt;\u0026thinsp;G, T\u0026thinsp;\u0026gt;\u0026thinsp;C, and T\u0026thinsp;\u0026gt;\u0026thinsp;G). Building upon these findings, we calculated the mutation frequencies of all genes across the risk groups and visualized genes with mutation frequencies exceeding 5% using maftools. To elucidate the functional implications of these mutated genes, over-representation analysis (ORA) was performed with the ClusterProfiler package(\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e), retaining pathways with p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05 as statistically significant.\u003c/p\u003e\u003cp\u003eSubsequently, to identify group-specific tumor driver genes, we utilized the OncodriveCLUST algorithm, visualizing the results with maftools and performing ORA with ClusterProfiler to determine enriched pathways associated with these driver genes. In parallel, co-occurring and mutually exclusive mutation patterns between the risk groups were explored using Fisher's exact test. The resulting gene interaction networks were visualized via the GGally package, providing insights into distinct mutational landscapes, which were further analyzed through ORA to identify significant pathways (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05).\u003c/p\u003e\u003cp\u003eTo comprehensively assess the mutational impact at the pathway level, we calculated pathway mutation burdens based on gene sets from the Reactome, KEGG, and Hallmark databases using maftools. The mutation burden across pathways was visualized with ggplot2, with particular attention to those affecting more than 10% of patients in each group. Finally, pathways identified from the aforementioned analyses were systematically annotated and classified. This integrative approach enabled a detailed evaluation of how gene mutations influence key biological processes, providing mechanistic insights into the differential risk stratification between the two groups.\u003c/p\u003e\u003cp\u003eFurthermore, to explore chromosomal-level alterations, we performed a differential analysis of copy number variations (CNVs) between the high- and low-risk groups to identify significantly altered chromosomal segments. First, significantly different CNV regions were identified based on GISTIC 2.0(\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e) analysis results, with statistical comparisons between groups conducted using the limma package. Following the identification of these regions, functional analysis was performed to elucidate their biological relevance. Genes located within the significantly altered chromosomal segments were extracted and subjected to ORA using ClusterProfiler, taking into account both amplification and deletion statuses. Pathways with p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05 were considered significant, enabling the identification of affected biological pathways and their potential functional consequences associated with CNV differences across risk groups.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec23\" class=\"Section2\"\u003e\u003ch2\u003e4.11 Functional Analysis Based on Transcriptomic and Proteomic Profiles\u003c/h2\u003e\u003cp\u003eTo investigate transcriptomic and proteomic differences between the high- and low-risk groups, we performed a comprehensive differential expression analysis to identify significantly altered genes and proteins. First, using the Wilcoxon rank-sum test, we applied the limma package to compare samples from different risk groups, selecting genes and proteins with a false discovery rate (FDR)\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and an absolute log₂ fold change (|log₂FC|)\u0026thinsp;\u0026gt;\u0026thinsp;1 as significantly differentially expressed.\u003c/p\u003e\u003cp\u003eBased on these differentially expressed genes (DEGs) and proteins, we conducted over-representation analysis (ORA) using the ClusterProfiler package to identify enriched pathways with p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05. The enriched pathways were subsequently clustered and visualized using the aPEAR package(\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e), allowing for a clearer understanding of functional groupings. To further explore biological processes associated with expression changes, we performed GSEA. DEGs were ranked according to their log₂FC values, and GSEA was carried out with ClusterProfiler, identifying pathways with FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05. The results were clustered and visualized using aPEAR, providing insights into pathway-level alterations between the risk groups.\u003c/p\u003e\u003cp\u003eIn addition to ORA and GSEA, we employed gene set variation analysis (GSVA) to assess pathway activity differences across patient transcriptomes. By leveraging gene sets from the Reactome, KEGG, and Hallmark databases, GSVA scores were calculated and filtered based on reliable pathways identified through GSEA. These scores were then visualized using the ComplexHeatmap package(\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e), facilitating the annotation and comparison of relevant pathways between the two risk groups.\u003c/p\u003e\u003cp\u003eThrough the integration of differential expression analysis, pathway enrichment, and variation analyses, this comprehensive approach provides a multidimensional understanding of the molecular differences underlying risk stratification in GBM patients.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec24\" class=\"Section2\"\u003e\u003ch2\u003e4.12 Tumor Microenvironment (TME) Analysis Based on Transcriptomic Data\u003c/h2\u003e\u003cp\u003eBased on the stratification of high and low-risk groups, we inferred the tumor immune microenvironment from the transcriptomic data and identified immune activity-related indicators that exhibited significant differences between the groups. First, we performed a differential analysis of immune-related pathways by selecting immune pathways from the GSEA results, categorizing their functions, and visualizing them using the R package ggplot2 to examine the differences in pathway activation between the two groups. Second, we evaluated the immune and stromal cell infiltration in tumor tissues using the ESTIMATE(\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e) algorithm and calculated the immune and stromal scores. We then conducted a Wilcoxon rank-sum test to analyze the differences in scores between the two risk groups. Third, we assessed the immune phenotypes, including antigen presentation (MHC-MHC molecules), effector cells (EC-effector cells), suppressor cells (SC-suppressor cells), and immune checkpoints (CP-immune checkpoints) from the transcriptomic data and calculated the scores for each immune activity(\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e). A Wilcoxon rank-sum test was performed to examine the differences in immune scores between the two risk groups. Fourth, we analyzed the Tumor Immune Dysfunction and Exclusion(\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e) by predicting the cancer patients' responses to immune checkpoint inhibitors (such as PD-1/PD-L1 inhibitors) and other immune-related parameters. Differences in immune-related indices between the two groups were assessed using the Wilcoxon rank-sum test. Finally, we performed a differential analysis of the Cancer Immune Cycle(\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e) by predicting the anti-cancer immune states at seven immune stages and calculating the corresponding scores. The Wilcoxon rank-sum test was applied to analyze the differences in anti-cancer immune state scores between the two risk groups.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec25\" class=\"Section2\"\u003e\u003ch2\u003e4.13 Heterogeneity Analysis Between Single-Cell and Spatial Transcriptomic Data\u003c/h2\u003e\u003cp\u003eThe R package seurat(\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e) was used to perform preprocessing such as single-cell data quality control and annotation to determine the cell subtypes contained in the data, and the RCTD algorithm(\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e) was used to map the single-cell subtypes to the patient's spatial transcriptome data. Based on the multimodal fusion risk stratification model, single-cell and spatial transcriptomic data were categorized into risk groups, and their heterogeneity was comprehensively analyzed. First, the multimodal cohort was divided into a training set and a validation set in a 7:3 ratio. In the training set, logistic regression and ROC analysis were used to select genes with P\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and AUC\u0026thinsp;\u0026gt;\u0026thinsp;0.7, followed by Lasso regression for feature selection and dimensionality reduction. Genes with non-zero coefficients were retained for model construction, and a classifier was built by combining 14 machine learning algorithms (e.g., AdaBoost, decision tree, GBDT, SVM, XGBoost). The classifier was applied to pseudo-bulk data from the single-cell transcriptomic cohort to assign risk groups. Next, the Wilcoxon rank-sum test was employed to analyze differences in cell subpopulation proportions between the two groups, identifying significantly different subpopulations. Differential gene expression analysis was then performed for the cell subpopulations, selecting genes with FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and absolute log2FC\u0026thinsp;\u0026gt;\u0026thinsp;2, followed by enrichment analysis using ClusterProfiler to identify significant pathways (P\u0026thinsp;\u0026lt;\u0026thinsp;0.05) and analyze functional characteristics of the same subpopulations across risk groups. Additionally, we employed the R package Scissor to integrate single-cell transcriptomic data with bulk RNA-seq and survival information, enabling the identification of single-cell subpopulations associated with patient prognosis(\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e). This approach allowed us to pinpoint individual cells whose transcriptional profiles are positively or negatively correlated with clinical outcomes.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec26\" class=\"Section2\"\u003e\u003ch2\u003e4.14 Identification of Therapeutic Targets Based on Multi-Omics Data\u003c/h2\u003e\u003cp\u003eBased on multi-modal risk stratification, we identified EIF3I as a potential therapeutic target associated with the high-risk group. Firstly, differential analyses were conducted at the copy number variation, transcriptomic, and proteomic levels according to the multi-modal risk stratification, and the characteristic genes of the high-risk group were selected, followed by identification of the common genes across these levels. These genes were then subjected to correlation analysis between copy number variation-transcriptome and transcriptome-proteome data, which helped identify genes with high correlations across different modalities. Subsequently, these genes were analyzed in the TCGA and CGGA databases to assess differential expression between GBM and low-grade gliomas at the transcriptomic level, selecting those with high expression in GBM. Finally, survival analysis was performed on these genes across multiple databases at various omics levels and based on the Log-rank P-value and univariate Cox P-value (both \u0026lt;\u0026thinsp;0.05), potential therapeutic targets related to the high-risk group were screened and ranked.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec27\" class=\"Section2\"\u003e\u003ch2\u003e4.15 Acquisition of GBM Single-cell RNA-seq Data\u003c/h2\u003e\u003cp\u003eTissue Dissociation and Single-cell Suspension Preparation: The tissue was transferred to a culture dish containing 1\u0026times; PBS (without RNase and Ca/Mg ions), placed on ice, and washed to remove surface blood clots and fats. The tissue was then minced into 0.5 mm\u0026sup2; pieces and washed again with PBS. Dissociation solvent (0.35% collagenase IV, 2 mg/ml papain, 120 Units/ml DNase I) was added, and the mixture was incubated in a 37\u0026deg;C water bath shaker (100 rpm) for 20 minutes. The reaction was terminated by adding PBS containing 10% fetal bovine serum, followed by gentle pipetting 5\u0026ndash;10 times. The cell suspension was filtered through a 70\u0026thinsp;\u0026minus;\u0026thinsp;30 \u0026micro;m cell strainer, and centrifuged at 300g for 5 minutes at 4\u0026deg;C to collect the cell pellet. The pellet was resuspended in 100 \u0026micro;l of 1\u0026times; PBS (0.04% BSA) solution, followed by the addition of 1 ml of 1\u0026times; red blood cell lysis buffer. The cells were incubated at room temperature or on ice for 2\u0026ndash;10 minutes. Afterward, cells were centrifuged for 5 minutes, and the pellet was resuspended in 100 \u0026micro;l of Dead Cell Removal MicroBeads reagent and incubated for 15 minutes at room temperature. Dead cells were removed using MS Columns, followed by another 5-minute centrifugation to collect the viable cells, which were then resuspended in 1\u0026times; PBS (0.04% BSA). After two centrifugation steps at 4\u0026deg;C, the final viable cells were obtained. Cell viability was assessed using trypan blue staining (\u0026gt;\u0026thinsp;85%) and the cell count was determined, with a concentration of 700\u0026ndash;1200 cells/\u0026micro;l. Chromium 10x Genomics Library Construction and Sequencing: According to the protocol of the 10X Genomics Chromium Single-Cell 3' Kit (V3), the single-cell suspension was loaded onto the 10x Chromium chip, with an expected capture of 8,000 single cells. cDNA amplification and library construction were performed according to standard protocols, and the libraries were sequenced on an Illumina NovaSeq 6000 system with paired-end sequencing (150 bp) at a depth of at least 20,000 reads per cell. The operation was performed by LC-Bio Technology (Hangzhou, China).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec28\" class=\"Section2\"\u003e\u003ch2\u003e4.16 Acquisition of GBM Spatial Transcriptomics Data\u003c/h2\u003e\u003cp\u003eOCT-embedded tissue blocks were sectioned into 10 \u0026micro;m thick slices measuring 6.5 mm x 6.5 mm, which were then adhered to Visium slides. Hematoxylin and eosin (H\u0026amp;E) staining was performed according to the 10x Genomics Visium Fresh Frozen Tissue Processing Protocol. Fluorescent stitching microscopy (Olympus Fluoview 1000) was used to capture H\u0026amp;E images. The tissue was removed, and library construction was carried out according to the 10x Genomics recommended protocol. Each spatial transcriptomics library was processed using the 10x Genomics Space Ranger software (version 2.0.0) and aligned to the reference genome (mm10 for mouse/GRCh38 for human). UMI counts were compiled for each probe site. To distinguish probe sites covering the tissue from background, only barcodes associated with the tissue-covered regions were retained. The filtered UMI count matrix was generated, followed by manual exclusion of probe sites that were not covered by tissue and undetected by Space Ranger. The UMI count matrix was further optimized.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec29\" class=\"Section2\"\u003e\u003ch2\u003e4.17 Experimental verification of biological functions of the targets\u003c/h2\u003e\u003cp\u003e(\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) Cell culture: The GBM cell lines U251(RRID:CVCL_0021), A172 (RRID:CVCL_0131) and U118(RRID:CVCL_0633) were obtained from Fenghui Biotechnology Co.Ltd. Cells were cultured in Dulbecco's modified Eagle's Medium (DMEM, Sigma) supplemented with 10% fetal bovine serum (FBS, Gibco) at 37.0\u0026deg;C under a humidified 5.0% CO2 atmosphere. All cell lines were tested for mycoplasma every 2 months.All cell lines were contamination free. (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) Cell transfection: All siRNAs and pLKO.1-puro shRNA plasmids (RRID:Addgene_8453) were synthesized by Tsingke Biotechnology Co.Ltd. The sequences of shRNAs and siRNAs are summarized in \u003cb\u003eData S2\u003c/b\u003e. Plasmids and siRNAs transfection was performed with jetPRIME reagent (Polyplus) according to the manufacturer\u0026rsquo;s instructions. Transfection efficiency was determined by reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR). (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e) RT-qPCR: Total RNA was extracted using Trizol (ABclonal), and 1000 ng of RNA was reversed transcribed into cDNA using the RT Master Mix for qPCR II (MCE). And qPCR was performed using the 2X Universal SYBR Green Fast qPCR Mix (ABclonal) and the CFX Opus 384 Real-Time PCR System (BioRad). Gene expression values were normalized to GADPH, which was used as an internal control. The primer sequences used in this study are listed in \u003cb\u003eData S2\u003c/b\u003e. (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e) Cell viability assay: Cell viability was measured using Cell Counting Kit 8 (CCK-8, Servicebio) according to the manufacturer's instructions. At 0, 24, 48, and 72h, absorbance was measured at 450 nm by using a microplate reader to estimate the cell viability. All experiments were performed in biological triplicate. (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e) EdU staining assay: The EdU assay was performed according to the manufacturer's instructions (BeyoClick\u0026trade; EdU Cell Proliferation Kit with Alexa Fluor 555; Beyotime, China). Briefly, cells were cultured in 24-well plates, treated with 500\u0026micro;L of medium containing 10 \u0026micro;M EdU at 37.0\u0026deg;C under a 5.0% CO2 atmosphere for 2h, fixed with 4% paraformaldehyde for 20 min, incubated with PBS containing 0.3% Triton X‐100 for 10 min, and treated with Click Reaction Solution for 30 min. The nuclei were counterstained with Hoechst 33342 for 10 min. Three randomly selected regions from each group were visualized using a fluorescence microscope (Leica, Wetzlar, Germany). (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e) Wound‐healing assay: GBM cells from different transfection groups were seeded in 6‐well plates at 1 \u0026times; 106 cells/well and grown overnight to form monolayers. A sterile tip was used to create a wound line (cell‐free area) across the culture plate surface, and the deplated cells were removed by washing with PBS. Then, cells were cultured in FBS‐free DMEM under a humidified 5.0% CO2 atmosphere at 37.0\u0026deg;C for 24 hours, and images of the wound lines were acquired using a phase‐contrast microscope (Zeiss, Germany). Each assay was performed in triplicate. The migration ability of GBM cell was estimated by measuring the scratch width. (\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e) Transwell invasion assay: For the invasion assay, the Transwell upper chambers (Corning) were pre‐coated with Matrigel (Corning) before cell seeding. After serum starvation for 24 hours, GBM cells from different transfection groups were seeded in the upper chambers at 105 cells/200 \u0026micro;L serum‐free DMEM, while 600\u0026micro;L of DMEM with 20% FBS was added to the lower chambers. After 24 hours, all non‐invaded cells on the filter side of the upper chamber were removed by using a cotton swab, and the polycarbonate membrane of the upper chamber was fixed with 4% paraformaldehyde for 30 min, rinsed three times with PBS, and stained with 0.3% crystal violet for 10 min. Invaded cells were counted under a microscope (Zeiss, Germany). All Transwell assays were repeated at least thrice.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec30\" class=\"Section2\"\u003e\u003ch2\u003e4.18 Statistical Analysis\u003c/h2\u003e\u003cp\u003eAll data processing, statistical analyses, and graphing were performed using R software (version 4.2.2). Normality was assessed using the Shapiro-Wilk test. The Wilcoxon rank-sum test was used to compare continuous variables between two groups. The correlation between two continuous variables was evaluated using Spearman or Pearson correlation coefficient. Survival analysis and Kaplan-Meier curve plotting were performed using the survival and survminer packages. The Benjamini-Hochberg method was used to correct the FDR value obtained from multiple comparisons of \u003cem\u003eP\u003c/em\u003e values. All statistical tests were two-sided, with \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05 indicating statistical significance.\u003c/p\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConceptualization: YDM, ZYZ, YCJ, ZCL\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eData curation: ZLW, ZYM, WCD, MKW, WWW, JY, HRL, WYL, YHY, TC, CYM, MMY, JF, SL, BY, YWZ, RKC, DYS, GY\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eFormal analysis: ZYZ, ZLW, ZQL,\u0026nbsp;QCS, YSZ,\u0026nbsp;JXD\u003c/p\u003e\n\u003cp\u003eProject administration: YDM, ZYZ, YCJ, XZL, ZCL, HRZ, DL\u003c/p\u003e\n\u003cp\u003eWriting \u0026ndash; original draft: ZYZ, ZLW, RL, DLP, JDL, YNQ\u003c/p\u003e\n\u003cp\u003eWriting \u0026ndash; review \u0026amp; editing: ZYZ, ZLW\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAuthors declare that they have no competing interests\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability Statement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe raw WES and RNA-seq data generated in this study have been deposited in the Genome Sequence Archive (GSA) database under accession code HRA006184, accessible at [https://ngdc.cncb.ac.cn/gsa-human/browse/HRA006184]. Additionally, the raw MS-based proteomics data are deposited in the Open Archive for Miscellaneous Data (OMIX) database under accession code OMIX005975, available at [https://ngdc.cncb.ac.cn/omix/release/OMIX005975]. Due to data privacy laws, the raw radiomics and pathomics data are protected and are not available. However, the processed WES, RNA-seq, MS-proteomics, radiomics, pathomics, scRNA-seq and ST data can be accessed via Zenodo at [https://doi.org/10.5281/zenodo.14897960]. For all other data from this study, inquiries can be made to the corresponding authors upon reasonable request.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding:\u0026nbsp;\u003c/strong\u003eNational Natural Science Foundation of China grant 82273493 (ZYZ)\u003c/p\u003e\n\u003cp\u003eNational Natural Science Foundation of China grant U20A2017 (ZCL)\u003c/p\u003e\n\u003cp\u003eNational Natural Science Foundation of China grant 82173090 (XZL)\u003c/p\u003e\n\u003cp\u003eNational Natural Science Foundation of China grant U1904148 (XZL)\u003c/p\u003e\n\u003cp\u003eScience and Technology Program of Henan Province grant 242102311107 (YCJ)\u003c/p\u003e\n\u003cp\u003eScience and Technology Research and Development Joint Fund of Henan Province grant 242301420014 (ZYZ)\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eQ. T. Ostrom et al., CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2016-2020. Neuro Oncol 2023;25(12 Suppl 2):iv1-iv99.\u003c/li\u003e\n\u003cli\u003eH. Yan et al., IDH1 and IDH2 mutations in gliomas. N Engl J Med 2009;360(8):765-773.\u003c/li\u003e\n\u003cli\u003eD. N. Louis et al., The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro Oncol 2021;23(8):1231-1251.\u003c/li\u003e\n\u003cli\u003eM. J. van den Bent et al., Primary brain tumours in adults. Lancet 2018;392(10145):432-446.\u003c/li\u003e\n\u003cli\u003eR. Mathur et al., Glioblastoma evolution and heterogeneity from a 3D whole-tumor perspective. Cell 2024;187(2):446-463.e16.\u003c/li\u003e\n\u003cli\u003eL. B. Wang et al., Proteogenomic and metabolomic characterization of human glioblastoma. Cancer Cell 2021;39(4):509-528.e20.\u003c/li\u003e\n\u003cli\u003eK. H. Kim et al., Integrated proteogenomic characterization of glioblastoma evolution. Cancer Cell 2024;42(3):358-377.e8.\u003c/li\u003e\n\u003cli\u003eH. S. Phillips et al., Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell 2006;9(3):157-173.\u003c/li\u003e\n\u003cli\u003eR. G. Verhaak et al., Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell 2010;17(1):98-110.\u003c/li\u003e\n\u003cli\u003eM. Ceccarelli et al., Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma. Cell 2016;164(3):550-563.\u003c/li\u003e\n\u003cli\u003eM. Wu et al., Multi-omics and pharmacological characterization of patient-derived glioma cell lines. Nat Commun 2024;15(1):6740.\u003c/li\u003e\n\u003cli\u003eQ. Wang et al., Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell 2018;33(1):152.\u003c/li\u003e\n\u003cli\u003eK. White et al., Identification, validation and biological characterisation of novel glioblastoma tumour microenvironment subtypes: implications for precision immunotherapy. Ann Oncol 2023;34(3):300-314.\u003c/li\u003e\n\u003cli\u003eF. Guan et al., Biological underpinnings of radiomic magnetic resonance imaging phenotypes for risk stratification in IDH wild-type glioblastoma. J Transl Med 2023;21(1):841.\u003c/li\u003e\n\u003cli\u003eX. Wu et al., Big data and artificial intelligence in cancer research. Trends Cancer 2024;10(2):147-160.\u003c/li\u003e\n\u003cli\u003eN. Jamshidi et al., Illuminating radiogenomic characteristics of glioblastoma multiforme through integration of MR imaging, messenger RNA expression, and DNA copy number variation. Radiology 2014;270(1):1-2.\u003c/li\u003e\n\u003cli\u003eF. Carrillo-Perez et al., Applications of artificial intelligence in dentistry: A comprehensive review. J Esthet Restor Dent 2022;34(1):259-280.\u003c/li\u003e\n\u003cli\u003eJ. Lipkova et al., Artificial intelligence for multimodal data integration in oncology. Cancer Cell 2022;40(10):1095-1110.\u003c/li\u003e\n\u003cli\u003eD. Tamborero et al., OncodriveCLUST: exploiting the positional clustering of somatic mutations to identify cancer genes. Bioinformatics 2013;29(18):2238-2244.\u003c/li\u003e\n\u003cli\u003eM. Kanehisa, S. Goto, KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28(1):27-30.\u003c/li\u003e\n\u003cli\u003eA. Liberzon et al., The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1(6):417-425.\u003c/li\u003e\n\u003cli\u003eM. Milacic et al., The Reactome Pathway Knowledgebase 2024. Nucleic Acids Res 2024;52(D1):D672-D678.\u003c/li\u003e\n\u003cli\u003eM. Reich et al., GenePattern 2.0. Nat Genet 2006;38(5):500-501.\u003c/li\u003e\n\u003cli\u003eJ. Yan et al., Image-based deep learning identifies glioblastoma risk groups with genomic and transcriptomic heterogeneity: a multi-center study. Eur Radiol 2023;33(2):904-914.\u003c/li\u003e\n\u003cli\u003eP. Karschnia et al., Development and validation of a clinical risk model for postoperative outcome in newly diagnosed glioblastoma: a report of the RANO resect group. Neuro Oncol 2025;27(4):1046-1060.\u003c/li\u003e\n\u003cli\u003eJ. Park et al., Transcriptome profiling-based identification of prognostic subtypes and multi-omics signatures of glioblastoma. Sci Rep 2019;9(1):10555.\u003c/li\u003e\n\u003cli\u003eV. M. Ravi et al., Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma. Cancer Cell 2022;40(6):639-655.e13.\u003c/li\u003e\n\u003cli\u003eGalbo PM Jr et al., Molecular Features of Cancer-associated Fibroblast Subtypes and their Implication on Cancer Pathogenesis, Prognosis, and Immunotherapy Resistance. Clin Cancer Res 2021;27(9):2636-2647.\u003c/li\u003e\n\u003cli\u003eC. Neftel et al., An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell 2019;178(4):835-849.e21..\u003c/li\u003e\n\u003cli\u003eK. Woroniecka et al., T-Cell Exhaustion Signatures Vary with Tumor Type and Are Severe in Glioblastoma. Clin Cancer Res 2018;24(17):4175-4186.\u003c/li\u003e\n\u003cli\u003eF. S. Ng et al., Progenitor-like traits contribute to patient survival and prognosis in oligodendroglial tumors. Clin Cancer Res 2012;18(15):4122-4135.\u003c/li\u003e\n\u003cli\u003eS. Oh et al., Integrated pharmaco-proteogenomics defines two subgroups in isocitrate dehydrogenase wild-type glioblastoma with prognostic and therapeutic opportunities. Nat Commun 2020;11(1):3288.\u003c/li\u003e\n\u003cli\u003eP. A. Yushkevich et al., User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage 2006;31(3):1116-1128.\u003c/li\u003e\n\u003cli\u003eJ. J. M. van Griethuysen et al., Computational Radiomics System to Decode the Radiographic Phenotype. Cancer Res 2017;77(21):e104-e107.\u003c/li\u003e\n\u003cli\u003eS. Kalra et al., Yottixel - An Image Search Engine for Large Archives of Histopathology Whole Slide Images. Med Image Anal 2020;65:101757.\u003c/li\u003e\n\u003cli\u003eD. R. Stirling et al., CellProfiler 4: improvements in speed, utility and usability. BMC Bioinformatics 2021;22(1):433.\u003c/li\u003e\n\u003cli\u003eH. Liu et al., Mime: A flexible machine-learning framework to construct and visualize models for clinical characteristics prediction and feature selection. Comput Struct Biotechnol J 2024;23:2798-2810.\u003c/li\u003e\n\u003cli\u003eA. Mayakonda et al., Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28(11):1747-1756.\u003c/li\u003e\n\u003cli\u003eS. Xu et al., Using clusterProfiler to characterize multiomics data. Nat Protoc 2024;19(11):3292-3320.\u003c/li\u003e\n\u003cli\u003eI. Kerseviciute, J. Gordevicius, aPEAR: an R package for autonomous visualization of pathway enrichment networks. Bioinformatics 2023;39(11):btad672.\u003c/li\u003e\n\u003cli\u003eZ. Gu, Complex heatmap visualization. Imeta 2022;1(3):e43.\u003c/li\u003e\n\u003cli\u003eK. Yoshihara et al., Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612.\u003c/li\u003e\n\u003cli\u003eP. Charoentong et al., Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18(1):248-262.\u003c/li\u003e\n\u003cli\u003eP. Jiang et al., Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24(10):1550-1558.\u003c/li\u003e\n\u003cli\u003eMellman I et al., The cancer-immunity cycle: Indication, genotype, and immunotype. Immunity 2023;56(10):2188-2205.\u003c/li\u003e\n\u003cli\u003eY. Hao et al., Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42(2):293-304.\u003c/li\u003e\n\u003cli\u003eD. M. Cable et al., Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol 2022;40(4):517-526.\u003c/li\u003e\n\u003cli\u003eD. Sun et al., Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol 2022;40(4):527-538.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"molecular-cancer","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"molc","sideBox":"Learn more about [Molecular Cancer](http://gsejournal.biomedcentral.com/)","snPcode":"12943","submissionUrl":"https://submission.nature.com/new-submission/12943/3","title":"Molecular Cancer","twitterHandle":"@SN_Oncology","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Multi-Omics, Risk Stratification, Molecular Heterogeneity, Glioblastoma","lastPublishedDoi":"10.21203/rs.3.rs-7510857/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7510857/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eMultimodal data integration reveals causal features often missed by single-modality analyses, offering a more comprehensive view of glioblastoma (GBM) complexity. We collected radiomic, pathomic, genomic, transcriptomic, and proteomic data from patients with IDH-wild-type GBM to construct a machine learning\u0026ndash;based risk stratification model. While sample sizes varied across modalities, 147 patients with complete data across all five omics layers were used for integrative analysis. This approach identified two clinically distinct subgroups. The low-risk group, linked to favorable outcomes, showed enhanced neurodevelopmental signatures, increased neuronal infiltration, and more oligodendrocytes. In contrast, the high-risk group, associated with poor prognosis, exhibited strong proliferative signals and hyperactive cell cycle pathways. Downstream multi-omics analysis identified PDIA4, EIF3I, and RFT1 as potential prognostic biomarkers and therapeutic targets in high-risk GBM. These findings underscore the utility of multimodal machine learning in refining prognostic models, characterizing tumor heterogeneity, and informing personalized treatment strategies.\u003c/p\u003e","manuscriptTitle":"Machine Learning Model on Multi-Omics Data Enables Risk Stratification and Identifies Molecular Heterogeneity and Therapeutic Targets in Glioblastoma","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-09-12 14:27:11","doi":"10.21203/rs.3.rs-7510857/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-10-07T10:57:58+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-09-30T22:46:12+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-09-22T02:53:52+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"133432238567625447594121178860158261081","date":"2025-09-05T00:09:01+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"94574374111945117627408711398992124875","date":"2025-09-04T17:02:50+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-09-04T16:10:52+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-09-03T02:20:01+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-09-03T02:19:44+00:00","index":"","fulltext":""},{"type":"submitted","content":"Molecular Cancer","date":"2025-09-01T17:46:06+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"molecular-cancer","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"molc","sideBox":"Learn more about [Molecular Cancer](http://gsejournal.biomedcentral.com/)","snPcode":"12943","submissionUrl":"https://submission.nature.com/new-submission/12943/3","title":"Molecular Cancer","twitterHandle":"@SN_Oncology","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"b2b31168-dcd1-4955-a8f5-64e3f4a157d0","owner":[],"postedDate":"September 12th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2026-03-09T16:12:12+00:00","versionOfRecord":{"articleIdentity":"rs-7510857","link":"https://doi.org/10.1186/s12943-026-02637-2","journal":{"identity":"molecular-cancer","isVorOnly":false,"title":"Molecular Cancer"},"publishedOn":"2026-03-07 15:59:58","publishedOnDateReadable":"March 7th, 2026"},"versionCreatedAt":"2025-09-12 14:27:11","video":"","vorDoi":"10.1186/s12943-026-02637-2","vorDoiUrl":"https://doi.org/10.1186/s12943-026-02637-2","workflowStages":[]},"version":"v1","identity":"rs-7510857","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7510857","identity":"rs-7510857","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}