Figure 1 displays the process of this research. The transcriptome expression profile (TPM) and associated clinical data for GBM were downloaded from public databases, including TCGA and GEO. We obtained a TCGA GBM cohort and two GEO GBM cohorts (GSE83300 cohort and GSE74187 cohort) for subsequent analyses. First, the original files were background adjusted and quantile normalized, and then the two datasets were combined after eliminating batch effects by applying the “Combat” algorithm. Excluding samples without OS information, 209 GBM samples were included. Table S1 provides detailed clinical information about the 209 GBM patients, including age, sex, overall survival time, survival status, and IDH1 mutation status.

Twelve CRGs were retrieved from previous publications. In consensus clustering analysis, the R package “ConsensusClusterPlus” was applied. The clustering criteria were as follows: the cumulative distribution function (CDF) curve is relatively flat rather than steep, the sample size of each subtype should not be too small, and the correlation within subtypes was enhanced after clustering, while the correlation between subtypes was weakened. Based on the expression profile of prognostic CRGs, samples were classified into distinct cuproptosis subtypes.

To identify the clinical significance of the classification above, analyses of the relationship between cuproptosis subtypes and age, sex, and prognosis were conducted, respectively. Using R packages “survminer” and “survival”, we conducted a Kaplan-Meier survival analysis to compare the difference of OS in different cuproptosis subtypes.

Based on the ssGSEA algorithm, we obtained the relative contents of 23 kinds of immune cells in every GBM sample. Furthermore, in light of the ESTIMATE algorithm, we obtained the ImmuneScore, StromalScore, and ESTIMATEScore of every GBM sample. Three of the most studied immune checkpoints were selected, including CTLA-4, PD-L1, and PD-1, and differentially expressed analyses of checkpoints between two cuproptosis subtypes were performed.

Based on the R package “limma”, we obtained DEGs between the two cuproptosis subtypes with screening conditions of fold change (FC) > 1.5 and adjusted P< 0.05. Then, a univariate Cox regression analysis was conducted to obtain the prognostic DEGs. Finally, in light of the prognostic DEGs, unsupervised cluster analysis was performed to categorize GBM samples into distinct gene subtypes.

GBM samples were randomly categorized into the training set and the testing set. Sample sizes of the two sets were about the same, with 105 GBM patients in the training set and 104 GBM patients in the testing set. We chose the training set to establish a risk score model. Briefly, using the R package “glmnet”, Lasso Cox regression analysis was conducted in light of the prognostic DEGs. We established the risk score model by analyzing the change trajectory of every gene and using the 10-fold cross-validation method.

The calculation formula is as follows: Risk score = Σ (Coefi × Exp) In this formula, Coefi indicates the risk coefficient, and Exp indicates the expression level. In the training set, samples were categorized into high-risk and low-risk clusters in light of the median. Then a Kaplan-Meier survival analysis was conducted to compare the OS difference between the high-risk and low-risk clusters. Next, using the R package “ggplot2”, PCA was performed to observe the distinction between the two risk clusters. Finally, the ROC curve was constructed to assess the prognostic prediction performance of this model. For further verification, samples in the testing set, GSE83300 cohort, and GSE74187 cohort were categorized based on the median risk score, respectively, followed by Kaplan-Meier survival analysis, PCA, and ROC curve construction.

The Ethics Committee permitted this experiment, and all included patients signed informed consent. Five pairs of GBM specimens and nearby non-tumor specimens were collected from the Department of Neurosurgery from 2017 to 2020. Samples were stored at -80°C immediately after resection and then transferred to liquid nitrogen.

Firstly, based on the manufacturer’s instructions, tissue RNA was extracted by Trizol (Invitrogen, USA). Secondly, reverse transcription was conducted by PrimeScript™ RT Mix (Takara) containing random and oligo primers to obtain cDNA. Finally, qPCR analysis was performed on a Light-Cycler 480 qPCR instrument with LightCycler^® 480 SYBR Master (Roche). Here are the PCR conditions: preincubation at 96°C for 6 minutes, 39 cycles of amplification with 9 s at 93°C, 19 s at 58°C, followed by an extension at 68°C for 16 s. The expression levels of target genes were obtained by a formula of 2^-ΔΔCt. Table S2 provides sequences of primers used in this research.

