All human specimens were meticulously selected based on stringent criteria to ensure the integrity and relevance of our research findings. The specimens comprised glioma tissues obtained from six patients undergoing surgical resection at the China-Japan Union Hospital. Only specimens from patients with histopathologically confirmed glioblastoma were included, ensuring the study focused on relevant disease pathology. Specimens were obtained from patients who had not received any form of chemotherapy or radiotherapy prior to surgery. This criterion was established to minimize the potential influence of prior treatments on the gene expression profiles of the tumor tissues. Human samples were collected with the approval of the China-Japan Union Hospital Ethics Committee. All patients signed informed consent under the declaration of Helsinki.

Glioma patients with the gene expression profile, clinical features, and survival information were acquired from the TCGA database. Testing sets were adopted from the CGGA database (13) and GlioVis database (Rembrandt and Gill cohorts) (14). Immunotherapy cohorts were downloaded from the GEO database (GSE35640 and GSE78220 cohorts) (15, 16) and one published study (17). The single-cell profiles of glioma were acquired from the GSE182109 accession in GEO database (18).

Go and KEGG databases were applied to perform function analysis in the clusterProfiler package (19). P value less than 0.05 is significant.

We manually collected cuproptosis-associated genes from literature (5) and GSEA database (GOBP_COPPER_ION_HOMEOSTASIS.v2022 and GOBP_RESPONSE_TO_COPPER_ION.v2022) (20). To assess the cuproptosis response in glioma, we analyzed differential cuproptosis-related genes between the glioma and normal tissues and established the cuproptosis score (CP-score) based on the differential genes using the ssGSEA (for bulk data) and Ucell (for single-cell data) algorithms.

To filter out the significant prognostic genes for GBM, we performed Cox and Lasso regression based on the differential cuproptosis-associated genes. Significant genes with (Cox regression< 0.05) were retained for further analysis. To improve the model’s robustness, we randomly divided the training set into internal training and testing sets at 50% and then tested the model in three external sets. Each potential model was evaluated using 1-, 3- and 5-year AUC. Finally, the cuproptosis scoring model (CP-model) was constructed using the mathematical methodology below:

ere n is the number of cuproptosis-related genes, exp is the gene expression profiles, and β is the Cox regression coefficient. GBM individuals were categorized by the CP-model. The predictive ability was assessed using survival analysis.

We assessed the hazard ratios (HRs) using Cox regression for CP-model and other clinical features. HRs and P-value were shown in the forest plot. We then constructed a nomogram to display the survival time of GBM individuals in indicated times. Calibration plot demonstrated the accuracy of nomogram. Finally, we evaluated the CP-model and five existing models using C-index, KM curves, and restricted mean survival time (RMST).

We applied CIBERSORT algorithm to assess the proportions of immune cell types (21). To assess immune response and tumor cell activity, we adopted 20 pathways from a published study (22) and calculated pathway activity using the ssGSEA with the gsva R package. ESTIMATE and TIDE algorithms were applied to evaluate tumor purity and T cell dysfunction (23, 24).

The GDSC database and pRRophetic package were leveraged to evaluate the drug sensitivity for each specimen in R (25). Therein, we calculated the IC50 with 10-fold cross-validation and default setting (26).

Glioma tissue samples were prepared as 4 μm sections. The slides were deparaffinized, rehydrated in a series of gradient ethanol, and recovered them by heating the tissue at 100°C in citrate buffer for 1 hour. The endogenous peroxidase activity was inhibited with 3% hydrogen peroxide for 10 minutes at room temperature. The slides were treated with primary antibodies overnight at 4°C. After washing with PBS, we then incubated it with goat serum at room temperature for 30 min. Each segment was rinsed with PBS and then incubated with DAB for 5 min. Two pathologists assessed the IHC staining. IHC was performed using antibodies against CD3 (ab16669, Abcam), FOXP3 (ab20034, Abcam), Tryptase (ab2378, Abcam), and CD163 (ab79056, Abcam).

The classification of glioblastoma patients into high-risk and low-risk groups was based on their CP-model score, which was computed from the expression levels of CP-model genes measured by quantitative PCR (qPCR). This approach allowed for a robust stratification of patients. qPCR assays were conducted to measure the expression levels of the genes comprising our CP-model. RNA was extracted from glioblastoma tissue samples using TRIzol reagent, and cDNA synthesis was carried out using a High-Capacity cDNA Reverse Transcription Kit. qPCR was performed using SYBR Green PCR Master Mix in a Real-Time PCR Detection System. The CP-model score for each patient was calculated based on the expression levels obtained from the qPCR results using the formula derived from our model. The ΔΔCt method was used to analyze the relative changes in gene expression normalized to an endogenous reference gene and relative to a control.

