In this cohort study, RNA sequencing data were obtained from the TCGA patients diagnosed with glioblastoma. Publically available level 3 transcriptome fragments per kilobase million (FPKM) data of the TCGA-GBM project (n = 166, Illumina HiSeq platform) were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov). In the present study, FPKM values were transformed into transcripts per kilobase million (TPM) values, which are more comparable between different samples (22). Raw counts data of these corresponding GBM samples were also downloaded for differential expression analysis. The GISTIC copy number and methylation status of splicing factors in the TCGA-GBM cohort were downloaded from cBioPortal for Cancer Genomics (https://www.cbioportal.org/). Curated clinical information of these glioblastoma patients was derived from a TCGA Pan-Cancer research (23).

AS events profile of 152 glioblastoma patients in the selected TCGA-GBM cohort, a matrix consists of the Percent Spliced In (PSI) value of each sample, was obtained from the TCGA SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/PSIdownload.jsp). PSI value means the ratio between exons inclusion and exclusion read counts and indicates the efficiency of certain splicing process (24). To ensure the reliability of subsequent analyses, we predefined stringent filtering criteria: the PSI values of AS events in all GBM samples with a standard deviation > 0.05, mean > 0.1, and the percentage of samples with PSI values > 80 were included.

Univariate Cox regression analysis was carried out to find out the prognostic relationship between AS events and overall survival of GBM patients. AS events with P-value less than 0.05 in the univariate Cox regression analysis were regarded as survival-related factors and were selected for the next analysis. The least absolute shrinkage and selection operator (LASSO) cox model using 10-fold cross-validation was performed on the top 20 significant survival-related AS events that have the highest prognostic value to select the most useful prognostic markers of all types of AS events. The optimal lambda value was estimated based on the minimum criteria. Then, these candidate survival-related AS events were further subjected to multivariate Cox regression analysis (25). The risk-score model of AS events for predicting survival of GBM was constructed based on the combination of an AS PSI value and the corresponding regression coefficient (β) obtained from the multivariate Cox regression analysis (26). The formula of calculating the risk score of each GBM patient was as follows:]

where N is the number of AS events in the signature, β_ASi is the regression coefficient, and PSI_ASi represents the PSI value of a certain AS event in a single sample. GBM patients were divided into a high- and low-risk group according to the median risk score. Kaplan-Meier survival analysis based on the log-rank test was used to evaluate the predictive value of the risk model. Hazard ratios (HRs) and 95% confidence intervals (95% CIs) were calculated both in univariate and multivariate Cox regression analyses. The following clinical-relevant factors were included in the multivariate Cox regression analysis: gender, age, MGMT methylation status, Karnofsky Performance Status (KPS) score, IDH mutation status, and risk score. Finally, time-dependent receiver operating characteristic (ROC) curves were performed using the “timeROC” R package (v0.4) to display the accuracy of this risk model in predicting the clinical outcomes of GBM patients at different times.

Based on all survival-related AS events, a modified clustering method, SNF-CC, which combining Similar network fusion (SNF) (27) with Consensus clustering (CC) (28), was applied to perform classification of all GBM patients in the cohort. SNF-CC algorithm was executed by “ExecuteSNF.CC” function implanted in “CancerSubtypes” R package (v1.12.1) (29). To guarantee a balance between high stability and low ambiguity, the detailed parameters were set as follows: clusterNum = 2, K = 20, alpha = 0.5, t = 20, maxK = 5, pItem = 0.8 and reps = 500. The consensus heatmap and cumulative distribution function (CDF) were used to select a more appropriate number of clusters. Furthermore, silhouette width, an index represents how similar a sample is matched to its identified cluster compared to other clusters, was validated in our study. The associations between AS events-based subtypes, clinicopathological features (age, overall survival status, KPS, IDH mutation, and MGMT methylation status), GBM molecular subtypes (Verhaak classification and EM/PM classification), and immune features (described below) were evaluated.

Gene set enrichment analysis (GSEA) analysis was conducted using the R package “fgsea” (v1.12.0) (30). We ranked the GBM samples according to their log2-fold change value (from high to low) derived from differential expression analysis between two distinct groups and checked if any signaling pathway or molecular hallmark was associated with these differentially expressed genes. GSEA gene sets (curated gene sets (C2), Gene Ontology (GO) gene sets (C5), and hallmark gene sets (H)) were downloaded from the Molecular Signatures Database (v7.0) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). To obtain robust results, gene-set permutations were performed 10,000 times, and enrichment P-values were adjusted by false discovery rates (FDR). FDR-adjusted P < 0.05 were considered as statistically significant. GO term enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for genes of interest were performed using metascape (v3.5), which is an easy-to-use web portal that provides a comprehensive analysis for the functional annotation of lists of genes (31).