The protein expression levels of model-related genes were verified by Western blotting. Samples were treated with RIPA lysis buffer. A BCA assay kit (Thermofisher, USA) was used to determine the protein concentration. An equal amount of protein (20 μg) was loaded into lanes. After being separated by electrophoresis, proteins were electrically transferred to a PVDF membrane (Millipore, USA). After blockaded with 5% milk, the membrane was incubated with the primary antibodies: anti-PDIA4 (SAB1404743, Sigma-Aldrich, St. Louis, MO, USA), anti-DUSP6 (SAB1410312, Sigma-Aldrich, St. Louis, MO, USA), anti-PTPRN (ab207750, Abcam, Cambridge, MA, USA), anti-PILRB (ab198267, Abcam, Cambridge, MA, USA), anti-CBLN1 (ab181379, Abcam, Cambridge, MA, USA), and anti-GAPDH (ab181602, Abcam, Cambridge, MA, USA), each was diluted at a ratio of 1:1000 and incubated over-night at 4 °C. After incubation with the corresponding horseradish peroxidase-linked secondary antibody at room temperature for 1 h, target proteins were developed by the enhanced chemiluminescence kit (Millipore, USA).

Three clinical characteristics, age, gender, and IDH1 mutation status, were selected, and samples were categorized based on each clinical characteristic. Then we conducted a differential analysis of risk score with the Chi-square test. We conducted univariate and multivariate regression analyses to verify if the risk score is an independent prognostic factor. Moreover, we conducted stratification analyses. Briefly, the samples were stratified according to each selected clinical characteristic. The stratified samples were further stratified in light of the median risk score, followed by Kaplan-Meier survival analyses to observe if the OS difference in each subgroup was statistically significant.

In light of the CIBERSORT algorithm, we obtained the relative contents of 22 types of immune cells in every GBM sample. Next, we conducted correlation analyses between risk scores and the contents of immune cells and the correlation analyses between the contents of immune cells and the expression of model-related genes. Moreover, in light of the ESTIMATE algorithm, we obtained the ImmuneScore, StromalScore, and ESTIMATEScore of every GBM sample. The differential analyses of ImmuneScore, StromalScore, and ESTIMATEScore between the high-risk and low-risk clusters were conducted. Moreover, we picked multiple immune checkpoint molecules followed by differential analyses of immune checkpoint expression levels between the two risk clusters. Finally, a correlation analysis between CSC indexes and risk scores was conducted. The CSC index is an indicator describing the degree of similarity between tumor cells and stem cells, which can be considered as the quantification of CSCs. Each CSC index ranges from low (zero) to high (one) stemness. The CSC index was calculated using the innovative one-class logistic regression (OCLR) machine learning algorithm (11).

Based on the R package “maftools”, we obtained the TMB of every GBM sample. Next, the differential analysis of TMB between the two risk clusters was conducted. Moreover, we calculated the IC50 values of chemotherapeutic drugs commonly used using the R package “pRRophetic”. It’s an R package predicting chemotherapeutic response based on tumor gene expression profiles (12). Next, differential analyses of the IC50 values between the high-risk and low-risk clusters were performed to assess the differences in therapeutic effects.

Using the R package “rms”, a predictive nomogram was constructed according to clinical features and risk score. A patient’s score equals the sum of the corresponding scores for every variable. Based on the total score, we can predict the survival rates of a patient at 0.5, 1.0, and 1.5 years. Moreover, ROC curves were constructed to obtain the AUC values of 0.5, 1.0, and 1.5 years to assess the performance of this nomogram. Finally, a calibration plot was established to observe the differences between the predicted survival events of 0.5, 1.0, and 1.5 years and the actual results, thus further assessing the performance of this nomogram.

R was used for most of the statistical analyses in this study.

In this study, a total of 12 CRGs were included. We conducted differentially expressed analyses of CRGs between normal and GBM samples and found that 10 CRGs were significantly up-regulated in GBM samples, including SLC31A1, CDKN2A, MTF1, LIPT1, FDX1, PDHB, PDHA1, LIAS, DLD, and DLAT. Conversely, 2 CRGs were significantly down-regulated, including ATP7B and GLS (Figure 2A).