To inspect the cuproptosis dysfunction in glioblastoma, we collected 65 cuproptosis-related genes from the published study and GSEA database and performed the differential analysis between glioblastoma and normal tissue. We found 31 differential expressed cuproptosis-related genes, including 26 up-regulated and 5 down-regulated genes ( Figure 1A ). To comprehensively analysis the correlation among these differential genes, we categorized them into two clusters based on hierarchical clustering and constructed a regulatory network. We observed a significant positive association of DLAT with DLD (cor = 0.703), and a negative association of MAP1LC3A with SORD (cor = -0.369) within the cluster A ( Figure 1B ).

We further constructed the CP-score to evaluate cuproptosis response in GBM. The CP-score was remarkably higher in the GBM patients versus the normal individuals ( Figure 1C ), and significantly associated with the advanced grade ( Figure 1D ). Rembrandt and Gill datasets were used to validate the results ( Figures 1E, F ). Furthermore, we performed enrichment analysis of the differential expressed cuproptosis-related genes and observed that they were enriched in multiple functions and pathways, such as copper regulation and homeostasis, metal transition pathways and mineral absorption ( Figures 1G, H ).

As tumor microenvironment (TME) are involved in tumorigenesis, we then assessed the CP-score with infiltrated immune cells in GBM. We found that CP-score were significantly correlated with several cell infiltrations ( Figure 1I ), such as positive correlation with M2 macrophages ( Figure 1J ), negative correlation with CD8⁺ T cells ( Figure 1K ).

To further assess the correlation between CP-score and TME in single-cell level, we acquired 60094 cells from GEO database under the GSE182109 accession number after quality control. We categorized these cells into 17 clusters ( Figure 2A ), and then annotated the cell types using the celltypist algorithm ( Figure 2B ). Moreover, the CP-score were calculated using the Ucell algorithm ( Figure 2C ). To distinguish each cell type, we labeled the representative markers of each cell type and demonstrated the top differential genes ( Figures 2D, E ). We observed that CP-score were significantly correlated with immune cell infiltration in GBM samples ( Figure 2F ).

Due to the importance of T cell and macrophages in the TME, we evaluated the CP-score in distinct T cell and macrophage subsets. We annotated the subtypes and observed the strong correlation among the cell diversity of CP-score which were in accordance with the bulk data ( Figures 2G–J ).

We further analyzed cuproptosis-related genes with glioblastoma prognosis using Cox regression. Twenty-five genes were associated with glioblastoma prognosis and ten-fold cross-validation was leveraged to calculate the deviance ( Figure 3A ). To identify the genes with the highest predictive value, we performed the Cox regression with LASSO penalty in the TCGA and internal training sets using 50% random sampling. We finally got five cuproptosis genes to establish the cuproptosis risk model (CP-model) as shown below:

We categorized GBM patients into two risk subgroups using the CP-model for the training and external sets using the same cut point obtained from the training set. As expected, we observed a shorter OS in the high-risk subtype from the TCGA training cohort as well as other five cohorts ( Figure 3B ). In addition, the predictive ability of CP-model was also evaluated using the time-dependent ROC analysis in the 1-, 3- and 5-year were all above 0.5, indicating the excellent predictive ability of CP-model ( Figure 3C ). Finally, the relative gene expression profile and risk score distribution were illustrated in Figure 3D . To monitor the robustness of the CP-model, the reliability of CP-model using the TCGA set was evaluated in three external sets using the same formula and cutoff from training cohort.

To evaluate the clinical application of the CP-model, we performed univariate and multivariate analyses with some clinical vicariates. We observed that the CP-model demonstrated a more reliable for glioblastoma prognosis compared with age and grade ( Figures 4A, B ). In addition, a nomogram was built to help clinicians to analyze the survival probability of GBM individuals by combining the CP-model and other independent factors ( Figure 4C ). In brief, these factors were weighted using a point system. Assigning values to the variables included drawing an upward-sloping straight line and then rescaling the total of the points for each variable. The total points were calculated by adding together all the individual variable scores. Illustrating a straight line to the outcome axis allowed us to calculate the 1-, 3-, and 5-year survival probabilities for GBM. The calibration plot revealed that the bias-corrected line was near the ideal curve, indicating high degrees of concordance between the predictive and actual points ( Figure 4D ). We further applied a time-dependent ROC curve to assess the survival probability, the AUC were 0.82, 0.86, and 0.79 in 1-, 3-, and 5-year, respectively ( Figure 4E ). The risk scores were a better predictor of outcomes than traditional factors ( Figure 4F ).

We chose five popular risk models from the literature to evaluate against our CP-model (27–31), and we plotted their ROC and KM curves to show how they compare ( Figures 5A, C ). C-index was calculated using the rms package in R. According to the findings, the C-index for CP-model was higher than other five predictive risk models ( Figure 5B ). The high-risk subgroup’s RMST across all six models is shown in Figure 5D . Overall, our results indicate that the CP-model was more accurate for predicting glioblastoma patients’ survival than individual prognostic markers and other existing models.

