The overall framework of the study is illustrated in Figure 1, and the specific process is described in the following sections. We obtained a cohort of GBM samples with accessible clinical information, genomics and transcriptomics data from The Cancer Genome Atlas (TCGA) data portal¹ and matched proteomics profile from The Cancer Proteome Atlas (TCPA) repository.² The corresponding H&E histopathological images were obtained from The Cancer Imaging Archive (TCIA).³ A total of 439 GBM patients were selected from TCGA based on the completeness of clinical records and image availability of high-quality histopathological images in TCIA, excluding cases with incomplete data. All included patients had corresponding genomic, transcriptomic, and proteomic data for a comprehensive multi-omics analysis. The GBM tissue microarrays (TMA) of 67 patients with clinical and follow-up data were purchased from Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). Clinical information of patients involved in TMA and TCGA cohorts is provided in Supplementary materials 2, 3.

To extract the quantitative features from whole-slide histopathological images, we applied the Openslide Python library (28) to segment the images into 1,000 × 1,000 pixel sub-images. Furthermore we randomly selected 50 sub-images on behalf of each patient to minimize selection bias and reduce computational load. Image feature extraction was conducted by CellProfiler (29),⁴ an open-source tool for biological-image analysis. The H&E-stained images were converted to grayscale for the extraction of features, which can be specifically categorized into 10 aspects including correlation, image area occupied, image granularity, image intensity, image quality, object intensity, object neighbors, object radial distribution, object size shape and texture. In particular, the textural features were calculated by CellProfiler to quantitatively present the perceived textures of histopathological images, thereby measuring the extent and nature of textures within objects in grayscale images. Through automatic identification and segmentation, these quantitative features objectively interpret the size, shape, spatial distribution, the texture of nucleus and the relationship of pixel intensities, etc. Afterwards, each sub-image was screened to exclude irrelevant features. Eventually, a total of 550 image features were extracted, with the average feature values of 50 representative sub-images of each slide calculated for subsequent analysis.

Two machine learning methods including least absolute shrinkage and selection operator Cox (LASSO-Cox) regression (R package “glmnet”) and support vector machines-recursive feature elimination (SVM-RFE) (R package “e1071”) were performed independently to identify potential informative image features related to prognostic prediction. LASSO-Cox regression applies L1 regularization, effectively reducing multicollinearity, selecting the most survival-associated features and mitigating overfitting by shrinking less relevant coefficients to zero (39). The SVM model can classify data points by maximizing the distance of the hyperplane with high accuracy, thus identifying predictive models or classifiers. SVM-RFE is a feature selection algorithm according to recursive feature deletion sequences with maximum interval principle. It ranks features based on their contribution to classification performance, iteratively eliminating the least informative ones. The integration of LASSO-Cox and SVM-RFE has been demonstrated to improve the model’s generalizability and predictive performance by reducing overfitting and enhancing feature selection reliability (40, 41). Eventually, the features within the intersection of the results by two algorithms were identified as the prognostic histopathological image features (PHIF).

We performed weighted gene co-expression network analysis (WGCNA) based on training set to investigate the association of the prognostic histopathological image features and corresponding gene expression, aiming to further understand the upstream biological mechanisms. WGCNA (42) has been applied to identify modules of genes with highly correlated expression by analyzing the connections between corresponding genes and converting the expression profile into the weighted network. Co-expressed gene networks may facilitate the identification of underlying biological processes, candidate biomarkers and certain clinical traits. Additionally, we applied Metascape for enrichment analysis to estimate the interlinkage between key modules.

In total we included 439 GBM patients with the matched information of histopathological images and other omics from TCGA portal. To minimize overfitting caused by high-dimensional image features, we initially employed XGBoost, GBDT, LASSO, and RF for feature selection and extracted 550 histopathological image features (HIFs) out of the segmented tumor tissue images. Subsequently, to evaluate the clinical practicability of the 550 HIFs, we employed eight algorithms (RF, GBDT, AdaBoost, LR, DT, SVM, NB, and KNN) as classifiers in predicting four common somatic mutations (ATRX, IDH, MGMT, and TERT) and five RNA-based molecular subtypes (classical, mesenchymal, neural, proneural, and G-CIMP). We systematically compared the predictive performances of all classifiers across multiple molecular features, and RF consistently achieved the highest predictive accuracy among the eight classifiers, independent of the feature selection method used. The AUC values for RF models showed superior classification ability across all tested molecular characteristics as shown in Figure 2 and Supplementary material 1. Therefore, we selected RF as a robustly performed algorithm for subsequent prognostic model construction. Additionally, the HIF models validated by GBDT and AdaBoost (ADABAG) also achieved a relatively accurate classification effect under different feature screening methods, which indicates the clinical practicability of HIFs in distinguishing the somatic mutations and molecular subtypes of GBM.