Next, analysis of somatic mutations in 12 CRGs revealed a low-frequency mutation in GBM samples. Of 461 GBM samples, there were 16 (3.47%) mutations in CRGs. Five CRGs (CDKN2A, MTF1, ATP7B, DLD, and GLS) had a mutation frequency of 1%, respectively, while three CRGs (FDX1, LIAS, and LIPT1) did not have any mutations (Figure 2B). Furthermore, we analyzed the copy number variation (CNV) in 12 CRGs and discovered frequent CNV in only 2 CRGs. CDKN2A and ATP7B showed an apparent CNV decrease (Figure 2D). Figure 2C displays the locations of 12 CNV on their respective chromosomes. CRG with CNV loss, such as ATP7B, the expression level was lower in the GBM sample than that in the normal sample, while the other CRG with CNV loss, CDKN2A, was significantly elevated in the GBM sample. Besides, other significantly differentially expressed CRGs showed very low frequencies of CNV.

A total of 209 patients from two GBM cohorts (TCGA-GBM and GSE83300 cohort) were merged. Detailed information on the 209 GBM samples was obtained. The prognostic significance of 12 CRGs was determined by Kaplan-Meier survival analyses and univariate Cox regression analyses. Five CRGs, including ATP7B, CDKN2A, DLD, MTF1, and SLC31A1, were identified as prognostic CRGs (Table S3). As shown in Figure 3A, a cuproptosis network systematically revealed the interactions and prognostic significances of CRGs in GBM.

Further, according to the expression profiles of the five prognostic CRGs, GBM samples were classified by consensus clustering algorithm. The consensus matrix heatmap showed that k = 2 is the optimal classification method, and the GBM samples were classified into CRGcluster A (number of samples = 119) and CRGcluster B (number of samples = 90) (Figure 3B). PCA indicated that the cuproptosis transcription profiles of the two subtypes were different (Figure 3C). Moreover, Kaplan-Meier curves showed that GBM samples in CRGcluster B had a longer OS than samples in CRGcluster A (chi-square test, P = 0.048; Figure 3D). By comparing the clinical characteristics of the two subtypes, no noticeable difference in age and sex was observed (Figure 3E). While the expression levels of most CRGs in CRGcluster A were higher than that in CRGcluster B (Figures 3E, F).

First, based on the ssGSEA algorithm, we obtained the relative contents of 23 kinds of immune cells in every GBM sample. And the relative contents of some kinds of immune cells, including MDSCs, CD56⁺ natural killer cells, macrophages, eosinophils, type 2 T helper cells, mast cells, monocytes, and CD56^- natural killer cells, significantly differed between the two subtypes (Figure 4A). As for the immune checkpoints, the expression level of PD-L1 in CRGcluster B was lower than that in CRGcluster A, while the expression levels of CTLA4 and PD-1 were higher than that in CRGcluster A (Figure 4B). Moreover, using the ESTIMATE algorithm, we obtained the TME scores of every GBM sample, including ImmuneScore, StromalScore, and ESTIMATEScore. ImmuneScore represents the content of the immune component, StromalScore represents the content of the matrix component, and ESTIMATEScore is the sum of the two. The differential analysis showed that TME scores were slightly higher in CRGcluster B, but no significant difference was observed (Figure 4C).

Using the R package “limma”, we identified 360 cuproptosis subtype-related DEGs. And via univariate Cox regression analysis, we identified 79 prognostic DEGs. Furthermore, according to the expression profiles of the 79 prognostic DEGs, GBM samples were classified by consensus clustering algorithm. The consensus matrix heatmap showed that k = 2 is the optimal classification method, and the GBM samples were classified into geneCluster A (number of samples = 96) and geneCluster B (number of samples = 113) (Figure 5A). Kaplan-Meier curves revealed that GBM samples in geneCluster B had a longer OS than in geneCluster A (chi-square test, P = 0.018; Figure 5B). By comparing the clinical characteristics of the two subtypes, no noticeable difference in age and sex was observed. While the expression levels of most prognostic DEGs in geneCluster A were higher than in geneCluster B (Figure 5C). The differentially expressed analysis of 12 CRGs between geneCluster A and geneCluster B is shown in Figure 5D.