We applied a heatmap to show the expressed profiles of five cuproptosis-related genes and clinicopathologic factors ( Figure 6A ). The expression of SNCA, SORD and LCAT were inhibited, and ANKRD9 and SLC31A were relatively higher in the higher risk patients. Furthermore, most of the patients in grade III/IV were categorized as high-risk subgroups, reflecting the feasibly predictive ability of the CP-model. Risk scores were also associated with age, dead status, and grade in glioblastoma ( Figures 6B–D ). In line with risk stratification, the risk score was correlated positively with ANKRD9 and SLC31A, but negatively with SNCA, SORD and LCAT ( Figures 6E, F ). Altogether, our results indicated that the CP-model was independent of conventional clinical characteristics.

To mechanically detect the difference between the two subtypes, we applied GSEA analysis based on GO and KEGG databases. We observed that some GO functions were activated in the high-risk patients, such as adaptive immune response, B cell and T cell-mediated activation ( Figure 7A ). However, some GO functions were inhibited, including multiple amino acid transport functions ( Figure 7A ). In addition, multiple disease-related pathways were activated for KEGG analysis, but AMPK and JAK-STAT signaling pathways were inhibited ( Figure 7B ).

Next, we analyzed the principal component analysis (PCA) using the full set of genes ( Figure 7C ), on all cuproptosis genes ( Figure 7D ), and five selected cuproptosis genes from the CP-model ( Figure 7E ). We showed that the two risk subgroups were distinguishable based on the expression patterns of the selected cuproptosis regulatory genes.

To assess the TME in glioblastoma patients, we calculated immune response and tumor cell activity using the ssGSEA and assessed the infiltrated immune cell using the CIBERSORT algorithm. We found that 17 pathways and 14 types of cells were dramatically altered between the two subgroups ( Figures 8A, B ). Moreover, the risk score was positively correlated with seven cells, but negatively correlated with four cells ( Figure 8C ). Finally, we validated our findings using IHC by targeting immune cell markers ( Figure 8D ). We confirmed that CD3 targeting T cells were accumulated in low-risk group patients. However, FOXP3 targeting Tregs, tryptase targeting mast cells and CD163 targeting M2 macrophage were collected in high-risk group patients.

We next evaluated the tumor heterogeneity using ESTIMATE. High-risk subgroup patients have significantly elevated estimate, stromal, and immune scores but have low tumor purity ( Figure 8E ). Moreover, TIDE and T cell dysfunction, but not exclusion, were substantially lower in the high-risk subgroup ( Figure 8F ). Patients with low-risk and high TIDE scores exhibited the more extended OS than the others ( Figure 8G ). High-risk score correlated with high intratumor heterogeneity, high proliferation activity, high leukocyte fraction and low CTA score ( Figure 8H ). Thus, our findings show that aberrant immunological infiltration and tumor heterogeneity in glioblastoma may serve as prognostic markers and immunotherapy targets and may have substantial therapeutic consequences.

Cancer immunity cycle is a set of sequential events that killing cancer cells (32). The high-risk subgroup patients had a low score in steps 2 3 and 5, compared with the low-risk ones ( Figure 9A ). Conversely, the high-risk group presented low score in steps 1, 4, and 6 ( Figure 9A ). These findings indicated that high-risk patients had a better ability to identify tumors and induced T cells to tumors, but a weaker ability to infiltrate immune cells and kill tumors (32).

Next, we assessed the immunotherapeutic prediction of the CP-model using the SubMAP algorithm. High-risk score patients would get a better effect from anti-MAGA-A3 treatment ( Figure 9B ). However, we do not observe significant responsiveness in anti-PD-1/PD-L1/CTLA4 treatments. We further detailed the immunotherapeutic responsiveness in a different cohort from GSE35640 (targeting MAGE-A3 immunotherapy), IMvigor210 (targeting PD-L1 immunotherapy) and GSE78220 (targeting PD-1 immunotherapy). Consistent with the SubMAP results, the risk score could be negatively associated with the anti-MAGA-A3 response. Higher score patients presented worse outcomes in the GSE35640 cohort ( Figure 9C ). In addition, the risk score was positively associated with ICIs, such as PD-1, PD-L1, CTLA4, HAVCR2, LAG3 and TIGHT ( Figure 9D ). Interestingly, the CP-model could predict the response to the combination treatment of anti-PD-1 and anti-CTLA4 ( Figure 9E ). Since chemotherapy is a common way of anticancer, we analyzed the CP-model with chemotherapy response. We found that Pazopanib, Parthenolide, NVP.TAE684, Cisplatin, Docetaxel and GSK269962A would be more effective in high-risk patients ( Figure 9F ).