The characteristics of cell senescence-associated genes
A total of 39 cell senescence-associated genes were included in this study through difference analysis (Additional file 9: Table. S1). We analyzed the regulatory relationships associated with the expression of 39 cell senescence-associated genes in GBM patients, for example, EZH2 positively regulates PTTG1 expression, while PTTG1 has a significant positive effect on TACC3 expression, and the result is presented in Fig. 2A. In the meantime, to deepen our understanding, we researched the somatic mutation prevalence of cell senescence-associated genes among GBM. Among them, a total of 13 out of 39 cell senescence-associated genes have somatic mutations and the MATK had the highest mutation rate, which is 2%, and mainly missense mutations (Fig. 2B). The analysis straightforwardly represented the differences in mutation profiles of different cell senescence-associated genes in GBM. Moreover, the investigation of 39 cell senescence-associated genes exhibited that CNV-related mutations were widespread. CKD4, SOX2, TACC3, CDK6, NUAK1, SOCS1 and MYC showed widespread CNV amplification, while EZH2, HK3, KL, IGFBP1, IGFBP3, VENTX, CAV1, SORBS2, PTTG1, GATA4, TLR3, SIX1, CBX7, AAK1, ALOX15B and TRPM8 had CNV deletions (Fig. 2C). Meanwhile, the localization of the 39 cell senescence-associated genes on TCGA-GBM 23 chromosomes was determined (Fig. 2D). The expression of most cell senescence-associated genes in GBM tissues was different from that in normal tissues. Additionally, we found that the expression of 31 cell senescence-associated genes significantly varied among three subtypes (CL, MES and PN) of GBM (Additional file 1: Fig. S1 and Additional file 2: Fig. S2). In an exploration of the effect of 39 cell senescence-associated genes on the overall survival (OS) of GBM patients, we found that downregulation of MYC, LIMA1 and SIX1 expression and upregulation of PTTG1, EZH2, AURKA, TACC3, VEGFA, CENPA, IGFBP3 and SOCS1 expression exhibited a statistically significant effect on the overall survival of GBM patients (Fig. 2E).
Furthermore, the IHC staining of PTTG1 and MYC was downloaded from the HPA. The result showed that the protein encoded by PTTG1 was higher expressed in glioma than normal tissues (#Antibody CAB008373). The level of MYC had lower expression in glioma than normal tissues (#Antibody CAB010307) (Fig. 3A). IHC staining was performed to detect the representative PTTG1 and MYC protein levels in gliomas and peritumor tissues (15 cases) of The Second Hospital of Shandong University. Consistent with the previous results, the protein level of PTTG1 was significantly increased in glioma in comparison to peritumor tissue, while the level of MYC was decreased in tumor (Fig. 3B).
To explore the effect of PTTG1 and MYC on glioma cell survival and invasion, we transfected cells with siRNA to downregulate PTTG1 and MYC. Western blot was used to verify the efficiency (Fig. 3C). Then EdU assays were performed. The downregulation of PTTG1 resulted in significant decreases in the percentage of EdU positive cells in U251 cells 48 h after transfection, while knockdown of MYC showed the opposite trend (Fig. 3D). In order to investigate the influence of PTTG1 and MYC on the invasion of glioma cells, we conducted a 3D collagen spheroid invasion assay. Silencing PTTG1 reduced the area invaded by U251 spheroids relative to controls, while knockdown of MYC increased the area. (Fig. 3E). These results suggested that PTTG1 and MYC silencing significantly influenced glioma cell activity in vitro.
Characterization and TME in different potential molecular subtypes
The unsupervised cluster analysis was used to classify the GBM patients which were the GlioVis dataset. 3 potential subtypes, gene. cluster1 (C1), gene.cluster2 (C2) and gene.cluster3 (C3) were acquired (Fig. 4A). The result of the Log-rank test showed the difference in GBM patient survival between 3 potential molecular subtypes (P = 0.0027) (Fig. 4B), which indicated that the GBM patients in C1, C2 and C3 subtypes were different and that confirmed the 3 subtypes distinction is reasonable. The mRNAsi and EREG-mRNAsi of the C3 cluster were closer to the 1 in the stemness index analysis, which suggests that GBM in C3 has a notable resemblance to stem cells (Fig. 4C). However, further research is necessary, as some studies have shown that higher indices appear to be directly associated with the degree of progression and poor prognosis for a variety of cancers, which is not consistent with the results of this study. In the Sankey diagram, the interrelationship between 3 subtypes and clinical typology, and G-CIMP/NON-G-CIMP was identified (Fig. 4D). In the exploration of the TME between the 3 subtypes, we found discrepancies in immunization checkpoints and immune cell infiltration. For instance, the expression of CD80, CD86, LDHA and PTPRC were abundant in C2 subtypes (Fig. 4E). And T cell regulatory, T cell CD8, B cells naive, dendritic cell activated, and other immune cells were also remarkably abundant in C2 subtypes. While NK cells activated and macrophages M2 were a significant expression in C1 subtypes (Fig. 4F). Meanwhile, The violin diagram demonstrates the differences in stromal score ESTIMATE score and immune score between the 3 potential subgroups (Fig. 4G). In the GOBP and KEGG analysis of the biological pathways of GBM, it was found that the senescence pathway plays a major role in GBM and the results are shown in Fig. 4H , I.
Construction of the risk score and Acquisition of clinical prognostic factors
The WGCNA operations were performed on the combined GBM dataset to obtain the key modules most relevant to the clinical features (Fig. 5A). A total of 228 differential cell senescence-associated genes were obtained, while co-expression modules were identified (Additional file 10: Fig. S10). According to the heatmap of module-trait relationships, the ME blue and ME turquoise modules demonstrate the highest pertinence with clinical features (Fig. 5B). Univariate Cox regression algorithm was exerted to, preliminary acquisition of 62 genes (Additional file 11: Fig. S11) associated with GBM prognosis and the HR and P values of the 228 cell senescence-associated genes were calculated, and the result is shown in Fig. 5C. Next, the LASSO algorithm and multivariate Cox regression analysis were applied to determine the prognostic gene set of GBM and ultimately found 11 gene sets (Fig. 5D and E). Finally, 11 cell senescence-associated genes were used to construct the risk score for predicting the survival and prognosis of GBM patients. The cohort of GlioVis-GBM patients was systematically randomized to distinguish the train set (n = 293) and the test set(n = 292). The Kaplan-Meier analyses evidenced that the GBM patients with high-risk scores in the train set (P < 0.0001) as well as the test set (P = 0.0053) corresponded with less favorable survival (Fig. 5F). The 1-, 2-,3-,5- year AUC of the train set were 0.710, 0.782,0.802 and 0.864, and of the test set were 0.576, 0.621, 0.683 and 0.602, respectively (Fig. 5G). In exploring the relationship between survival status and risk score, the survival rate of GBM patients gradually decreased as the risk score increased in both the train set as well as the test set (Fig. 5H). The multivariate Cox regression analyses were utilized to assess clinical independent prognostic factors for GBM patients. The results show that age was an independent prognostic factor for GBM patients in the train set as well as the test set (Additional file 3: Fig. S3 and (Additional file 4: Fig. S4).
The Hybrid Nomogram and TME
The mutual relationships regarding the three potential subtypes, clinical typology, the G-CIMP or NON-G-CIMP and risk scores were reflected by the Sankey diagram (Fig. 6A). The hybrid nomogram can be used in the clinical administration of GBM patients due to its stable and accurate characteristics (Fig. 6B). The nomogram was also integrated into ROC to assessed the survival time and survival situation of the GBM patients. The 1-, 2-,3-,5- year AUC of the train set were 0.741, 0.800,0.857, and 0.849, and of the test set were 0.623, 0.725, 0.732 and 0.693 (Fig. 6C). The calibration curves were applied to evaluate and validate the ROC, and were shown in Fig. 6D. Meanwhile, the differentially expressed genes of the risk factors were displayed in Fig. 6E. In the valuation of the liaison between the risk score and different underlying subtypes of GBM patients, the Wilcoxon test revealed a remarkable discrepancy in the risk score between the C1, C2 and C3 and the statistical difference between the risk score and G-CIMP or NON-G-CIMP also was appearance. And also, the Wilcoxon test proved a variance of age in high-risk and low-risk scores (Fig. 6F). In the analysis of TME, several immune-related factors in the low-risk score were found to be notably abundant (Fig. 6G). About the immune checkpoint, CD40LG, CD8A, JKA1, JAK2, LDHB and others were more pronounced high expression in low-risk scores as well as LDHA and YTHDF1 were expressed in high-risk scores significantly (Fig. 6H).
Validation of the risk score
Apply the same arithmetic algorithm to construct the risk scoring model whose samples were obtained from CGGA datasets. Then, the Kaplan-Meier analyses ferreted out the expression of GBM patients with high-risk scores corresponded with poorer survival (P = 0.0075) (Fig. 7A). It is shown that the percentage weights of clinical prognostic factors in GBM patients with a high-risk score and low-risk score are in Fig. 7B. For example, in the high-risk score group, the percentage of GBM patients who were no-methylated was 60.2%. While in the low-risk score group, only 40% of GBM patients were no-methylated. After, we drew the hybrid nomogram according to the CGGA dataset (Additional file 5: Fig. S5). The 1-, 2-,3-,5year AUCs which were combined risk score or nomogram score, were 0.564, 0.650,0.627, 0.685, and 0751, 0.747, 0.758, 0.717, respectively (Fig. 7C). We also passed the correction curve for verification (Fig. 7D). For further verification, we also constructed risk scoring models based on TCGA datasets and the same result was obtained by the Kaplan-Meier test (Fig. 7E). Furthermore, the map tools package was leveraged to the tumor somatic mutations presented in high- and low-risk scores respectively. The result shows that the low-risk score presented a wider range of TMB than the high-risk score (Fig. 7F). Meanwhile, the tumor somatic mutations of several genes are rarely observed in the low-risk group but frequently observed in the high-risk group, for example, PTEN. PTEN, which is the main negative regulator of the pI3k/Akt pathway and was a tumor suppressor gene with phosphatase activity closely related to tumorigenesis [32, 33]. But interestingly, TP3, the somatic mutation rate was as high as 63% in the low-risk group, which is a tumor suppressor gene with the function of regulating cell division and proliferation. The association of TP53 mutations with the development of a variety of tumors has been well documented [34, 35]. Therefore, further research is needed.
Pan-cancer analysis and drug sensitivity
Pan-cancer analysis was used to assess similarities and differences in risk score models between different tumor types. We systematically evaluated TMB, MSI as well as the expression of CD274 among pan-cancer. The risk score was proactively correlated with TMB in BRCA, COAD, LGG, PAAD, STAD and THYM (P < 0.05), while the inverse correlation with TMB in KIRC, KIRP, LAML and UVM (P < 0.05) (Fig. 8A). For MSI, a positive correlation in STAD, DLBC, COAD, HNSC and THCA as well as a negative correlation in CHOL and KIRC, was defined (P < 0.05) (Fig. 8B). Additionally, the risk score was positively relevant to CD274 expression in ACC, COAD, HNSC, LGG, SKCM, THCA and negatively relevant with CD274 content in BRCA, CESC, HNSC, KIRC, LAML, LUSC, OV and PCPG(P < 0.05) (Fig. 8C). In addition, the mutual relationship between risk score and several immune cell infiltration as well as stemness indices were calculated, respectively (Additional file 6: Fig S6 and (Additional file 7: Fig. S7, (Additional file 8: Fig. S8). Using the immunotherapy cohort of advanced urothelial cancer (IMvigor210 cohort) to evaluate the impact of risk scores on predicting immunotherapy sensitivity. The Log-rank test also showed that GBM patients with a high-risk score were associated with poorer survival conditions (Fig. 8D). Then, integration risk scores were analyzed with immune checkpoint blockade (ICB) treatment studies. The results showed that the proportion of GBM patients of the high-risk score group in the response groups (CR and PR) was notably lower than in the low-risk score group, while the percentage of patients in the no/limited response groups (SD and PD) showed the contrary phenomena, pointing that the risk score could prove the response of GBM patients to ICB therapy. However, in exploring immune phenotypes in high- and low-risk scores, the desert phenomenon was more notable in the low-risk score, while the inflamed was seen more in the high-risk group (Fig. 8E). From combined risk score and tumor neoantigen burden correlation analysis, the GBM patients with low-risk scores together with a high neoantigen burden exhibited the most stretched-out survival time and the GBM patients with high-risk scores in connection with a low neoantigen burden had the worst survival situation (Fig. 8F). The Spearman correlation analysis was shafted to measure the value of the risk score to anticipate drug sensitivity for multiple types of cancer. Finally, 119 drugs for which the risk score and drug sensitivity significantly correlated, were obtained from the GDSC database. Subsequently, we selected the 50 most representative drugs for mapping. The risk score was the most significantly negatively sensitive to 5 drugs, including AZD5991, YK.4.279, Alisertib, Vinblastine and Eg5_9814; and the most significantly positively correlated with sensitivity to 5 drugs, including LJI308, AMG.319, CZC24832, PLX.4720 and PFI3 (Fig. 8G). Among them, the strongest drug sensitivity is LJI308. The study shows that LJI308 is a powerful selective inhibitor of RSK, which can inhibit the growth and proliferation of cancer stem cells [36]. In addition, the signal path targeted by the selected drugs was discovered. The relationship between drug sensitivity and risk score targeted the cell cycle, mitosis, microtubule, DNA replication and apoptosis regulation signaling were positive. On the contrary, drug sensitivity with negatively related to the risk score targeting the PI3K/mTOR signaling and chromatin histone methylation (Fig. 8H). In summary, the establishment of a risk score will be beneficial in exploring the facility and effective treatment strategies for GBM.