To assess the correlation between histopathological image features (HIFs) and the prognosis of GBM patients, we conducted survival analyses based on individual HIFs. We first assigned the patients into two groups in line with the median value of each HIF (higher than median vs. lower than median) for survival analyses. Afterwards, we carried out univariate Cox analyses based on all HIFs to identify protective prognostic imaging factors, and the top 20 features significantly correlated with the overall survival (OS) was demonstrated in Figure 3A. The four most significant HIFs, with the smallest p-value included one Zernike shape feature (Median_Cells_AreaShape_Zernike_5_5) and three cell texture features (Mean_Cells_Texture_Contrast_3_45, Mean_Cells_Texture_DifferenceEntropy_3_45 and StDev_Cells_Texture_SumAverage_3_0). In particular, Zernike features are a series of 30 shape features based on Zernike polynomials, ranging from order 0 to order 9, which have been frequently extracted for representing the shape parameters in cell nucleus. Cell texture features quantify the correlations between nearby pixels in the regions of interest, which suggests that the global modes of cell nuclei and cytoplasm are all related to clinical survival outcomes. The Kaplan–Meier survival curves of four image features indicated significant differences between groups with high-value and low-value features, demonstrating the feasibility of HIFs in predicting the survival of GBM patients (Figure 3B).

Additionally, according to the expression level of the four predictive features mentioned above, we evaluated the sub-images of high-expressed and low-expressed prognostic features. We utilized TCGA internal validation and TMA external validation cohorts to assess the robustness of the predictive models and reduce the potential overfitting to the specific characteristics of the initial dataset. These validation steps serve as important safeguards against overfitting and bias, which enhances the reliability of our models across diverse datasets. The patients were identified as high-risk and low-risk groups based on the median value of risk scores, and the representative histopathological sub-images showed visible differences in TCGA and TMA external validation cohorts (Figure 3C). The image processing involving cell recognition and segmentation was conducted by CellProfiler, and different cell types were also outlined.

To develop a more accurate predictive model for overall survival (OS) in GBM patients, we estimated the prognostic value of genetic profiles and further incorporated the HIFs with genomics data. Patients were randomly assigned into training (n = 136) and validation (n = 135) sets. To enhance the stability of the measurement, we estimated the mutation status of genes in training set and included the 100 most common somatic mutations in the prognostic model to reduce the dimensionality of the genomics data. The top 15 genes with the most frequent alterations are presented in Figure 4A. Based on the HIFs and 100 mutations we constructed prognosis-relate models in the training set. We applied time-dependent ROC in the validation set since it is more appropriate to represent time-to-event outcomes in the prognostic models compared to the classical ROC curve analysis approach (43). As illustrated in Figures 4C–E, the AUCs for histopathological image features (HIF) model exceeded those of genomics (G) model in 1-year (0.715 vs. 0.634), 3-year (0.813 vs. 0.723) and 5-year (0.829 vs. 0.692) respectively. Moreover, the integrated model of HIF and genomics (HIF + G) reached a better predictive capacity in 3-year and 5-year (AUC = 0.826 and 0.834) than the former two single-omics models. According to the median value of risk score acquired from each model, the patients were then divided into high-risk and low-risk groups. The HIF model and integrative model (HIF + G) showed more accurate prognostic performance (HR = 3.86, 95%CI: 2.67–5.30, p < 0.001, Figure 4) as depicted in Kaplan–Meier curves (Figure 4B).

To further validate the predictive power of the prognostic model, we implemented an external verification using the TMA-GBM cohort. Patients in the external validation set were also divided into high-risk and low-risk groups according to the median risk score. The Kaplan–Meier survival curve revealed a significant difference in survival probability between the groups (p = 0.039, Figure 4F). The 1-year, 3-year and 5-year AUCs of time-dependent ROC were 0.716, 0.712, and 0.703, respectively (Figure 4G). The results thus verified the prognostic capacity of the HIFs in GBM patients.

Transcriptomics can serve as an approach for a comprehensive understanding of the interconnection between the genome, proteome, and cellular phenotype by analyzing the RNA transcripts that reflect the underlying genotype. Based on the training set, we involved 100 whole expressed mRNA genes to decrease the dimensionality and further build the transcriptomics predictive model of OS. The patients were categorized into short-term group (deceased, 12 months ≥ OS ≥1 month) and long-term group (OS ≥60 months) according to the clinical survival status (4, 5, 7, 8). In addition, we applied Metascape for pathways enrichment in the short-term survival group based on the mRNA sequencing data (Figure 5A). Regulation of insulin-like growth factor (IGF) transport and uptake by insulin-like growth factor binding proteins (IGFBPs) has been proven to modulate essential cellular processes and be implicated in certain disorders including malignant, metabolic and immune diseases (44, 45). Previous studies have reported the potential effect of IGF in biological processes associated with tumor growth and invasion inhibition in GBM (46), which may suggest a new effective target for anti-cancer treatment strategies.

