Construction and assessment of the NRG-based signature
We first utilized the GTEx and TCGA databases to screen 1085 genes associated with glioma prognosis, which were then intersected with 2299 genes related to neddylation, resulting in a total of 108 NRGs associated with glioma prognosis (Fig. 2A). Genes exhibiting abnormal expression were filtered using |log2FC|=2 criteria (Fig. 2B). Hazard ratios of several NRGs such as DNAJB6, SMN1, and RPN2 were presented using forest plots (Fig. 2C). Next, multivariate Cox regression analysis and the LASSO regression algorithm were employed to identify the most suitable combination of NRGs for prognostic prediction (Fig. 2D, E). Consequently, 7 NRGs, including TOP2A, F2R, UST, HSPA1B, LGALS3BP, UROS, and OSBPL11, were selected for inclusion in the construction of the prognostic signature (Table 1).
Construction and assessment of the signature in TCGA
Using the prognostic signature, risk scores were computed for each sample, stratifying patients into low- and high-risk groups based on the median risk score of the training set. The distribution of risk scores and expression levels among glioma patients in the training set is depicted in Fig. 3A. Notably, the scatter plot distribution of survival time indicates a positive correlation between poorer prognosis and higher risk scores. Furthermore, analysis of risk scores, survival time, and expression levels across both testing and entire cohorts yielded results consistent with those observed in the training set (Fig. 3B, C). Kaplan-Meier analysis revealed significantly shorter OS times in the high-risk group compared to the low-risk group across all subsets (Fig. 3D-F). Interestingly, Kaplan-Meier analysis indicated significantly longer PFS among high-risk glioma patients compared to low-risk patients (Fig. 3G-I). These results underscore the robustness of our signature in predicting prognosis among glioma patients.
Assessment of the signature in CGGA
We further evaluated our signature in the CGGA325 and CGGA693 glioma datasets. Patients with high-risk scores exhibited significantly shorter survival times. A gene heatmap depicted the expression characteristics of the seven genes modeled across all patients (Fig. 4A, B). The prognostic analysis also aligns with the results from TCGA, showing significantly shorter OS in patients classified as high-risk group (Fig. 4C, D).
Stratified analysis of glioma
To ascertain the relationship between our prognostic risk score and clinical features, we conducted a correlation analysis of commonly used clinical characteristics, including age, gender, WHO grade, 1p/19q co-deletion status, IDH1 mutant status and ATRX mutant status. The results indicate a significant correlation between age, glioma grade, IDH1 mutation status, and 1p/19q codeletion status with our risk score (Fig. 5A-F). To determine the independent prognostic factors associated with glioma patient outcomes, we conducted univariate and multivariate Cox regression analyses on potential predictive factors. The results of the analysis revealed significantly increased hazard ratios for age, glioma grade, 1p/19q codeletion status, IDH1 mutation status, ATRX mutation status, and risk score (Fig. 5G). However, the multivariate regression analysis revealed that ATRX mutation status (P = 0.281) did not reach statistical significance, while other factors such as age (HR = 1.045, 95% CI = 1.030–1.060, P < 0.001), grade (HR = 2.110, 95% CI = 1.553–2.867, P < 0.001), 1p/19q codeletion status (HR = 0.352, 95% CI = 0.199–0.622, P < 0.001), IDH1 mutation status (HR = 1.837, 95% CI = 1.146–2.944, P = 0.011), and risk score (HR = 1.012, 95% CI = 1.003–1.020, P = 0.007) demonstrated statistical significance (Fig. 5H).
Nomogram for predicting glioma prognosis
The prognostic charts forecasting the 1-year, 3-year, and 5-year survival rates for the TCGA cohort displayed AUC (area under the curve) values of 0.870, 0.890, and 0.807, respectively (Fig. 6A). In the CGGA325 cohort, the AUC values for the 1-year, 3-year, and 5-year survival rates were 0.773, 0.851, and 0.864, respectively (Fig. 6B). Similarly, in the CGGA693 cohort, the AUC values for the 1-year, 3-year, and 5-year survival rates were 0.697, 0.735, and 0.726, respectively (Fig. 6C). The signature (AUC = 0.890) appeared to provide more precise prognosis prediction compared to clinical factors such as age (AUC = 0.804), grade (AUC = 0.772), 1p/19q codeletion status (AUC = 0.621), and IDH1 mutation status (AUC = 0.785) (Fig. 6D). Additionally, we computed the concordance index for these clinical variables spanning from 1 to 5 years, and the risk score exhibited the largest AUC, underscoring the importance of our prognostic signature (Fig. 6E).
Next, we developed a nomogram incorporating IDH1 status, risk score, grade, 1p/19q status, and age to forecast the prognosis of glioma patients (Fig. 6F). The nomogram exhibited predictive accuracy with a concordance index of 0.861 (95% confidence interval: 0.835–0.888) (Fig. 6G).
Functional enrichment analysis
After data compilation, we conducted GO, KEGG, and GSEA functional analyses separately. The top three relevant biological processes (BP) in the graphene oxide analysis were organelle fission, nuclear division, and embryonic organ development. Simultaneously, cellular components (CC) and molecular functions (MF) were associated with collagen-containing extracellular matrix, endoplasmic reticulum lumen, chromosomal region, as well as extracellular matrix structural constituent, glycosaminoglycan binding, and peptidase regulator activity, respectively (Fig. 7A). Moreover, KEGG pathway analysis revealed significant pathways such as “focal adhesion” and “pathways in cancer”, further elucidating the potential regulatory mechanisms underlying tumors (Fig. 7B). Finally, GSEA analysis was performed on both high-risk and low-risk groups to compare differences in enriched signaling pathways. Pathways such as anterior posterior pattern specification, embryonic, focal adhesion, graft-versus-host disease, and systemic lupus erythematosus seemed to be more associated with the high-risk group (Fig. 7C), while pathways including glutamate receptor signaling pathway, positive regulation of excitatory postsynaptic potential, regulation of neuronal synaptic plasticity, neuron projection membrane, and neurotransmitter receptor complex appeared to be more associated with the low-risk group (Fig. 7D).
Correlation analysis of tumor immunity and molecule alteration
We utilized the ESTIMATE algorithm and found that the high-risk group exhibited higher StromlScore, ImmuneScore, and EstimateScore compared to the low-risk group (Fig. 8A). Subsequently, we examined the composition of various immune cells and observed a higher proportion of macrophages (including macrophage M2) in the high-risk group (Fig. 8B, C). Furthermore, we conducted a comprehensive analysis of immune-related functional disparities. Significant differences were observed between the high- and low-risk groups in functions related to APC (antigen-presenting cell) co-inhibition, CD8 + T cells, HLA (human leukocyte antigen), and other immune-related functions (Fig. 8D).
We compared the differences in Tumor Mutational Burden (TMB) between the high- and low-risk groups and found that the total TMB was significantly higher in the high-risk group than in the low-risk group (Fig. 8E). Notably, the low TMB and low risk score group (L-TMB + low risk) exhibited the highest survival probability (Fig. 8F, G). The mutation landscape plot revealed that the high-risk group had fewer IDH1 mutations but more mutations in EGFR and PTEN, all of which are molecular features associated with poor prognosis (Fig. 8H). Finally, we screened for potential drugs that could aid in the treatment of glioma (Figure S1).
Analysis of gene expressions at the single-cell level
Further, we performed scRNA-seq analysis to further analyze the expression of the gene signature at a single cell level. The expression characteristics of the CGGA scRNA-seq dataset after quality control and filtering are showed in Fig. 9A. The top 1500 highly variable genes are shown in Fig. 9B. After PCA analysis, the top 20 PCs with P < 0.05 were selected for further analysis. Totally, 15 clusters of cells were visualized by the t-SNE dimensionality reduction algorithm (Fig. 9C). Moreover, 5 major cell types (astrocytes, monocytes, macrophages, T cells, and epithelial cells) were identified (Fig. 9C). As Fig. 9E illustrated, genes like LGALS3BP were expressed in most cell types, whereas HSPA1B was mostly expressed in macrophages and monocytes. Gene like UST and OSBPL11 was expressed very low in specific cell types.
HSPA1B promotes the proliferation, migration and invasion of glioma cells and macrophage polarization
We investigated the correlation between HSPA1B expression and the malignant phenotype of glioma cells. Utilizing small interfering RNA (siRNA), we downregulated the expression levels of HSPA1B in U251 and T98G cells (Fig. 10A). Subsequent CCK8 and colony formation assays revealed that the proliferation capacity of U251 and T98G cells was inhibited upon downregulation of HSPA1B expression (Fig. 10B, C). Transwell and wound-healing assays demonstrated that the invasive and migratory abilities of U251 and T98G cells were suppressed following HSPA1B knockdown (Fig. 10D, E). Moreover, immunofluorescence staining analysis showed that the fluorescence intensity ratio of CD163 and CD68 was considerably lower in the si-HSPA1B group than in the control group in both U251 and T98G cell (Supplementary Fig. 2A-B). The statistical analyses were revealed in (Supplementary Fig. 2C-D). These results showed that HSPA1B promotes tumor-associated macrophage polarization in glioma cells.