This study uses glioma data from the GLASS consortium. After preprocessing, gene expression TPM values and corresponding clinical outcomes are obtained for 180 IDH wild-type recurrent samples from 158 patients. Approximately 11.1% of the samples are right-censored.

Genes with more than 80% zero expression values and samples with missing data are removed. A five-fold cross-validation is performed to identify markers associated with survival. For each fold, the training and testing sets are normalized separately using log 2 (x + 1) transformation followed by z-score normalization. Due to the high dimensionality of the dataset (over 20,000 genes), a univariate log-rank test is applied to select the top 1,000 survival-related genes for subsequent analysis. Expression levels are dichotomized into high and low groups based on the median values.

As previously described (Liu et al., 2025), we apply one-dimensional hierarchical clustering to group samples for each gene based on their normalized expression values, identifying clusters of patients with the top 10% highest absolute log 2 expression. These clusters form the nodes of a gene-gene interaction network, where edges are defined by shared patients between nodes, enabling the model to focus exclusively on gene-level interactions. To minimize the influence of indirect or complex relationships within the same gene, we randomly select one sample cluster per gene during each graph sampling iteration and repeat this process 1,000 times. We then apply a random walk algorithm to identify paths of connected nodes whose associated patient groups exhibit significantly different survival outcomes, as determined by a log-rank test comparing patients within the path to all others. Finally, we use Fisher’s exact test to identify prognostic markers, testing the null hypothesis that the presence of specific markers is independent of the significant survival-associated paths. The graphical illustration is as Figure 1.

Due to the censored properties of time-to-event outcomes, such as overall median survival time, conventional methods like ANOVA are not appropriate for analyzing survival data in IDH wild-type gliomas. To address this, we employ three regularized Cox proportional hazards (PH) models, specifically Ridge Cox, Lasso Cox, and Elastic Net Cox, for feature selection to identify prognostic biomarkers in IDH wild-type gliomas. These models are well-suited for handling high-dimensional genomic data and account for censoring by modeling the hazard function over time.

The Ridge Cox model applies L2 regularization, which shrinks regression coefficients toward zero to reduce overfitting while retaining all features, making it effective for correlated predictors (Hoerl and Kennard, 1970; Tibshirani, 1997). The objective function for Ridge Cox regression given by Equation 1: L β = ℓ β + λ ∑ j = 1 p β j 2 , (1) where ℓ ( β ) is the Cox partial log-likelihood, β = ( β 1 , … , β p ) ⊤ are the regression coefficients, and λ > 0 is the regularization parameter.

In contrast, the Lasso Cox model uses L1 regularization to perform both shrinkage and feature selection by setting some coefficients to exactly zero, thus identifying a sparse set of predictors (Hoerl and Kennard, 1970; Tibshirani, 1996). The Lasso Cox objective function is as Equation 2: L β = ℓ β + λ ∑ j = 1 p | β j | . (2)

The Elastic Net Cox model combines L1 and L2 regularization, balancing the strengths of Ridge and Lasso to handle correlated features while selecting a subset of relevant biomarkers (Zou and Hastie, 2005; Simon et al., 2011). The Elastic Net Cox objective function is as Equation 3: L β = ℓ β + λ α ∑ j = 1 p | β j | + 1 − α ∑ j = 1 p β j 2 , (3) where α controls the balance between L1 and L2 penalties.

