Abstract Glioblastoma (GBM), a highly heterogeneous and aggressive brain tumor, presents significant clinical challenges due to its frequent recurrence and poor prognosis. In this study, we employed high-dimensional weighted gene co-expression network analysis (hd-WGCNA) and single-cell transcriptomic analysis to investigate the molecular heterogeneity of GBM. We identified functional gene modules associated with tumor cell subpopulations exhibiting highly malignant traits, particularly linked to proteasome dysregulation. Intercellular communication analysis revealed extensive interactions between malignant tumor subpopulations and tumor microenvironment (TME), highlighting critical crosstalk with tumor-associated macrophages (TAMs) and T cells. Using machine learning, we developed risk scores based on these malignant gene modules, which effectively stratify GBM patients by prognosis and treatment response, particularly in relation to immunotherapy. Furthermore, we systematically evaluated pathway enrichment, genomic variations, and drug response differences across risk groups. Finally, we validated the oncogenic role of PSMC2, a key gene in the proteasome complex, demonstrating its role in promoting GBM progression through cell proliferation, invasion, and epithelial-mesenchymal transition (EMT). Our findings provide novel insights into GBM heterogeneity, prognosis, and therapeutic strategies, suggesting PSMC2 as a potential therapeutic target. Keywords: Prognosis, Glioblastoma, Tumor microenvironment, PSMC2 Subject terms: Cancer microenvironment, CNS cancer Introduction Glioblastoma (GBM) is the most common malignant primary brain tumor and a significant cause of mortality. The median overall survival for glioblastoma is approximately 15 months^[30]1,[31]2. Despite aggressive treatment modalities, glioblastoma exhibits treatment resistance and is highly prone to recurrence, with tumors always recurring at or near the resection cavity within the radiotherapy field^[32]3. A defining feature of glioblastoma is its remarkable heterogeneity, which contributes to its aggressive behavior and treatment evasion. This heterogeneity manifests in various cell populations, each with distinct genetic and phenotypic characteristics, influencing tumor growth and response to therapies^[33]4. The tumor microenvironment (TME) plays a crucial role in shaping this heterogeneity, characterized by an immunosuppressive landscape populated by pro-tumor immune cells such as exhausted T cells, tumor-associated microglia, tumor-associated macrophages (TAMs), and myeloid-derived suppressor cells (MDSCs)^[34]5,[35]6. The interactions between the diverse tumor cell populations and the TME promote tumor progression and resistance to therapy, yet the underlying mechanisms remain poorly understood. Further exploration of the relationship between GBM heterogeneity and the TME is essential for developing more effective treatment strategies. Emerging studies have explored the prognostic classification of GBM using transcriptomic or genomic data. These classifications are based on core biomarkers or enriched pathway signals, each associated with distinct clinical features^[36]7–[37]9. Nonetheless, significant discrepancies persisted among previous classifications. A meaningful and clinically applicable molecular taxonomy to guide clinical decisions remained elusive, as the cellular and molecular mechanisms underlying treatment responses in the subtypes remained unclear. PSMC2 (proteasome 26S subunit ATPase 2) is a 19S regulatory component of the 26S proteasome that regulates the selective degradation of intracellular proteins^[38]10. Existing studies have shown that PSMC2 is upregulated in tumors such as hepatocellular carcinoma, ovarian cancer, glioma, breast cancer, and melanoma, and is associated with poor prognosis^[39]11–[40]14. Limited research indicates that PSMC2 is involved in malignant pathways, including the PI3K/AKT/mTOR pathway, apoptosis resistance, migration, and cell cycle progression^[41]15–[42]18. This indicates the important role of PSMC2 in tumor progression. However, its role in GBM has not been sufficiently studied. In this study, we utilized hd-WGCNA to identify gene modules associated with highly malignant cell subpopulations by integrating single-cell sequencing data. These subpopulations reflect the intrinsic heterogeneity of glioblastoma (GBM), which contributes to its aggressive nature and treatment resistance. Additionally, we constructed a cell communication network between the gene modules defining these highly malignant subpopulations and the tumor microenvironment (TME). Using a machine learning approach, we developed prognostic features related to these malignant subpopulations and systematically assessed their predictive value for patient outcomes in GBM. We also explored the role of these features within the context of the TME, particularly their ability to predict responses to immunotherapy and their correlation with immune cell infiltration. Furthermore, we characterized different risk groups of GBM patients based on pathway enrichment, genomic variation, and differential drug responses, highlighting the diverse molecular landscape of the tumor. Finally, through a series of analyses and experiments, we validated the oncogenic role of PSMC2 in driving malignancy in GBM. Materials and methods Public data sources acquisition and pre-process The gene expression matrix and the corresponding clinical data for GBM were collected from the Gene Expression Omnibus (GEO, [43]https://www.ncbi.nlm.nih.gov/geo/), the Genotype-Tissue Expression (GTEx, [44]https://gtexportal.org), The Cancer Genome Atlas (TCGA, [45]https://portal.gdc.cancer.gov), and the Chinese Glioma Genome Atlas (CGGA, [46]http://www.cgga.org.cn), including TCGA-GBM and CGGA-GBM. Raw data were converted to transcripts per kilobase million (TPM) and log2-transformed (TPM + 1). Single cell sequencing analysis Single-cell sequencing data from SCP503, which includes 23 tumor specimens from 7 GBM patients, were collected and preprocessed. For quality control, the raw gene counts matrix was filtered according to the original research, normalized using the Seurat R package^[47]19, and pre-processed with the SCTransform function. The UMAP method was used to reduce the dimensionality of the data. The cells were annotated according to the original metadata^[48]20. In order to investigate the heterogeneity of tumor cells, tumor cells originating from tumor samples were further extracted. The R package “hdWGCNA” was employed to identify gene functional modules associated with tumor heterogeneity, with the optimal soft threshold set to 10, as the average connectivity of the topological network is most suitable when this threshold is set to 10. The R package “CytoTRACE” was utilized to define the differentiation potential of tumor cells^[49]21, while the R package “InferCNV” was used to calculate copy number variations in Tumor cells. Subsequently, CNV scores were further computed. CellCHAT (version 1.1.0) was employed to identify communication networks between cells. The interaction information between ligands and receptors was sourced from CellChatDB^[50]22. Development of risk model To further investigate the potential prognostic value of genes within the brown module, we identified optimal prognostic biomarkers among these genes using the Lasso Cox regression model, implemented with the “glmnet” R package^[51]23, in the TCGA cohort. Subsequently, we constructed the risk scores using the regression coefficients derived from the Lasso Cox regression. To select the most suitable lambda value and avoid overfitting, we applied ten-fold cross-validation. The risk scores were formulated as follows: graphic file with name 41598_2024_83571_Article_Equa.gif GBM patients were stratified into high- and low-risk groups based on the median risk scores. The Kaplan–Meier method was utilized to evaluate survival disparities between these groups. Furthermore, time-dependent receiver operating characteristic (ROC) curve analysis was performed to assess the precision of the risk scores in forecasting GBM prognosis. Additionally, we validated the model in the CGGA cohort to ensure its robustness and predictive accuracy. Exploration of tumor microenvironment infiltrations The ssGSEA algorithm was used to quantify the relative abundance of immune infiltration. The ‘IOBR’ R package serves as an integrated tool for TME analysis^[52]24. The Estimate scores for each sample in the TCGA-GBM cohort, as well as the immune cell estimations using the Cibersort deconvolution, were computed using the ‘deconvo_tme’ function of ‘IOBR’. Specific parameters were set to their default values without any modifications. Additionally, the Tumor-Associated-Macrophage (TAMs) signature and Immune Checkpoint Blockade (ICB) resistance score signatures were obtained from the ‘IOBR’. These signatures were employed in estimating the relative scores for each risk group using the ‘ssGSEA’ method. Drug sensitivity analysis and Immunotherapy response predictions To analyze drug sensitivity, we used the “oncoPredict” R package selecting “GDSC2” as the training data. The parameter “batchCorrect” was set to “eb”. In brief, we predicted the corresponding drug sensitivity values for patients in the TCGA-GBM cohort using drug sensitivity data and gene expression information from “GDSC2”. Additionally, we calculated the Pearson correlation between risk scores and drug sensitivity values to explore the association between risk scores and drug response. We set the threshold for the absolute value of the correlation to be greater than or equal to 0.3 to identify associated drugs. For immunotherapy analyses, due to limited research on immunotherapy for glioblastoma, we downloaded microarray and survival data, along with corresponding response information, from “IMvigor210CoreBiologies” to validate the predictive power of the risk scores for treatment response. Pathway enrichment analysis The Gene Ontology (GO) analysis of differentially expressed genes (|logFC| ≥ 1, P < 0.05) among different risk populations were performed using the “genekitr” R package. The Gene Set Enrichment Analysis (GSEA) of different risk groups was based on the hallmark gene set using the “genekitr” R package. The gene set variation analysis (GSVA) for KEGG ([53]http://www.kegg.jp/kegg/kegg1.html) pathways^[54]25–[55]27, based on the TCGA-GBM dataset, was performed using the “gsva” R package. The hallmark gene set enrichment of single-cell sequencing data was measured using the single-sample rank-based scoring approach (singscores). The R package ‘aPEAR’ was used to analyze Gene Ontology Biological Process (GOBP) pathway enrichment for genes co-expressed with PSMC2 in the TCGA-GBM dataset^[56]28. Cell culture Human GBM cells (U87-MG, LN-229) and HEK-293FT cells were obtained from the American Type Culture Collection (ATCC, USA), subcultured, and frozen in our laboratory. GBM cells were cultured in DMEM medium (Vivacell, China) supplemented with 10% fetal bovine serum (Vivacell) and 1% penicillin–streptomycin (Genview, China) at 37 °C with 5% CO2. shRNA transduction The PLKO.1-based lentiviral shRNA technique was used to construct the GBM cells with stably expressed the PSMC2 shRNA or negative control shRNA. The shRNA plasmids were purchased from Tsingke Biotechnology (China). The shRNA plasmid, along with packaging plasmids, was transfected into the HEK-293T cells to produce lentiviral particles. The GBM cells were infected with lentiviruses expressing the shRNA constructs. The following shRNA sequences were used: PSMC2 shRNA#1,5ʹ-CAACGTAAAGCAGTTTGCCAA-3ʹ; shRNA#2,5ʹ-AAGCAAGTTGAAGATGACATT-3ʹ. MTT assay GBM cells suspension was diluted to 1000–2000 cells per well and seeded into 96-well plates. OD values were taken on days 1, 3, 5, and 7. 20 µL MTT (5 mg/mL, Genview) was added to to the designated well, and incubated for 4 h in 5% CO[2] at 37 °C. After incubation, the medium was removed, followed by the addition of 200 μL of DMSO (Sigma, USA) , and the mixture was shaken horizontally for 10 min. The absorbance at 540 nm was then measured to assess cell viability. Edu immunofluorescence In a 24-well plate, 2 × 10^4 cells were seeded and incubated for 24 h. Then, the cells were then treated with 1-ethynyl-2ʹ-deoxyuridine (Edu, Beyotime, China) for 4 h. The cells were washed three times with PBS and fixed for 30 min in 4% paraformaldehyde (PFA, Beyotime). The cells were then permeabilized for 15 min with 0.3% Triton X-100 (Sigma). Following permeabilization, the cells were incubated in the dark for 30 min with 100 µL Click reaction solution. After gentle washing, the cells were incubated in the dark for 20 min with a 200 µL 1 µg/mL DAPI (Beyotime) solution. Following three additional washes, 200 µL PBS was added, and the cells were ovserved and imaged using an inverted fluorescent microscope. Clonal formation assay The GBM cells were plated into 6-well plates at a density of 1000 cells per well with complete medium. Cells were cultured for 7–10 days and then fixed in 4% PFA (Beyotime) and subsequently stained with Crystal Violet stain. (Beyotime, China). Migration and invasion assay Cells were seeded in 24-well Millicell hanging chambers with a pore diameter of 8.0 µm (Corning, USA). The cells were grown in serum-free medium in the upper chamber and 10% FBS in the lower chamber. For the invasion assay, the upper chamber was coated with Matrigel (BD,USA), whereas no coating was applied for the migration assay. All cells were cultivated at 37 °C for 24 h for migration and 48 h for invasion. After incubation, the cells were stained with Crystal Violet (Beyotime) for microscopic imaging. Wound healing assay GBM cells in the logarithmic growth phase were seeded into a six-well plate at a density of 1 ×10^6 cells per well. Once the cells reached 90% confluence, a scratch was made using a 300-400 µm pipette tip to ensure uniform wound width as the start. The cells were then washed twice with PBS to remove debris and cultured in serum-free medium. The wound area was subsequently monitored and imaged at designated time points. Soft agar colony formation assay In each well of a 6-well plate, add 1 mL of 2× medium supplemented with 0.6% agarose (Sigma) and wait for it to harden into basic agar. The constructed upper layer includes 0.3% agar, and 10,00–2000 cells. After 2–3 weeks of incubation at 37 °C in 5% CO[2], the colonies were stained with 200 µL MTT (5 mg/mL, Genview) per well for 30 min, and the number of colonies was counted. Quantitative real-time PCR Total RNA was extracted using the RNAiso Plus protocol (Tiangen, China). M-MLV reverse transcriptase was then used to synthesize cDNA from 1 µg of RNA in each sample (Promega, USA). The amplification reaction conditions were set according to the manufacturer’s instructions (Promega). The qRT-PCR was performed in triplicate using the Roche LightCycler 96 real-time PCR equipment (Roche, Switzerland). B[2]M expression was served as a control. The primers used for real-time PCR are listed in Table [57]S1. Western blot Total proteins were extracted using strong RIPA cleavage buffer including the serine protease inhibitor phenylmethane sulfonyl fluoride (PMSF, Beyotime). The protein content was measured using a BCA protein detection kit (Dingguo, China). Protein samples were separated using 10% or 12% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE), and the protein was transferred to a PVDF membrane (Merck). At room temperature, the membrane was blocked with either 5% bovine serum albumin (BSA, Dingguo) or 5% skim milk in TBST for approximately 2 h, followed by incubation overnight at 4 °C with the primary antibody diluted to the recommended concentration. The membrane was then incubated with an HRP-conjugated secondary antibody (Life Technology, USA) at 37 °C for 2 h. Finally, protein bands were visualized using the ECL Western blot detection system with ECL reagent (Beyotime). The antibodies used for Western blot analysis were as follows: PSMC2 ([58]HA500404, Huabio, China), α-TUBULIN (11224-1-AP, Proteintech, China), MMP2 (10373-2-AP, Proteintech, China), and MMP9 (10375-2-AP, Proteintech, China). Statistical analysis Two continuous variables conforming to a normal distribution were compared using the t test. The Kruskal–Wallis test was used to compare differences among two or more groups. Survival curves were generated for each subgroup in the dataset using the Kaplan–Meier method, and statistical differences were assessed using log-rank tests. All statistical analyses were carried out using R version 4.3.0 ([59]https://www.r-project.org/). P-values were calculated on a two-sided basis. Correlations between variables were assessed through Pearson correlation analysis. Results Integrative scRNA-seq analysis reveals gene modules associated with malignant phenotype scRNA-seq is a powerful technique for investigating tumor heterogeneity. scRNA-seq data from SCP503 were collected and preprocessed. Principal Component Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction were conducted for data visualization (Fig. [60]1A). Based on the original study, cells were classified into immune, normal brain, and tumor categories (Fig. [61]1B). To explore heterogeneity within GBM tumor cells and their associated functional modules, hd-WGCNA was employed with a soft threshold of 10 (Fig. [62]1C), corresponding to optimal average connectivity. This analysis identified a total of 14 gene modules (Fig. [63]1D), each displaying distinct enrichment patterns in tumor cells. Notably, the green, red, and purple modules showed minimal expression variation across tumor cells, while the brown, turquoise, blue, and green-yellow modules exhibited varied expression patterns (Fig. [64]1E). Fig. 1. [65]Fig. 1 [66]Open in a new tab Integrative scRNA-seq analysis reveals gene modules associated with malignant phenotype. (A) Uniform Manifold Approximation and Projection (UMAP) graph of different cell clusters. (B) UMAP graph of different cell lineages. (C) The top-left panel illustrates the soft power threshold used to select a scale-free topology model with a fit greater than or equal to 0.8. The remaining three panels display the mean, median, and maximum connectivity of the topological network, respectively, for various minimum soft thresholds, reflecting the network’s connectivity. The average connectivity of the topological network is most stable when the lowest soft threshold is set to 10. (D) hd-WGCNA dendrogram plot of 14 modules. (E) Feature plots of the corresponding harmonized Module Eigengene scores (hMEs) in tumor cells. (F) t-SNE plot of CytoTRACE-predicted cell order (more to less differentiated), indicated by colour from blue to red. (G) InferCNV analysis revealing copy number variations (CNV) in tumor cells. (H) Heatmap of Pearson correlations between different modules and Cytoscores, CNV scores, and recurrence. (I–K) Gene Ontology (GO) analysis of the brown module genes. (P < 0.05, **P < 0.01, ***P < 0.001). Our next objective was to investigate tumor heterogeneity at the single-cell level, with a focus on characteristics associated with multiple malignant traits. CytoTRACE, a method for inferring cell trajectories through gene expression states, was employed to assess cell differentiation potential^[67]21. Using CytoTRACE, we visualized the heterogeneous differentiation potential of GBM cells (Fig. [68]1F). Genomic instability, a hallmark of malignancy that often increases with tumor progression, is characterized by copy number variations (CNV)^[69]29. We used InferCNV to analyze CNV in GBM cells at the cellular level (Fig. [70]1G), calculating CNV scores for each cell. We then correlated Cytoscores (calculated by CytoTRACE), CNV scores, and recurrence phenotype with different modules, identifying the brown module as closely associated with tumor dedifferentiation, CNV scores, and recurrence (Fig. [71]1H). Gene Ontology (GO) enrichment analysis further elucidated the biological significance of the genes within the brown module. Gene Ontology Cellular Component (GOCC) analysis highlighted associations with the proteasome, endopeptidase, and peptidase (Fig. [72]1I). Biological Process (GOBP) analysis revealed significant involvement in aerobic respiration, purine-related processes, and proteasome-related processes (Fig. [73]1J). Gene Ontology Molecular Function (GOMF) enrichment indicated significant enrichment in proteasome-activating activity (Fig. [74]1K). Exploring glioblastoma heterogeneous cellular communication in tumor microenvironment Next, we aimed to explore the association between the brown module and the tumor microenvironment (TME). Immune cells were further subdivided into macrophages and T cells based on classical markers (Fig. [75]2A,B). We observed significant enrichment of brown module harmonized module eigengenes (hMEs) in different tumor cell clusters, particularly in clusters “19”, “5”, and “13” (Fig. [76]2C). Therefore, these tumor cells were further defined as brown module-associated tumor cells, exhibiting characteristics of tumor dedifferentiation, recurrence, and genomic instability. Using CellCHAT, we constructed communication networks among different cell types. Our analysis revealed extensive interactions of brown module associated-tumor cells with normal brain cells, macrophages, and T cells. Importantly, brown module-associated tumor cells showed heightened cellular communication intensity compared to other tumor cells, particularly in interactions with macrophages (Fig. [77]2D). Cell signalling patterns revealed distinct outgoing communication pathways among different cell types, with brown module-associated tumor cells exhibiting the highest level of signal activation. Signals activated in brown module-associated tumors included PTN, CypA, CD99, COLLAGEN, THY1, NOTCH, among others, which are associated with malignant signalling pathways related to tumor proliferation, invasion and stemness maintenance. This reflects the malignancy of brown module-associated tumors (Fig. [78]2E). Fig. 2. [79]Fig. 2 [80]Open in a new tab Exploring glioblastoma heterogeneous cellular communication in tumor microenvironment. (A) UMAP graph of different cell types. (B) Feature plots of classical cell type markers. (C) Boxplot showing brown module hMEs across different tumor cell clusters. Clusters 19, 5, and 13 exhibit relatively high brown module hMEs. (D) Cell communication network diagram. The left panel represents the number of interactions, while the right panel represent the interaction weights/strength. Line thickness represents the relative strength of the signal. (E) Outgoing signalling patterns Heatmap analysed by CellCHAT. The colour of heatmap represents the relative strength of signals, with green indicating signal enrichment. Development and validation of risk scores based on brown module genes Next, we aimed to construct a prognostic model based on the brown module for risk stratification in GBM. Lasso–Cox regression was employed to determine the optimal λ value and further identified 13 risk genes (Fig. [81]3A,B), which were used to construct the risk scores. The Kaplan–Meier survival curve was utilized to validate the risk stratification capability of the risk scores in glioblastoma (Fig. [82]3C), showing that patients with risk scores above the median exhibited poorer prognosis (Fig. [83]3D). Time-dependent ROC curves were generated to assess the predictive ability of the risk scores for survival (Fig. [84]3E). The results indicated stable performance of the prognostic model, with 1-year, 3-year, and 5-year AUC values of 0.73, 0.80, and 0.87, respectively. Additionally, we validated the significant prognostic efficacy of the risk model for overall survival in the CGGA-GBM cohort (Fig. [85]3F), with AUC values of 0.81, 0.78, and 0.75 at 1-year, 3-year, and 5-year, respectively (Fig. [86]3G). Fig. 3. [87]Fig. 3 [88]Open in a new tab Development and validation of risk scores based on the brown module genes. (A) Lasso coefficient profiles of prognostic-related brown module genes. (B) The tuning parameter (λ) selection in the Lasso model was determined through tenfold cross-validation using minimum criteria. (C) Survival status, time distribution, and model genes expression of patients in the high- and low-risk scores groups. (D) Kaplan–Meier survival curves for glioblastoma patients stratified by different risk types in the TCGA-GBM cohort. (E) Time-dependent receiver operating characteristic (ROC) curve analysis at 1, 3, and 5 years showing the area under curve values (AUC) for overall survival of TCGA database. (F) Kaplan–Meier survival curves for glioblastoma patients according to different risk types in the CGGA cohort. (G) Time-dependent ROC curve analysis at 1, 3 and 5 years showing AUC for overall survival in the CGGA database. Immune landscape of different risk groups To assess the correlation between risk scores and the tumor microenvironment (TME) of GBM, we calculated the numbers of immune cells within patients’ TME using the ssGSEA method, revealing widespread immune cell infiltration in the high-risk group (Fig. [89]4A). Further analysis revealed positive correlations between estimate scores, immune scores, stromal scores, and risk scores, while a negative correlation was observed with tumor purity. Additionally, xCell stromal scores, TME scores, immune scores, and risk scores showed positive correlations (Fig. [90]4B). Interestingly, the analysis indicated that individuals in the high-risk group, despite belonging to the immune-hot subtype, had poorer clinical prognosis. It was speculated that although immune cells infiltrated extensively in the high-risk group, normal anti-cancer immune functions may have been depleted. Furthermore, the correlation between risk scores and various immune cell types was assessed using the CIBERSORT algorithm. As expected, risk scores showed positive correlation with NK cell resistance, regulatory T cells, and neutrophils, while displaying negative correlation with NK cell activation (Fig. [91]4C). Immune checkpoint markers PD-L1 and PD-1 were also significantly upregulated in the high-risk group (Fig. [92]4D,E), further confirming immune suppression in high-risk group. Fig. 4. [93]Fig. 4 [94]Open in a new tab Immune landscape of different risk groups. (A) Immune cell heatmap based on TCGA-GBM datasets. Red signifies higher immune cell infiltration, while blue indicates a lack of immune cell infiltration. (B) Heatmap of the correlation between the Estimate scores, xCell scores, and risk scores. (C) A scatter plot shows the correlation between immune cell infiltration significance and the risk scores, computed using the “Cibersort” method on TCGA-GBM datasets. (D,E) Boxplot of PD-1 or PD-L1 log2(TPM) expression levels of different risk groups. The numerical values at the top of the panel represent P-values (F,G) Boxplot of TAMs signature scores and ICB resistance signature scores of different risk groups. The signatures were calculated by ssGSEA algorithm. (H) Kaplan–Meier survival curves to different risk groups in the ICB treatment cohort (IMvigor210). (I) Boxplot illustrating the risk scores for different response groups. (J) Barplot depicting the proportion of PD-L1 responses in different risk groups. (P < 0.05, **P < 0.01, ***P < 0.001, Complete Response (CR), Partial Response (PR), Stable Disease (SD), Progress Disease (PD)). Tumor-associated macrophages (TAMs) play a crucial role in constituting the immune microenvironment responsible for tumor progression and immune suppression. In the TME of glioblastomas, TAMs play a key regulatory role. The TAMs scores for different risk groups were further calculated, revealing higher TAMs scores in the high-risk group (Fig. [95]4F). This indicates potential extensive infiltration of TAMs in the high-risk group. The efficacy of immune checkpoint blockade (ICB) therapy, an emerging and promising treatment approach, was further evaluated in terms of the response of different risk groups to treatment. A significantly higher ICB resistance scores were observed in the high-risk scores group compared to the low-risk scores group, indicating a potential resistance to ICB therapy in individuals with high-risk scores (Fig. [96]4G). These predictive results were subsequently validated in the IMvigor210 cohort (Fig. [97]4H–J). Enrichment analysis reveals biological evidence of poor prognosis for high-risk patients To investigate the biological mechanisms associated with the risk scores, we performed Gene Ontology (GO) enrichment analysis based on differentially expressed genes (DEGs) in different risk groups. Our analysis revealed the upregulation of cytokine activity and pathways related to the extracellular matrix in the high-risk group, which are crucial in modulating tumor metastasis (Fig. [98]5A–C). The Gene Set Enrichment Analysis (GSEA) of hallmark gene sets was further performed to explore the differences in signalling pathways enrichment between high- and low-risk scores groups. As a result, cancer-related pathways were significantly activated in high-risk group, such as Epithelial mesenchymal transition (EMT), KRAS signalling up, Glycolysis, TGF-beta signalling (Fig. [99]5D). Furthermore, Gene Set Variation Analysis (GSVA) was employed to compute the enrichment status of each patient for KEGG pathways. We observed that in high-risk patients, pathways closely associated with tumor progression, such as ECM receptor interaction and JAK-STAT signalling, were significantly activated, while they were suppressed in low-risk scores patients (Fig. [100]5E). In summary, our pathway analysis indicates the activation of pathways related to tumor metastasis, malignant progression, and immune regulation in high-risk populations. This provides a unique biological context for adverse prognosis in these individuals. Fig. 5. [101]Fig. 5 [102]Open in a new tab Enrichment analysis reveals biological evidence of poor prognosis for high-risk patients. (A–C) Gene-ontology analysis of differentially expressed genes in different risk groups patients. Salmon colour represents pathways upregulated in the high-risk group, while blue indicates pathways downregulated in the high-risk group. The horizontal axis represents the significance of pathway enrichment. (D) Gene set enrichment analysis (GSEA) comparing high-risk group to low-risk group. The display includes the top 6 significantly enriched in Hallmark database. (E) Heatmap of the activity of the top 20 KEGG pathways enrichment in different risk groups. Orange represents pathway activation, while blue represents pathway inhibition. Mutation and drug sensitivity analysis Different risk groups defined by risk scores exhibit distinct genomic mutation characteristics (Fig. [103]6A). IDH mutations, TP53 mutations, and MGMT promoter methylation are generally considered favourable prognostic factors in GBM^[104]30–[105]32. Consistently, our analysis reveals an enrichment of TP53 and IDH1/2 mutations in the low-risk group (Fig. [106]6B). Additionally, patients with wild-type IDH display significantly elevated risk scores (Fig. [107]6C). Moreover, the risk scores are notably increased in recurrent patients. Furthermore, patients with lower levels of MGMT promoter methylation also exhibit higher risk scores (P = 0.1) (Fig. [108]6D). Overall, our analysis establishes a robust association between risk scores and known prognostic features. Next, we explored the association between risk scores and drug sensitivity. Using bioinformatic methods, we predicted the correlation between risk scores and IC50 values of 198 drugs. Our analysis revealed a negative correlation between risk scores and 7 drugs, as well as a positive correlation with IC50 values of 25 drugs, including Vorinostat and Carmustine (|R| > 0.3) (Fig. [109]6E,F). This suggests distinct drug treatment strategies based on the risk scores. Fig. 6. [110]Fig. 6 [111]Open in a new tab Mutation and drug sensitivity analysis. (A) Mutations of top 15 genes on the basis of TCGA-GBM cohort. Each column represents each sample, and each row represents each gene. Mutation types are marked by unique color. The left panel represent the high-risk group and the right panel represent the low-risk groups. (B) Forest plots of mutations bias in high-risk group and low-risk group. (C) Boxplot of risk scores for different IDH1/2 mutation status. (D) Boxplot of risk scores for different sample types (Top panel) and different MGMT promoter statuses (Bottom panel). (E) Spearman correlation plot between risk scores and drug IC50 values, where red indicates a positive correlation with risk scores and blue indicates a negative correlation (Top panel). The bottom panel depicts drug sensitivity across different risk subgroups. (F) Boxplot of IC50 values in different risk groups. The top panel represents Vorinostat, while the bottom panel represents Carmustine (*P < 0.05, **P < 0.01, ***P < 0.001). Integrative single-cell sequencing and bulk RNA-sequencing analysis reveal the role of PSMC2 in promoting malignancy PSMC2 plays a crucial role in the risk scores. Analysis based on TCGA-GBM data revealed significant upregulation of PSMC2 in tumor tissues (Fig. [112]S1A), and Kaplan–Meier survival curves indicated its association with poor prognosis (Fig. [113]S1B). Further analysis at the cellular level unveiled predominantly high expression of PSMC2 in tumor cells (Fig. [114]7A). Furthermore, tumor cells were subset for CytoTRACE analysis, revealing extensive enrichment of PSMC2 within undifferentiated tumor cells (Fig. [115]7B). Additionally, we constructed a pseudotime trajectory of tumor cells, where our analysis highlighted a gradual increase in PSMC2 expression along the pseudotime trajectory, emphasizing its critical role in tumor progression (Fig. [116]7C,D). Supporting our analysis, pathway enrichment results revealed an extensive overlap between PSMC2 expression and pathways related to EMT and the maintenance of cancer stemness pathways, including WNT BETA CATENIN SIGNALING and HEGHOG SIGNALING (Fig. [117]7E,F, Fig. [118]S1C). Furthermore, using aPEAR on the TCGA dataset, we conducted Gene Ontology Biological Process (GOBP) pathway enrichment analysis of differentially expressed genes between glioblastoma patients with high and low PSMC2 expression. We observed significant upregulation of pathways associated with cell division in patients with high PSMC2 expression, while pathways related to neurons were downregulated (Fig. [119]7G). In summary, our analysis suggests the potential role of PSMC2 in promoting the malignancy of GBM. Fig. 7. [120]Fig. 7 [121]Open in a new tab Integrative single-cell sequencing and bulk RNA-sequencing analysis reveal the role of PSMC2 in promoting malignancy. (A) UMAP-density plot of PSMC2 expression levels in different cell lineages of glioblastoma. The greener colour indicates higher PSMC2 expression. Cells within the dashed line represent the same cell lineage. The chart includes annotations for cell types. (B) t-SNE plot of CytoTRACE predicted cell order (more to less differentiated), indicated by color from blue to red(left panel). t-SNE plot of PSMC2 expression levels(right panel). (C) Pseudotime trajectory plot of tumor cells. (D) The expression curve of PSMC2 across pseudotime. (E) UMAP-density plot of PSMC2 expression levels in tumor cells. (F) UMAP-density plot of pathway activity in tumor cells. (G) Gene Ontology Biological Process (GOBP) enrichment network generated by “aPEAR”. The nodes represent significant pathways, while the edges denote the similarity between them. They are coloured based on the normalized enrichment scores (NES). PSMC2 is involved in cell proliferation and colony formation in GBM cells To investigate the biological function of PSMC2 in GBM, we generated GBM cell lines with stable knockdown of PSMC2 using two different shRNA sequences (Fig. [122]8A). Cell viability was assessed using the MTT assay, revealing a significant inhibition of GBM cell proliferation upon PSMC2 knockdown (Fig. [123]8B). Additionally, colony formation assays demonstrated a notable reduction in the clonal capacity of PSMC2-depleted GBM cells compared to control cells (Fig. [124]8C,D). Further experiments with EDU staining showed a significantly lower EDU-positive rate in PSMC2 knockdown GBM cells compared to the control group (Fig. [125]8E,F). To evaluate the contribution of PSMC2 to GBM tumorigenesis, we performed a soft agar assay, which revealed a marked decrease in the size of colonies formed by PSMC2-depleted cells (Fig. [126]8G,H). All of the above results indicate that PSMC2 is closely related to GBM cell proliferation and colony formation. Fig. 8. [127]Fig. 8 [128]Open in a new tab PSMC2 is involved in cell proliferation and colony formation in GBM cells. (A) PSMC2 expression after PSMC2 knockdown was analyzed by Western blot and RT-qPCR in GBM cells. (B) Viability of GBM cells with PSMC2 knockdown. (C) Colony formation ability of GBM cells with PSMC2 knockdown. (D) The quantification of colony-formation ability was detected after knock down PSMC2 in GBM cells. (E) EdU staining after knockdown PSMC2 in GBM cells. (F) Quantification of EdU-positive cells was quantified following PSMC2 knockdown in GBM cells. (G) Colony formation ability after PSMC2 knockdown was assessed by soft agar assays in GBM cells. (H) Size statistics of colonies in panel (G). Scale bar = 200 μm. All data represent the mean ± SD from at least three independent experiments. Student’s t-test was performed to analyze significance; *P < 0.05, **P < 0.01, ***P < 0.001. PSMC2 is required for cell migration and invasion in GBM cells As one of the most malignant types of glioma cells, GBM cells showed powerful migration and invasion capabilities. In this study, we investigated the impact of PSMC2 knockdown on the migration and invasion abilities of GBM cells. The results from the wound healing assay demonstrated that at 0, 12, 24, and 48 h, the lateral migration rate of the introduced wound in the PSMC2 knockdown group was significantly lower compared to the control group (Fig. [129]9A,B). Moreover, the Transwell invasion assay revealed that PSMC2 knockdown significantly inhibited the invasion of GBM cells compared to the control group (Fig. [130]9C,D). Consistent with these findings, western blot analysis showed that PSMC2 knockdown led to decreased expression of the mesenchymal markers MMP2 and MMP9 (Fig. [131]9E). These results indicate that PSMC2 plays a critical role in promoting GBM cells migration and invasion, likely through its involvement in regulating the epithelial-mesenchymal transition (EMT) phenotype. Fig. 9. [132]Fig. 9 [133]Open in a new tab PSMC2 is required for cell migration and invasion in GBM cells. (A) The migration ability of GBM cells after knockdown of PSMC2 was determined by wound healing assay. (B) Cell healing rate for wound healing assay. (C) Migration and Invasion of GBM cells after PSMC2 knockdown. (D) Quantification of migration and invasion in panel (C). (E) Western blot analysis was performed to detect the expression of metastasis-related proteins after knockdown of PSMC2. Scale bar = 100 μm. All data represent the mean ± SD from at least three independent experiments. Student’s t test was performed to analyse significance; *P < 0.05, **P < 0.01, ***P < 0.001. Discussion Glioblastoma (GBM) represents the most prevalent intracranial tumor. Treatment modalities for glioblastoma include surgical resection, radiotherapy, and chemotherapy. However, tumor recurrence and progression are pervasive. Glioblastoma exhibits significant heterogeneity, leading to substantial variations in patient prognosis and treatment responses. Consequently, investigating individualized therapeutic strategies and understanding the heterogeneity of glioblastoma have emerged as a crucial research focus^[134]3,[135]5,[136]33. In our present study, we employed hdWGCNA to investigate tumor heterogeneity in GBM for the first time. By associating functional modules with tumor recurrence, dedifferentiation, and CNV, we identified a key gene module defining pro-malignancy traits. Further analysis revealed that genes within this module were significantly enriched in proteasome-related pathways, suggesting a crucial role of the proteasome in GBM progression. Supporting our analysis, in a siRNA screening analysis that identified relevant genes for GBM survival, 22% (12/55) were components of the 20S and 26S proteasome subunits^[137]34. Besides, proteasome inhibitors such as Bortezomib and MG-132 induce the accumulation of p21 protein and arrest cell cycle progression at the G2/M phase, highlighting the potent therapeutic potential of targeting the proteasome in GBM^[138]35. Emerging evidence increasingly suggests that the crosstalk between tumor cells and the tumor microenvironment plays a pivotal role in the progression of GBM. In our study, we identified tumor cells associated with the brown module as linked to multiple malignant traits. Additionally, these cells exhibited dysregulation of the proteasome pathway, highlighting the critical role of proteasome in GBM progression. Our cell communication analysis further supports extensive interaction between brown module-associated tumor cells and TAMs. The CellCHAT signalling network we depicted reveals the unique cellular signals of brown module-associated tumor cells, including the enrichment of signals such as CD99, PTN, and CypA. High expression of CD99 has been reported to be associated with glioma recurrence, immune checkpoint blockade (ICB) therapy resistance, and increased TAMs^[139]36. PTN secreted by glioblastoma stem cell induces tumor cell proliferation and promotes GBM progression through the PTN-PTPRZ1 axis, activating signaling pathways like β-catenin, PI3K-AKT-mTOR, and RAS-RAF-MEK-ERK^[140]37. In summary, the signaling network we constructed offers valuable insights into the crosstalk between heterogeneous GBM cells and the tumor microenvironment. Furthermore, risk scores constructed based on the brown module genes can effectively reflect the risk stratification of GBM. Through systematic analysis of the relationship between risk scores and the tumor microenvironment, we established the association between risk scores, TAMs, and the immunosuppressive microenvironment. Patients defined by risk scores exhibit a unique immune landscape, reflecting a preference for immune therapy response in the low-risk population. Additionally, we systematically analyzed pathway enrichment patterns, mutation features, and drug responses in patients with different risk scores, providing insights into understanding risk, tumor prognosis, and treatment selection. PSMC2 emerged as a key gene in the risk scores, with the highest kMEs in the risk model, suggesting a central role for PSMC2 in GBM. PSMC2 is a component of the proteasome complex and is essential for the assembly of 19S and 26S proteasomes^[141]10. Previous studies have demonstrated that knocking down PSMC2 promotes apoptosis of GBM cells^[142]17. Here, we found that the expression of PSMC2 is higher in GBM tumor tissues than in non-tumor tissues, and the high expression of PSMC2 is associated with poor prognosis in GBM patients, and promotes the tumor progression of GBM through cell proliferation, colony formation, and cell migration. Epithelial-mesenchymal transition (EMT) is a crucial process in tumor initiation and development, serving as a key factor in malignant tumor cell metastasis and influencing various aspects such as chemo sensitivity^[143]38, Studies have indicated that interfering with the EMT process is a vital strategy for inhibiting tumor progression. In our study, PSMC2 inhibits the proliferation and migration of GBM cells by decreasing the expression of mesenchymal markers MMP2 and MMP9, suggesting that PSMC2 could be used as a therapeutic target for GBM. This study is subject to several limitations. First, despite utilizing multiple analytical methods to analyze GBM heterogeneity, we lack experimental validation. Besides, the risk model we developed lacks validation with internal cohorts. In conclusion, our study unveils the complex landscape of gene modules associated with pro-malignancy phenotype in GBM, providing insights into tumor heterogeneity, prognosis, TME, and potential therapeutic strategies. Our results also suggest that PSMC2 is a promising candidate for further study of GBM and potential therapeutic targets. Supplementary Information [144]Supplementary Information.^ (1.6MB, pdf) Author contributions Jingsong Yan designed and conducted bioinformatics analyses. Nini Zhou, Manya Xiong performed the experiments. Shunqin Zhu provided technical supports and research ideas. Nini Zhou, Jingsong Yan and Manya Xiong analyzed the experimental data. Nini Zhou and Jingsong Yan wrote the manuscript. Shunqin Zhu revised the manuscript. Funding This work was supported by Chongqing Postgraduate Innovation and Entrepreneurship Project (No. [145]CYS22199); The Eaglet Program of Young Innovative Talents Cultivation ([146]CY230235); Chongqing Municipal Training Program of Innovation and Entrepreneurship for Undergraduate (Project S202310635405). Data availability All of the sequencing data of this study were collected from TCGA ([147]https://portal.gdc.cancer.gov), GEO ([148]https://www.ncbi.nlm.nih.gov/geo/), CGGA ([149]http://www.cgga.org.cn) and Broadstitute—Single-Cell—Portal ([150]https://singlecell.broadinstitute.org/single_cell). The corresponding data set is provided in the method. Code can be provided on reasonable request. Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Nini Zhou and Jingsong Yan contributed equally to this work. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-024-83571-5. References