The fractions of diverse infiltrated immune cells in tumor samples were estimated using the ssGSEA method implemented in the “GSVA” R package (v1.34.0, https://bioconductor.org/packages/release/bioc/html/GSVA.html). As an extension of Gene Set Enrichment Analysis (GSEA), ssGSEA evaluate the enrichment score of a certain gene set in every single sample (32). Twenty-four kinds of immune cell marker genes derived from a previously published research were integrated into immune cells specific gene sets (33). Markers associated with cells of the innate immune system, including natural killer (NK) cells, NK CD56dim cells, NK CD56bright cells, dendritic cells (DCs), immature DCs (iDCs), activated DCs (aDCs), neutrophils, mast cells, eosinophils, and macrophages, as well as those associated with cells of the adaptive immune system, including B, T central memory (Tcm), CD8+ T, T effector memory (Tem), T follicular helper (Tfh), Tγδ, Th1, Th2, Th17, and Treg cells, were included in the gene sets list. Finally, ssGSEA captured a numeric matrix containing enrichment scores of different immune cells across all GBM samples.

Immune infiltration features in GBM samples can be inferred by calculating tumor immunological indexes. According to a pan-cancer study aiming at developing an effective biomarker for predicting immunotherapy response (34), immune infiltration score (IIS), T cell infiltration score (TIS), antigen presentation machinery score (APS), and tumor immunogenicity score (TIGS) were computed. Briefly, IIS indicated total immune infiltration level in the tumor sample and was calculated as the mean of standardized infiltration scores for all kinds of immune cells obtained from the GSVA algorithm. Similarly, TIS was calculated using T cell subsets, including CD8+ T, T helper, T, Tcm, Tem, Th1, Th2, Th17, and Treg cells. APS was calculated with GSVA using APM related genes (PSMB5, PSMB6, PSMB7, PSMB8, PSMB9, PSMB10, TAP1, TAP2, ERAP1, ERAP2, CANX, CALR, PDIA3, TAPBP, B2M, HLA-A, HLA-B, and HLA-C). Tumor burden (TMB) was the number of non-synonymous mutations per megabase (MB). Here, we used 38 MB as the estimate of exome size and defined the TMB as the quotient of non-synonymous mutations and 38MB. Then, TIGS was calculated by multiply APS and TMB value. Predicted abundance of neoantigens in GBM samples was obtained from a previous study (35).

The subclass mapping method (SubMap) (36) was used to compare the gene expression matrices of different GBM subtypes with the expression profiles of several cancer types treated with ICB therapies, including transcriptomic data of 47 melanoma patients who received immunotherapy targeting CTLA-4 and PD-1, and NanoString data of 49 patients with four cancer types treated with anti-PD1 therapy (37, 38). This step was performed using the SubMap module (v3) on the GenePattern website (http://genepattern.broadinstitute.org/). Parameters of SubMap analysis were set as default: num marker genes = 100, num perm = 100, and num perm fisher = 1,000. P-values were adjusted by Bonferroni correction.

A list of splicing factors (SFs) was obtained from the SpliceAid2 database (http://www.introni.it/splicing.html). SpliceAid2 database embodies experimental curated splicing factors to help researchers understand the tissue-specific alternative splicing. The SF gene expression profile of GBM patients was retrieved from TCGA RNA-seq data. Overall survival-related SFs were determined using univariate Cox regression analysis. The correlation of SFs and AS events was evaluated by Pearson correlation analysis, and the regulatory network plot was generated in Cytoscape software (v3.8.0).

Multivariable Cox proportional hazards regression analysis was applied with the following clinical-relevant covariates: gender, age, MGMT methylation status, Karnofsky Performance Status (KPS) score, IDH mutation status, and risk score. A combined nomogram was generated by the “regplot” R package (v1.0) as a quantitative tool for predicting the likelihood to die of each patient. The concordance index (C-index) was calculated to assess the consistency between model prediction and actual clinical outcomes of patients. The calibration plot was generated to evaluate the accuracy of the prediction for 1- and 2- year overall survival using this nomogram by the “rms” R package (v5.1-4). Additionally, decision curve analysis (DCA) was applied to evaluate the performance of the nomogram by the “rmda” R package (v1.6).

The intersections of seven types of AS events in GBM were plotted using the “UpSetR” R package (v1.4.0) (39). The “glmnet” package (v3.0-2) was used to conduct the LASSO Cox regression model analysis. Kaplan-Meier survival curves were visualized to discover the difference of clinical outcomes between groups using “survival” and “survminer” R package (v0.4.6). 3D-PCA plot was generated using “mixOmics” (6.10.9) (40) and “rgl” (v0.100.5) (https://cran.r-project.org/web/packages/rgl/index.html) R package.

Restricted mean survival (RMS) represents the loss in average life expectancy for patients. We performed RMS time calculation based on the univariable Cox proportional hazards regression of overall survival (41). RMS and time ratio were estimated using the “survRM2” R package (v1.0-2).

ESTIMATE is a widely used method that can infer the fraction of immune and stromal cells in tumor samples using the gene expression profile (42). Tumor purity was defined as the percentage of malignant cells in a solid tumor sample. Here, the purity of each GBM sample was assessed by package “ESTIMATE” (v1.0.13, https://bioinformatics.mdanderson.org/estimate/rpackage.html). Furthermore, the tumor purity was further validated using the “TPES” R package (v1.0.0), which is a computational method for estimating tumor purity from single-nucleotide variants. Here, the tumor purity indexes of the TCGA-GBM samples calculated using the “TPES” algorithm were derived from the supplementary files of published literature (43).

All statistical analyses were performed using R (version 3.6.1, www.r-project.org), with chi-square or Fisher’s exact test for patients’ characteristics in different clusters, a Mann-Whitney U test for testing the differences in means of continuous data, and log-rank test for survival analysis. For all hypothetical tests, a two-sided P-value < 0.05 was considered to be statistically significant.

Integrated mRNA AS events profile was analyzed for 152 GBM patients from the TCGA cohort. Under the stringent filtering criteria, a total of 15,907 AS events from 6,491 genes were identified ( Table S1 ). These AS events consisted of seven different splicing patterns: Alternate Acceptor site (AA), Alternate Donor site (AD), Alternate Promoter (AP), Alternate Terminator (AT), Exon Skip (ES), Retained Intron (RI), and Mutually Exclusive Exons (ME). Visualization of specific numbers and intersections of all types of AS events in GBM was generated in Figure 1A . Among all these types, ES and AP splicing modes accounted for almost 60% of the total number (9,478/15,907). Notably, a single gene may have more than one type of splicing patterns. We detected 2,184 genes contain two or more types of AS events in GBM samples. For instance, PDE4DIP has six types of AS events, including AA, AD, AP, AT, ES, and RI.

To explore the correlations between AS events and the overall survival of GBM patients, we performed univariate Cox regression analysis using PSI values of all AS events and survival information of patients. Eventually, 775 AS events of 593 genes were confirmed to be prognostic in GBM (P-value < 0.05). Among them, we noticed that some genes have more than one type of prognostic AS events. For example, AA and ES splicing modes of ATG4D were all associated with overall survival. The UpSet plot diagram was illustrated to visualize the interactions among these AS events ( Figure 1B ). A more intuitive landscape of AS events and their prognostic value in GBM was shown using the circos plot ( Figure 1C ).

Numerous survival-related AS events were identified after performing univariate Cox regression analysis. In order to get the AS events with the greatest potential prognostic value, LASSO Cox regression analysis was carried out, and 18 AS events were finally selected ( Figures 2A, B ). We also validated whether the original genes of these 18 AS events are survival-related factors when applying the univariate Cox regression analysis to the corresponding gene expression profiles, and we found that only two of them remain statistically significant, including SERGEF and FAM86B1 (P-value < 0.05) ( Figure 2C ). Besides, only SERGEF and its splicing isoform, SERGEF-14562-AD, shared similar prognostic value, while FAM86B1 and FAM86B1-82719-AD had the reverse patterns in prognostic value. As revealed above, we identified and characterized that isoform-based analysis can capture some meaningful transcripts that were not available for gene-level analysis.

Then, we developed a formula to compute the risk score for each GBM patient based on the PSI values and regression coefficients of AS events. The risk score of each GBM sample ranged from -13.782 to -7.970, and the median score of -10.569 was defined as the cut-off to separate the whole GBM cohort into high- and low-risk groups. The efficacy of the risk model for predicting overall survival of GBM patients was evaluated by performing Kaplan–Meier survival analysis. The result indicated that the risk score-based signature was significantly associated with the prognosis of GBM patients. Patients in the high-risk group had worse clinical outcomes than those in the low-risk group (HR = 3.34, 95% CI: 2.23–5, P < 0.001) ( Figure 2D ). The performance of the risk model was assessed in terms of the RMS time ratio between different risk groups, and significant shorter RMS time were observed in high-risk group compared with their low-risk counterparts (RMS_low-risk/RMS_high-risk = 1.793, 95% CI: 1.514–2.123, P < 0.001) ( Table S2 ). Time-dependent ROC curves evaluated the performance of this risk model in predicting outcomes, and the area under the ROC curve was 0.81 at 1-year, 0.92 at 2-year, and 0.91 at 3-year, respectively ( Figure 2E ).

The risk score distribution, survival status, and AS events profile of this signature were shown in Figure 3A . More patients were found dead in the high-risk group than the low-risk group, and the overall survival time of high-risk patients was much shorter than patients of low-risk. We further calculated the different usage of the transcript isoforms in the alternative splicing signature in the different risk groups. The PSI value represents the frequency of the specific transcript isoform that occurred in the alternative splicing process. Therefore, we compared the PSI values of these AS events in the risk model signature to reflect the differential usage of the mRNA isoforms. The abundance of the risky AS events with HR greater than one in the high-risk group were more than the low-risk group, while the opposite pattern occurred in the protective AS events ( Figure 3B ). These results indicated that the substantial differences in survival between high- and low-risk groups might result from the isoform switching.

Verhaak’s GBM molecular classification was developed based on the gene expression profile of the TCGA GBM cohort and has been widely used for GBM research (9). The Verhaak subtype information of each sample in the present study was obtained from UCSC Xena (https://xenabrowser.net/). Next, we investigated the inter-tumor heterogeneity of risk score by examining the relationship of risk score with different Verhaak’s subtypes. Higher scores were found in the mesenchymal subtype compared to the classical and proneural subtype ( Figure 3C ). IDH mutation and MGMT methylation status are well-documented GBM molecular markers that can predict the overall survival of patients (44). In general, IDH wild-type and MGMT unmethylated status usually indicate a worse prognosis. We separated GBM patients into different groups based on these two biomarkers to explore the association between risk score and IDH mutation/MGMT methylation status ( Figure 3D ). We noticed that the risk score was significantly higher in the IDH wild-type group than the IDH mutation group (P < 0.001). A similar result was also found in MGMT unmethylated group (P < 0.01). Furthermore, the HRs for the AS events-based signature in the univariate and multivariate Cox regression analyses were 2.718 (P < 0.001, 95% CI: 2.186–3.38) and 2.496 (P < 0.001, 95% CI: 1.802–3.458), respectively ( Table S3 ). Thus, the risk score model was proved to be an independent prognostic factor for GBM patients.

Patients without detailed survival information were excluded before performing the SNF-CC algorithm, and then PSI values of AS events in all samples were scaled via z-score standardization. The built-in Cox model function of the “CancerSubtypes” package enabled us to filter crucial survival-related splicing patterns specific to GBM. The performance of this clustering method was assessed by clustering heatmap, CDF curves, and silhouette width. The results demonstrated that SNF-CC achieved adequate robustness when all GBM patients were categorized into two distinct clusters, with Cluster1 consists of 45 patients and Cluster2 of 106 patients ( Figure 4A , Figure S1 ). PCA analysis supported the effectiveness of this clustering method in defining GBM sub-grouping based on survival-related AS events ( Figure 4B ). Importantly, statistically significant differences in outcomes of patients were observed in these two clusters, and Cluster2 patients showed worse prognoses than those of Cluster 1 (HR = 2.07, 95% CI: 1.44–2.97, P < 0.001) ( Figure 4C ). A higher proportion of patients with IDH wild-type (fractions, 99.03% vs 81.82%, P = 3.06e-4) and MGMT unmethylation (fractions, 64.37% vs 44.44%, P = 0.047) was observed in Cluster2 compared with Cluster1 ( Figure 4D ). Also, we found a significant difference in Verhaak’s subtype among two clusters (P = 4.97e-11), and Cluster2 had more mesenchymal GBMs than Cluster1 (fractions, 55.24% vs. 15.91%) ( Figure 4E ). The EM/PM subtype is another glioma molecular classification based on the coexpression modules of EGFR and PDGFRA (45). The prognosis of EM glioma was much worse than PM glioma. Here, we noticed a significant relationship between EM/PM subtype and the novel classification scheme we developed (P = 1.23e-5), and Cluster2 had more EM^high samples than Cluster1 (fractions, 71.43% vs 50%).

The remarkable differences in both clinicopathological and molecular characteristics between these two clusters suggested that different functional annotations and signaling pathways might exist. GSEA manifested a remarkable activation in cancer-associated signaling pathways was significantly enriched in Cluster2 compared with Cluster1, including epithelial-mesenchymal transition (NES = 2.27, FDR = 5.1e-4), P53 (NES = 1.52, FDR = 4.3e-3), IL6 JAK STAT3 signaling (NES = 2.35, FDR = 5.1e-4), and interferon-gamma response pathways (NES = 2.5, FDR = 5.1e-4) ( Figure 5A , Table S4 ). Besides the cancer hallmark biological processes, Cluster2 was also involved in various immune-related responses, such as inflammatory response (NES = 2.55, FDR = 2e-3), adaptive innate immune response (NES = 2.59, FDR = 2e-3), innate immune response (NES = 2.49, FDR = 4.3e-3), and immune effector process regulation (NES = 2.34, FDR = 2e-3). Overexpressed genes in Cluster2 (log2-fold change > 1 and adjusted P-value < 0.05) were calculated as input for functional annotation analyses using Metascape. The GO and KEGG enrichment analysis indicated that Cluster2 had enriched genomic biological processes, including human immune response, cytokine-cytokine receptor interaction, T cell activation, cytokine production, IL-10 signaling, and regulation of leukocyte migration ( Figure 5B ).

The intratumoral microenvironment is a complex that consists of the tumor and non-tumor cells, such as stromal and immune cells (46). These non-malignant cells can regulate the progression of tumorigenesis via the cross-talk with malignant cells in GBM (47). We found that both IIS and TIS were much higher in Cluster2 than in Cluster1, which represented higher relative fractions of total immune cells (P < 0.001), and T cell subsets (P < 0.01) were infiltrated in tumors of Cluster2 ( Figure 5C ). Recently, immunotherapies targeting immune checkpoints blockade (ICB) had been proved to exhibit critical anti-tumor functions by promoting anti-tumor immune responses and inhibiting immunosuppressive effects. Clinical outcome of ICB therapy has been reported to be closely related to neoantigen abundance in some cancers (48, 49). Here, a significant depletion of the neoantigen burden was found in Cluster2 (P < 0.05). Besides, we also noticed a relatively higher level of APS in Cluster2 compared with Cluster1 (P < 0.001). In glioma, decreasing of neoantigens was tightly associated with intact APM function, increased infiltration of immune cells, and active immune processes (50). This may explain the phenomenon that Cluster2 contained higher APS, IIS, and TIS while harboring lower numbers of neoantigens.

To further validate the association between immune system processes and these two subtypes, the ssGSEA method was performed to estimate the differences in the detailed immune infiltration of 24 types of immune cells between these clusters ( Figure 6A ). We found that Cluster2 exhibited a higher abundance of pro-tumor immune cell types, including mast cells, immature dendritic cells (iDCs), CD56^dim natural killer cells, macrophages, neutrophils, and regulatory T cells (Tregs) ( Figure 6B ). Meanwhile, anti-tumor immune cells, such as CD8+ T cells and B cells, were significantly significantly depleted in Cluster2 samples. These results indicated that the enrichment of numerous immune-related pathways in Cluster2 might result from the differences in recruitment or differentiation of diverse immune cells in tumors.

Wang et al. pointed out that TIGS is an ideal index to predict clinical response to ICB therapy (33). Therefore, a higher TIGS level in Cluster2 suggested more patients may benefit from ICB therapy than those of Cluster1. The expression levels of immune checkpoints have been proposed to be useful references for selecting patients receiving immunotherapy. Thus, we further investigated the expression levels of several well-known immune checkpoints (TIM-3, TIGIT, LAG-3, PD-1, CTLA-4, PD-L1, PD-L2, and IDO1) between two distinct clusters. Most of them were more highly expressed in Cluster2 than in Cluster1 (P < 0.05) ( Figure 6C ). Finally, SubMap analysis manifested that samples in Cluster2 shared a higher similarity with the expression profile of melanoma patients who were responsive to PD-1 inhibitor treatment (P = 0.004, Bonferroni P = 0.032) ( Figure 6D ). Moreover, the same procedure performed on another cancer cohort containing 49 baseline tumors in four cancer types, GSE93157, also achieved similar results (P = 0.03, Bonferroni P = 0.12) ( Figure 6D ). These findings confirmed that GBM samples in Cluster2 may benefit from anti-PD-1 therapy compared with those of Cluster1.

We investigated the association between AS events-based risk score and tumor purity in GBM using both ESTIMATE and TPES algorithms, and we found a stable significant negative correlation among them (ESTIMATE: Pearson correlation coefficient (R) = -0.2, P = 0.015; TPES: R = -0.24, P = 0.0098) ( Figure 7A ). Meanwhile, higher stromal and immune score were observed in high- than the low-risk group (P < 0.05) ( Figure S2 ). We also explored the relationship between risk score and glioma mediated immunosuppression. The gene sets of glioma-related immunosuppressive factors, including immunosuppressive cytokines and checkpoints, tumor-supportive macrophage chemotactic and skewing molecules, immunosuppressive signaling pathways, and immunosuppressors were extracted from previous reports (51, 52). The risk score was positively correlated with the expression of most genes ( Figure S3 ). Pearson correlation analysis also found that risk score was positively correlated with immunosuppressive genes including TIMP1 (R = 0.450), BIRC3 (R = 0.435), ICAM1 (R = 0.361), CCL2 (R = 0.359), and RAB27A (R = 0.352) ( Figure 7B ). After that, ssGSEA analysis was carried out to assess the enrichment score of each gene set for every single GBM sample. Positive correlations between risk score and each gene set score were found ( Figure 7C ), indicating that AS events-based signature played a vital role in immunosuppression in GBM microenvironment.

SFs play an important role in regulating the process of exons inclusion and introns exclusion during alternatively splicing pre-mRNA. Changes of SFs lead to the production of diverse splicing patterns of genes, including some oncogenic isoforms, and thus promote or inhibit tumorigenesis (53). In order to depict the potential regulatory network between SFs and prognostic AS events, univariate Cox regression analysis was carried out to identify the survival-related SFs in GBM. 16 SFs were regarded as core SFs (P < 0.05), including CELF6, HSPB1, HSPA5, TIA1, PQBP1, TIAL1, EEF1A1, FAM50B, NUDT21, SMNDC1, HNRNPC, PRMT5, JUP, MYEF2, LSM2, and HNRNPLL. Significant correlations with |R| > 0.4 and P < 0.05 were shown in the network map ( Figure S4 , Table S5 ). A total of 14 prognostic SFs were correlated with 116 AS events. Several most significant correlations were presented in Figure 8A .

Posttranscriptional modification of SFs can also influence alternative splicing, such as phosphorylation and methylation (54, 55). In this study, the methylation levels of EEF1A1, FAM50B, PRMT5, and SMNDC1 promoters were negatively associated with their mRNA expression levels ( Figure 8B ). Next, we investigated the relationships between copy number alteration (CNA) and expression levels of prognostic SFs, and the expression levels of 10/14 SFs were associated with CNA events ( Figure 8C ).

To develop a quantitative tool for predicting the prognosis of GBM patients, we established a nomogram by integrating clinicopathological risk factors and AS events-based signature based on the multivariable Cox proportional hazards model ( Figure 9A ). The point scale in the nomogram was utilized to generate point to these variables, and the risk of death of each GBM patient was qualified by accumulating total points of all variables. The risk score was found to have the most excellent weight among all these variables, which was consistent with the result of the previous multivariable Cox regression analysis. The C-index of this nomogram reached 0.774 (95%CI: 0.743–0.805). The result of the calibration plot further confirmed the significant consistency between predicted and observed actual clinical outcomes of GBM patients ( Figure 9B ). The decision curve analysis showed a higher overall net benefit using the nomogram than either the “treat all” or the “treat none” approach within a range of threshold probabilities > 10% ( Figure 9C ). Moreover, the same result was found when compared with the base model, which contained age, gender, KPS, IDH status, and MGMT methylation status of patients. All these findings demonstrated that the AS events-based nomogram is an optimal model for predicting the prognosis of GBM patients in clinical management.