According to cuproptosis subtype-related DEGs, a prognostic risk score model was developed. Figure 6A displays the distribution of GBM samples by different classification methods. Firstly, patients with GBM were randomly categorized into training and testing sets. Sample sizes of the two sets were about the same, with 105 GBM samples in the training set and 104 GBM samples in the testing set. Secondly, we conducted LASSO regression analysis of the 79 prognostic DEGs and obtained seven candidate genes for the risk score model (Figures 6B, C), followed by multivariate Cox analysis of the seven candidate genes, and finally obtained five target genes (PDIA4, PILRB, DUSP6, CBLN1, and PTPRN), including four high-risk genes (PDIA4, PILRB, DUSP6, and PTPRN) and one low-risk gene (CBLN1). Calculation formula: Risk score = (0.2988 × PDIA4) + (0.1705 × PILRB) + (0.2448 × DUSP6) + (0.3055 × PTPRN) + (-0.2095 × CBLN1).

The differential analyses indicated that the risk scores of geneCluster A and CRGcluster A were higher than that of geneCluster B and CRGcluster B, respectively (Figures 6D, E). In the training set, according to the median, GBM samples were classified into a low-risk cluster (n = 53) and a high-risk cluster (n = 52). With the increase in risk score, the OS of GBM patients gradually decreased, and the deaths steadily increased (Figure 6F). PCA clearly distinguishes the two risk clusters (Figure 6G). Survival curves revealed that GBM samples in the low-risk cluster had a longer OS than in the high-risk group (chi-square test, p< 0.001; Figure 6H). Furthermore, AUC values of 0.5-, 1.0-, and 1.5-year based on this model were 0.643, 0.709, and 0.751, respectively (Figure 6I).

Next, this model was verified in the testing set, GSE83300 cohort, and GSE74187 cohort, respectively (Figure 7). According to the median risk score, GBM samples of the testing set, GSE83300 cohort, and GSE74187 cohort were categorized, respectively. With the increase in risk score, the OS of GBM patients gradually decreased, and the deaths steadily increased (Figures 7A, E, I). PCA clearly distinguishes the two risk clusters (Figures 7B, F, J). In the testing set, Kaplan-Meier curves revealed that GBM in the low-risk cluster had a longer OS than in the high-risk group (chi-square test, p = 0.018; Figure 7C). In the GSE83300 cohort, although the result of survival analysis did not reach statistical significance, GBM samples in the low-risk cluster tended to prolonged OS (chi-square test, p = 0.121; Figure 7G). In the GSE74187 cohort, Kaplan-Meier curves revealed that GBM in the low-risk cluster had a longer OS than in the high-risk group (chi-square test, p = 0.007; Figure 7K). In the testing set, AUC values of 0.5-, 1.0-, and 1.5-year survival rates based on this model were 0.610, 0.671, and 0.708, respectively (Figure 7D). In the GSE83300 cohort, the AUC values of 0.5-, 1.0-, and 1.5-year survival rates based on this model were 0.676, 0.731, and 0.718, respectively (Figure 7H). In the GSE74187 cohort, the AUC values of 0.5-, 1.0-, and 1.5-year survival rates based on this model were 0.619, 0.731, and 0.741, respectively (Figure 7l).

Using RT-qPCR assay, the mRNA expression levels of the five model-related genes in five pairs of GBM and adjacent tissues were detected. Compared with the paired adjacent specimens, the expression levels of PDIA4, DUSP6, and PTPRN were upregulated; however, the expression levels of PILRB and CBLN1 were downregulated in GBM specimens (Figures 8A–E). To validate the protein expression levels of the five model-related genes, a WB assay was conducted. Consistently, compared with the paired adjacent specimens, the protein expression levels of PDIA4, DUSP6, and PTPRN were upregulated, and the protein expression levels of PILRB and CBLN1 were downregulated in GBM specimens (Figures 8F, G).

The risk score in the IDH1 mutant cluster was lower than in the IDH1 wild-type cluster (P< 0.0001; Figure 9A). Results of univariate and multivariate regression analyses revealed that both prognostic risk score and IDH1 mutation status were independent prognostic factors (Figures 9B, C). In addition, stratification analyses were performed to evaluate if the model had broad applicability. Kaplan-Meier curves revealed that GBM in the low-risk cluster always had a longer OS than in the high-risk group, and p< 0.001 in the subcluster of age younger than 60, p = 0.005 in the subcluster of age older than 60, p< 0.001 in subcluster of male, p = 0.015 in subcluster of female, p< 0.001 in subcluster of IDH1 mutation, and p = 0.003 in subcluster of IDH1 wild-type (Figures 9D–I).