As demonstrated in the validation set, the transcriptomics model (RNA) displayed a good predictive performance for OS (1-year AUC = 0.751, 3-year AUC = 0.795 and 5-year AUC = 0.809), which were about equal to the HIF model (1-year AUC = 0.722, 3-year AUC = 0.815 and 5-year AUC = 0.835). Furthermore, we incorporated the transcriptomics and image features as the integrated model (HIF + RNA), which achieved the highest accuracy with the 1-year, 3-year and 5-year AUC increased to 0.769, 0.831 and 0.848 (Figures 5C–E). Additionally, Kaplan–Meier survival analyses also revealed significant differences in survival outcomes between the two groups, with the integrative HIF + RNA model presenting the most notable prognostic value (HR = 7.15, 95%CI: 4.51–10.41, p < 0.001, Figure 5B).

To improve the prognostic prediction of GBM we also incorporated proteomics profile from TCPA portal for further analysis through the reverse phase protein array (RPPA), a high-throughput proteomics method that can assess protein expression and activation states in abundant samples using small amounts of material. In total we involved 179 eligible protein profiles in the proteomics model based on the validation set. The integration of image features and proteomics features (HIF + P) achieved the highest AUCs in 1-year, 3-year and 5-year compared with the proteomics model (0.752 vs. 0.743, 0.835 vs. 0.813, 0.854 vs. 0.818) or the HIF model alone (Figures 6A–C). As shown in the survival analyses, patients in the high-risk group were significantly related to poor OS, and the integrated model (HIF + P) attained the best performance in prognosis prediction among the three models (HR = 6.35, 95%CI: 4.05–9.20, p < 0.001, Figure 6D).

According to the previous analyses, the histopathological image features have presented certain effectiveness in prognostic prediction for GBM patients, and histopathology + omics models have also indicated enhancement in predictive performance and accuracy than the single-omics models. Therefore, we expect to explore the prognostic capacity of a multi-omics predictive model incorporating all the omics features (HIF, genomics, transcriptomics, and proteomics). Based on the validation set, the multi-omics model achieved a 1-year AUC of 0.820, 3-year AUC of 0.926 and 5-year AUC of 0.878, representing an improvement over the HIF + genomics, HIF + transcriptomics and HIF + proteomics models (Figure 7A). Kaplan–Meier survival analysis illustrated a significant difference in survival between high-risk and low-risk groups (HR = 13.14, 95% CI: 7.95–25.95, p < 0.001, Figure 7B). Furthermore, the multi-omics model demonstrated superior net benefit in survival prediction compared to the other models (Figure 7C).

In order to identify the histopathological image features with higher prognostic value for OS, LASSO-Cox regression and SVM-RFE were performed independently. These combined approaches help mitigate the risk of overfitting and ensure the robustness of selected features across different selection frameworks. Previous studies (39–41) have demonstrated that the combination of LASSO and SVM-RFE enhances the reliability of prognostic feature identification in cancer research. A total of five imaging features involved in prognosis were selected via LASSO-Cox regression model, and SVM-RFE selected 12 imaging features with the most significant predictive ability. Ultimately, three overlapped features were identified as prognostic histopathological image features (PHIF), including StDev_Cells_AreaShape_FormFactor, StDev_Cells_AreaShape_Orientation and Mean_Cells_Texture_InfoMeas1_MaskedHematoxylin_3_90 (Figures 8A,B). Representative sub-images and detailed information of patients with high expressed and low expressed PHIF were displayed in Figure 8C and Supplementary material 4.

To explore the upstream genetic mechanisms, we employed WGCNA to construct a gene co-expression network in the training set and identify the gene clusters significantly correlated with the PHIF in GBM samples. Module-trait correlation analysis showed that the red module (219 genes) and turquoise module (868 genes) were significantly associated with the three prognostic image features of GBM among the six identified gene co-expression modules (Figure 9A). Therefore, we defined the red and turquoise module as the key modules of significant prognostic relevance for subsequent research.

Subsequently, we performed an enrichment analysis to explain the biological interpretations of the gene expression profile in the two modules. Genes in the red module were significantly related to several biological processes and pathways such as defense response to other organism, myeloid leukocyte activation, leukocyte cell–cell adhesion, activation of immune response and response to bacterium (Figure 9B). The results indicated that these genes may be involved in immune function, a crucial aspect of tumor immunology, which plays an important role in tumor initiation and progression. The genes in the turquoise module were primarily enriched in categories related to cell morphogenesis involved in differentiation, regulation of neuron differentiation and nervous system development, synapse organization and signaling (Figure 9C). These findings implied that turquoise module genes may have potential association with central nervous system pathways and cerebral function, which may correspond to tumorigenesis and progression in GBM.