Pathway enrichment analysis was conducted using ToppGene Suite (https://toppgene.cchmc.org/prioritization.jsp), a web-based tool for identifying biological pathways enriched in gene lists (Chen et al., 2009). The analysis integrates multiple pathway databases including Reactome, KEGG, and WikiPathways, and uses hypergeometric testing with Benjamini-Hochberg correction to determine statistically significant enrichments. Pathways with FDR-corrected p-values <0.05 were considered significant.

To evaluate the selected prognostic biomarkers for IDH wild-type gliomas, we employ several specialized survival analysis methods tailored to handle right-censored survival data.

The Cox PH model assumes a linear relationship between predictors and the logarithm of the hazard function, providing a robust framework for assessing the impact of selected biomarkers on survival outcomes (Cox, 1972). The hazard function is modeled as Equation 4: h t ∣ X = h 0 t   exp X β (4)

Where h 0 ( t ) is the baseline hazard function, X represents the covariates, and β are the estimated coefficients. As mentioned, the Elastic Net Cox model extends the Cox PH model by incorporating L1 and L2 regularization.

Random Survival Forests (RSF), an extension of traditional random forests, construct an ensemble of survival trees using bootstrap samples, effectively capturing complex, non-linear relationships between predictors and survival outcomes while accommodating censoring (Ishwaran and Kogalur, 2010; Ishwaran et al., 2008). RSF predicts the cumulative hazard function based on ensemble averaging of Nelson-Aalen estimators from individual trees.

Additionally, Support Vector Machines (SVM) adapted for survival data extend the standard SVM framework by integrating censoring information into the optimization process, modeling the relationship between high-dimensional predictors and time-to-event outcomes (Belle et al., 2011; Pölsterl et al., 2015). The survival SVM formulation focuses on maintaining the correct ranking of survival times among pairs of patients.

To assess model performance, we used two primary metrics:

Concordance Index (C-index): The C-index measures the model’s discriminative ability by calculating the proportion of all possible pairs of patients where the patient with the higher predicted risk experiences the event before the patient with the lower predicted risk (Harrell et al., 1996). Values range from 0.5 (random prediction) to 1.0 (perfect discrimination), with higher values indicating better model performance. The C-index is calculated as Equation 5: C = ∑ i j I y i < y j ⋅ I f ̂ x i > f ̂ x j ⋅ δ i ∑ i j I y i < y j ⋅ δ i (5)

Where y i is the observed time for patient i , f ̂ ( x i ) is the predicted risk, δ i is the event indicator, and I ( ) is the indicator function.

Integrated Brier Score (IBS): The IBS assesses both discrimination and calibration by measuring the average squared difference between observed outcomes and predicted probabilities over a range of time points (Graf et al., 1999; Gerds and Schumacher, 2006). Lower IBS values indicate better model performance. The IBS is defined as Equation 6: IBS = 1 max t i ∫ 0 max t i B S t   d t (6)

where B S ( t ) is the Brier score at time t , calculated using inverse probability of censoring weights to account for censored observations.

Mean C-index and IBS values over 10 replicated runs were used to measure predictive accuracy and calibration, providing a robust assessment of model performance.

PIs were calculated for each patient based on the coefficients of the selected features to quantify individual risk levels. For each survival model, PIs were computed as follows:

Cox PH and Elastic Net Cox: The PI is calculated as the linear predictor from the Cox model as Equation 7, PI cox = X β (7) where X represents the patient’s covariate values and β are the estimated coefficients (Royston and Altman, 2013).

Random Survival Forest: The PI is derived from the ensemble mortality score in Equation 8, PI RSF = 1 B ∑ b = 1 B CHF b t | X (8) where CHF b is the cumulative hazard function estimated from the b th tree, and B is the total number of trees (Ishwaran et al., 2008).

Survival SVM: The PI is based on the decision function value as Equation 9, PI SVM = ∑ i = 1 n α i k X , X i + b (9) where α i are the learned weights, k ( X , X i ) is the kernel function, and b is the bias term (Pölsterl et al., 2015).

Risk stratification was performed by dividing patients into risk groups based on PI scores using two different approaches:

Two-group stratification: Patients were divided into equal-sized low-risk and high-risk groups (0.5, 0.5) based on the median PI value.

Three-group stratification: Patients were divided into low-risk (bottom 1/3), intermediate-risk (middle 1/3), and high-risk (top 1/3) groups based on PI tertiles.

For each stratification approach, Kaplan-Meier survival curves were generated for each risk group, and log-rank tests were performed to evaluate the statistical significance of survival differences between groups (Bland and Altman, 2004). Mean log-rank p-values were computed across the 10-fold cross-validation to summarize model effectiveness in stratifying patients based on survival risk. Additionally, HRs with 95% confidence intervals were calculated to quantify the relative risk between groups (Therneau and Grambsch, 2000).

To evaluate the prognostic impact of selected gene expression markers, we employed an AFT model with a log-normal distribution, following the framework described by Kalbfleisch and Prentice (Kalbfleisch and Prentice, 2002). Unlike the Cox PH model, which focuses on the hazard function, the AFT model directly models the survival time and assumes that covariates act multiplicatively on survival time, effectively accelerating or decelerating the time to event (Klein and Moeschberger, 2003). Specifically, the model expresses the log of survival time T as Equation 10: log T = X β + σ W (10)

Here, X represents the matrix of covariates (i.e., gene expression markers), β denotes the regression coefficients, σ is a scale parameter, and W ∼ N ( 0,1 ) is a standard normal error term. The corresponding hazard function for the log-normal AFT model is as Equation 11: h t = ϕ log t − X β σ σ t 1 − Φ log t − X β σ (11) where ϕ and Φ are the standard normal probability density and cumulative distribution functions, respectively (Collett, 2015). HRs were calculated from the model coefficients using the formula in Equation 12: HR = exp − β σ (12)

This transformation allows us to interpret each marker’s effect on survival time in terms of relative risk (Wei, 1992; Carroll, 2003). Specifically, a HR > 1 indicates an increased risk (shorter survival time), while a HR < 1 indicates a decreased risk (longer survival time).

The log-normal distribution was selected over other parametric options (such as Weibull or log-logistic) based on Akaike Information Criterion (AIC) (Akaike, 1974). We estimated the model parameters using maximum likelihood estimation, with confidence intervals for HRs derived via the delta method (Hosmer et al., 2008).

To ensure reliability, we selected markers that appeared in at least three out of five cross-validation folds using our random-walk-based netSurvival approach, yielding 29 genes. For fair comparison, we also selected the top 29 genes from each regularized feature selection method (Elastic Net, Lasso, and Ridge) based on the sum of ranks across all folds. As shown in the Venn diagram (Supplementary Figure S1), these four methods collectively identified 67 unique genes.

Notably, netSurvival identified 22 unique genes not detected by the regularized methods, indicating that the network-based approach provides a distinct marker profile. This difference likely arises because netSurvival explicitly considers interactions between genes through network topology, whereas regularized Cox models (Elastic Net, Lasso, and Ridge) treat genes as independent features. This fundamental difference in methodology allows netSurvival to uncover biologically relevant markers that may be missed by independence-based approaches.

As shown in Figure 2, netSurvival achieves performance comparable to other methods in both patient risk stratification (measured by integrated Brier score, IBS) and prediction accuracy (measured by C-index). While netSurvival does not demonstrate superior predictive performance, its value lies in (1) identifying unique biologically relevant markers through gene interaction networks, (2) providing enhanced stability with more consistent performance across all five cross-validation folds (indicated by smaller error bars), and (3) offering complementary insights when integrated with traditional regularized methods for comprehensive biomarker discovery.

The Table 1 below presents p-values for group comparisons stratified by PIs. While our method’s p-values are not the lowest across all models compared to other approaches, the selected markers remain promising. Although these markers may not be the optimal predictors of patient risk, their competitive performance warrants further exploration.

We performed the pathway enrichment analysis using all 67 genes. Three significant pathways were identified–FGFR3 signaling in chondrocyte proliferation and terminal differentiation, nanoparticle-mediated activation of receptor signaling, and MYCN to transcriptional activation–each of which shows mechanistic or therapeutic relevance to IDH-wildtype gliomas. The FGFR3 signaling pathway, involving genes such as FGFR3, MAPK1/3, and SOX9, is particularly notable because FGFR3–TACC3 fusions are recurrent oncogenic events in IDH-wildtype glioblastomas, driving proliferative and mitogenic signaling through constitutive FGFR3 activation and MAPK/ERK cascades (Broggi et al., 2022; Gao et al., 2020). The nanoparticle-mediated receptor signaling pathway corresponds to the growing role of receptor-targeted nanocarrier systems–including EGFR-, LRP1-, and TfR-mediated transport–in glioma therapy, enabling blood–brain barrier penetration and tumor-specific drug delivery via receptor interactions (García et al., 2025). Finally, the MYCN transcriptional activation pathway reflects MYCN-driven oncogenic transcriptional programs, which define an aggressive subset of IDH-wildtype gliomas characterized by high proliferation, metabolic reprogramming, and poor prognosis (Borgenvik et al., 2021; Zhou et al., 2023). Together, these pathways underscore the intersection between receptor tyrosine kinase signaling, transcriptional dysregulation, and advanced therapeutic modalities as key axes in IDH-wildtype glioma biology.

We here check the union set of genes from four different feature selection approaches. The features obtained from all methods can be obtained from Supplementary Table S1. To find out markers with reference support, we implement a few-shot prompting through Grok 3 model. All references and findings obtained from the prompting process were rigorously verified through manual literature review to ensure their validity and relevance. The roles of 15 gene candidates with strong or moderate evidence are listed below.

The AFT lognormal model was employed to assess the association between gene expression levels and survival outcomes, with results presented as maximum likelihood parameter estimates. The model identified several genes with statistically significant effects on survival time (p < 0.05), including FN1, PTK2, CCND2, HIF3A, RAD51L3_RFFL, MAX, and EIF4B. The HRs derived from the model provide insight into the impact of these genes on survival. Specifically, FN1 (HR = 1.4034, p = 0.0183) and HIF3A (HR = 1.4907, p = 0.0002) were associated with poorer survival outcomes, indicating that higher expression of these genes accelerates time to event, with a 40.34% and 49.07% increased hazard, respectively. Similarly, EIF4B (HR = 1.5432, p = 0.0013) exhibited a strong adverse effect, with a 54.32% increased hazard for higher expression. Conversely, PTK2 (HR = 0.7601, p = 0.0314), CCND2 (HR = 0.7781, p = 0.0133), RAD51L3_RFFL (HR = 0.7928, p = 0.0164), and MAX (HR = 0.7897, p = 0.0418) were associated with improved survival, with HRs below 1 indicating that higher expression of these genes is linked to a slower progression to the event, reducing the hazard by approximately 23.99%, 22.19%, 20.72%, and 21.03%, respectively. These findings highlight the differential prognostic roles of these genes, with FN1, HIF3A, and EIF4B acting as risk factors, while PTK2, CCND2, RAD51L3_RFFL, and MAX confer a protective effect on survival.

Since recurrence is a key contributor to poor prognosis, we examined whether these 67 markers are also associated with tumor recurrence. As shown in Figure 3, the expression levels of five genes, KIF5C, LICNC00632, B4GALNT3, HIF3A, and RAD51L3-RFFL, differ significantly between primary and recurrent groups in both paired and unpaired Wilcoxon tests. Notably, four of these five genes have strong or moderate evidence supporting their functional roles in glioma, as summarized in Table 2, and two (HIF3A, RAD51L3-RFFL) show significant HRs in the Cox PH model (Table 3).

Pathway enrichment analysis using ToppGene identified three significant pathways providing biological context for these markers: FGFR3 signaling, which drives proliferation through FGFR3–TACC3 fusions and MAPK/ERK activation; nanoparticle-mediated receptor signaling, relevant to emerging therapeutic strategies for blood–brain barrier penetration; and MYCN transcriptional activation, defining an aggressive molecular subtype with poor prognosis.

Our AFT lognormal model revealed that FN1, HIF3A, and EIF4B are significantly associated with poorer survival outcomes, with HRs suggesting substantial increases in risk (40.34%, 49.07%, and 54.32%, respectively). In contrast, PTK2, CCND2, RAD51L3-RFFL, and MAX demonstrated protective effects, with higher expression correlating with improved survival outcomes. These findings align with emerging literature on the biological functions of these genes in glioma pathology.

FN1 (Fibronectin 1) promotes glioma cell invasion and migration by interacting with integrin receptors, particularly α 5 β 1, contributing to its role as a negative prognostic marker (Serres et al., 2014). Similarly, HIF3A (Hypoxia Inducible Factor 3 Alpha) likely contributes to the hypoxic adaptation that characterizes aggressive gliomas, promoting treatment resistance and accelerated disease progression (Monteiro et al., 2017). EIF4B (Eukaryotic Translation Initiation Factor 4B), with the highest hazard ratio among our markers, suggests that dysregulated protein synthesis may be a key driver of aggressive disease behavior in IDH wild-type gliomas.

Conversely, the protective associations observed with PTK2, CCND2, RAD51L3-RFFL, and MAX suggest potential compensatory or tumor-suppressive functions in glioma. PTK2 (Protein Tyrosine Kinase 2) has traditionally been considered pro-oncogenic in many cancers, but recent evidence suggests context-dependent functions that may explain its association with improved survival in our cohort (Sulzmaier et al., 2014). CCND2 (Cyclin D2) may reflect a more differentiated cellular state in glioma, consistent with its favorable prognostic association (Koyama-Nasu et al., 2013).

The identification of five genes (KIF5C, LINC00632, B4GALNT3, HIF3A, and RAD51L3-RFFL) with differential expression between primary and recurrent tumors provides critical insight into the molecular mechanisms driving recurrence, the key challenge in the management of IDH wild-type glioma. Notably, four of these genes have strong or moderate evidence supporting their functional roles in glioma progression. The dual significance of HIF3A and RAD51L3-RFFL, in both survival outcomes and recurrence patterns, underscores their potential as both prognostic biomarkers and therapeutic targets.

KIF5C promotes glioma progression by enhancing cell proliferation, migration, and invasion through the AKT signaling pathway and cytoskeletal organization (Joyce and Lindsay, 2024). In contrast, B4GALNT3 (Beta-1,4-N-Acetyl-Galactosaminyltransferase 3) may influence invasive properties by glycosylating cell surface proteins (Che et al., 2014). The long non-coding RNA LINC00632 represents an emerging class of regulatory molecules with potential roles in gene expression modulation in glioma (Wang et al., 2015). These recurrence-associated biomarkers may provide valuable insights for developing targeted approaches to prevent or delay tumor recurrence, the principal driver of mortality in glioma patients.

Our study benefits from several methodological strengths. First, the focus on recurrent IDH wild-type gliomas addresses a critical clinical challenge, as recurrent tumors represent the most therapeutically resistant phase of disease. Second, the integration of multiple feature selection approaches, including netSurvival, provides robust identification of candidate biomarkers. Third, the application of several survival analysis models enhances the reliability of our findings through methodological triangulation.

However, several limitations must be acknowledged. Despite the GLASS consortium providing one of the largest collections of longitudinal glioma data, our sample size remains modest given the molecular heterogeneity of these tumors. Further validation in independent cohorts is essential. Additionally, while gene expression data provides valuable insights, integration with other molecular data types (e.g., methylation profiles, proteomic data) could yield more comprehensive biomarker signatures. Finally, our study does not address potential treatment-induced alterations in gene expression patterns, which may confound the interpretation of recurrence-associated biomarkers.

The identified biomarkers hold considerable potential for improving clinical management of IDH wild-type gliomas. In the near term, these markers could enhance prognostic stratification, allowing more precise risk assessment and treatment planning. For example, patients with elevated expression of high-risk markers (FN1, HIF3A, EIF4B) might benefit from more aggressive treatment approaches or closer surveillance. In contrast, while those with favorable expression patterns might be candidates for treatment de-escalation strategies to minimize toxicity.

Beyond prognostication, our findings may inform the development of therapeutic interventions. The negative prognostic association of HIF3A suggests that targeting hypoxia adaptation pathways may be particularly valuable in high-risk patients. Similarly, the protective association of PTK2 raises intriguing questions about its context-dependent functions in glioma that warrant further investigation for potential therapeutic exploitation.

Future research should focus on several key areas: (1) functional validation of these biomarkers through in vitro and in vivo models to establish causal relationships with disease progression; (2) integration with spatial transcriptomic approaches to understand the tumor microenvironmental context of these markers; (3) development of clinically applicable assays for biomarker assessment in routine pathology; and (4) investigation of pharmacological approaches targeting the pathways involving these biomarkers.