Results of correlation analyses between the risk scores and the contents of immune cells indicated that they were positively correlated in regulatory T cells, follicular helper T cells, neutrophils, resting NK cells, and M0 macrophages, while negatively correlated in eosinophils, M2 macrophages, monocytes, and activated NK cells, (Figure 10A). Moreover, correlation analyses between the contents of immune cells and the expression levels of 5 model-related genes showed obvious correlations between some kinds of immune cells and specific genes (Figure 10B). Furthermore, GBM samples in the low-risk cluster had a significantly lower StromalScore and ESTIMATEScore than samples in the high-risk group (Figure 10C). Finally, differentially expressed analyses of immune checkpoints between the high-risk cluster and the low-risk cluster were conducted. Results revealed that the expression levels of 15 immune checkpoints were lower in the low-risk cluster than in the high-risk cluster, for example, PD-1 and PD-L1 (Figure 10D).

Results of the correlation analysis between the CSC indexes and risk scores indicated that they were negatively correlated (R = -0.35, P = 6.5e-06); that is, the stem cell characteristics of GBM patients with low-risk scores were more significant (Figure 11A). It is generally believed that tumors with high TMB would respond better to immunotherapy and thus have a better prognosis. Based on mutation data of the TCGA GBM cohort, the differential analysis demonstrated that TMB was lower in the high-risk cluster than in the low-risk group (Figure 11B). Next, the result of the correlation analysis between TMB and risk score indicated that they were negatively correlated (R = -0.048, p = 0.049; Figure 11C). Interestingly, there was no significant correlation in the low-risk group (R = 0.044, p = 0.12; Figure 11D) and a negative correlation in the high-risk group (R = -0.012, p = 0.04; Figure 11E). To determine the specific distribution of somatic mutations, we constructed the waterfall diagrams of the two risk score groups, and the ten most frequently mutated genes in the high-risk cluster were PTEN, EGFR, TP53, TTN, NF1, MUC16, PIK3CA, LRP2, RYR2, and SPTA1 (Figure 11F), while, in the low-risk cluster were TP53, PTEN, TTN, EGFR, MUC16, ATRX, SPTA1, FLG, IDH1, and RYR2 (Figure 11G). The mutation frequencies of PTEN, EGFR, and NF1 in the low-risk cluster were lower than in the high-risk cluster, while the mutation frequencies of ATRX, IDH1, TP53, MUC16, and PIK3R1 in the low-risk cluster were higher than in the high-risk cluster (Figures 11F, G). Furthermore, common drugs or compounds were selected to detect the association between drug sensitivities and risk scores. Results of differential analyses of IC50 values between the two risk clusters indicated that IC50 values of bryostatin, midostaurin, mirdametinib, ponatinib, and tipifarnib in the high-risk cluster were lower than in the low-risk cluster. In comparison, the IC50 values of afatinib and elesclomol in the high-risk cluster were higher than in the low-risk cluster. Results suggested that drug sensitivity was somewhat associated with risk score (Figures 11H–N).

To improve the applicability of this model in the clinic, a nomogram containing clinical parameters (age, sex, and IDH1 mutation status) and risk score was developed (Figure 12A). Based on the nomogram, the AUC values of 0.5-, 1.0-, and 1.5-year ROC curves for the prediction of OS were 0.716, 0.727, and 0.763 in the training set, respectively (Figure 12B), and were 0.741, 0.775, and 0.734 in the testing set, respectively (Figure 12C). Moreover, calibration curves showed that the predicted results were very close to the ideal results in the two sets (Figures 12D, E). Next, we detected the predictive performance of IDH1 mutation status alone. Based on IDH1 mutation status alone, the AUC values of 0.5-, 1.0-, and 1.5-year ROC curves for the prediction of OS were 0.663, 0.609, and 0.607 in the training set, respectively (Figure 12F), and were 0.679, 0.693, and 0.686 in the testing set, respectively (Figure 12G). Results indicated that the nomogram has excellent predictive performance and is better than IDH1 mutation status alone.