Ubiquitin-proteasome system-mediated ubiquitination modification patterns and characterization of tumor microenvironment infiltration, stemness and cellular senescence in low-grade glioma

Landscape of ubiquitination modification regulators in low-grade glioma

Twenty ubiquitination regulators (5 E1 activating enzymes, 5 E2 conjugating enzymes, 5 E3 ligases, and 5 DUBs) were comprehensively investigated in this research. The entire design of this research is described in Supplementary Figure 1A. The dynamic, invertible course of ubiquitination modification adjusted by regulators as well as their underlying biological roles for proteins were depicted in Figure 1A. Differential expression analysis of 20 ubiquitination regulators revealed that all of these regulators were markedly different in normal vs. LGG samples (Normal=400, LGG=523; Figure 1B). We next described the occurrence of copy number variations and somatic mutations of 20 ubiquitination regulators in LGG. Among the 506 patients, the mutation occurrence of 20 ubiquitination regulators was 1.38% (7 mutations; Figure 1C). USP44 had the highest mutation occurrence, followed by USP51, whereas E2-conjugating enzymes had no mutations in LGG patients. Additional investigation exhibited a strong mutation co-occurrence correlation between UBA7 and ATG7, as well as TRIM21 and ATG7, along with TRIM21 and UBA7 (Supplementary Figure 1B). The analysis of CNV alteration occurrence revealed a widespread CNV alteration in the ubiquitination regulators (Figure 1D). Figure 1E depicts the site of CNV alterations in ubiquitination regulators on chromosomes. We were able to successfully differentiate LGG samples from normal samples depending on the expression of the 20 ubiquitination regulators (Figure 1F). To determine if the aforementioned genetic variations impacted the expression of ubiquitination regulators in LGG patients, we explored the correlation between ubiquitination regulator mRNA expression levels and CNV and discovered that the expression levels of most ubiquitination regulators positively correlated with their copy numbers (Supplementary Figure 1C). Tumorigenesis is a complicated process, and CNV changes cannot completely account for the differences in ubiquitination regulator expression. Other factors such as methylation and transcription can influence gene expression.

Figure 1. Landscape of ubiquitination modification regulators in LGG. (A) Graphic abstract of this research. (B) Differential expression of ubiquitination modification regulators. (C) Mutation occurrence of 20 ubiquitination modification regulators. (D) CNV alteration occurrence of ubiquitination modification regulators. (E) Chromosomal sites of altered CNV in the ubiquitination modification regulators. (F) Principal component analysis of the 20 ubiquitination modification regulators in normal and LGG patients.

The above investigations showed remarkable heterology in the expression of ubiquitination regulators between normal and LGG samples, suggesting that an unbalance in the expression of ubiquitination regulators plays a critical role in the progression of LGG.

Ubiquitination modification patterns mediated by 20 regulators

Univariate Cox analysis confirmed the prognostic values of 20 ubiquitination regulators in low-grade glioma patients (Supplementary Figure 2A). The comprehensive landscape of ubiquitination regulator interactions and regulator relations, as well as the prognostic meaning for LGG patients, was illustrated with the ubiquitination regulator network (Figure 2A). To investigate the correlation among the regulators, we estimated pairwise connections between the expression of 20 ubiquitination regulators in LGG and discovered that positive associations outnumbered negative correlations (Supplementary Figure 2B). We discovered that not only the ubiquitination regulators in the homo-functional class exhibited a remarkable connection in expression, but there was also a strong relation between E1, E2, E3, and UDBs. Notably, we also observed that all ubiquitination regulators were significantly linked to immune cell infiltration and biological processes related to the regulation of the tumor microenvironment (Supplementary Figure 2C, 2D). Thus, cross-talk among the regulators of E1, E2, and E3 and UDBs could be crucial for the generation of distinct ubiquitination modification patterns and immune characteristics among LGG patients.

Figure 2. Biological characteristics of two ubiquitination modification modes. (A) The interplay between ubiquitination modification regulators in low-grade glioma. (B) Survival analyses for the two ubiquitination modification modes based on 523 patients with low-grade glioma from TCGA-LGG cohort. P<0.001. (C) GSVA analysis exhibits the activation states of biological pathways in two ubiquitination modification modes. (D) Abundance of Immune cell in the two ubiquitination modification modes. (E) The relationship between stemness index and the ubiquitination modification regulators.

Next, we used consensus clustering to classify the ubiquitination modification patterns of LGG based on the expression levels of 20 ubiquitination regulators. After unsupervised clustering, 326 LGG patients were found in UPS-cluster A, while the remaining 197 patients were found in UPS-cluster B (Supplementary Figure 3A–3C and Supplementary Table 1). Survival analysis of the two ubiquitination modification patterns indicated that UPS-cluster A has a better prognosis (Figure 2B). Furthermore, the CGGA-LGG dataset (as the validation cohort) was used to conduct the clustering analysis with the same method to validate the efficacy of unsupervised clustering. Notably, the same outcomes were gained, demonstrating the efficacy of our clustering (Supplementary Figure 3D–3F). Notably, there was a clear distinction in transcriptome expression patterns of ubiquitination modification regulators between the USP-clusters in the TCGA-LGG and CGGA-LGG cohorts (Supplementary Figure 3G–3H).

Different patterns of ubiquitination modification related to immune infiltration

To evaluate the biological behavior of the different ubiquitination modification patterns, we used GSVA analyses (Figure 2C). Samples in UPS-cluster A exhibited marked enrichment of pathways linked to propanoate metabolism and butanoate metabolism, while samples in UPS-cluster B showed significant enrichment of the cell cycle, nucleotide excision repair, DNA replication, and mismatch repair pathways.

Numerous investigations have demonstrated that ubiquitination modifications perform an essential role in the formation of TME and the DNA damage response [14, 15]. As a result, we first compared the differences in immune and stromal scores between ubiquitination modification patterns. The UPS-cluster B group had higher immune and stromal scores than the UPS-cluster A group (Supplementary Figure 4A, 4B). Subsequently, we identified distinctions in the DNA damage response-related biological processes among the ubiquitinated modification clusters. These biological processes, such as DNA damage repair, DNA replication, and mismatch repair, were markedly boosted in the UPS-cluster B group compared to the UPS-cluster A group (Supplementary Figure 4C). Furthermore, when comparing the UPS-cluster B group to the UPS-cluster A group, both pro- and anti-tumor immune signatures were upregulated, indicating that the ubiquitination modification pattern has a dual effect on anti-tumor immunity (Supplementary Figure 4D). The UPS-cluster B group had a higher abundance of immune cell infiltration compared to the UPS-cluster A group (Figure 2D). However, the UPS-cluster B group with higher immune cell infiltration had a worse prognosis. We speculated that the UPS-cluster B group’s DNA damage response-related phenotype suppressed immune cells’ anti-tumor capacity. Furthermore, the large amount of stroma in UPS-cluster B may confine immune cells around tumor cells to the periphery and prevent them from infiltrating into the tumor core, resulting in immune exclusion. Therefore, the UPS-cluster A group presented an immune-inflamed phenotype, and the UPS-cluster B group presented an immune-exclude phenotype. Given that biological processes like DNA damage repair and replication are linked to the CSC phenotypes and cellular senescence phenotypes, and the UPS plays a critical role in the regulation of the cell-cycle, we looked into the regulatory patterns of ubiquitination modification on the stemness feature and cellular senescence.

Distinct patterns of ubiquitination modification associated with stemness feature

Mounting studies have shown that ubiquitination is essential for the self-renewal, maintenance, differentiation, and tumorigenesis of CSCs [16]. Stemness indices (si) describe the differentiation status of stem cells. We obtained four stemness indices (mRNAsi, EREG-mRNAsi, mDNAsi, and EREG-mDNAsi) from previous studies to comprehensively analyze the stemness characteristics of LGG stem cells. Further analysis revealed marked correlations between regulator expression and stemness indices (Figure 2E). Moreover, the stemness of cancer stem cells is significantly linked to the prognosis of LGG patients (Supplementary Figure 4E). Distinctions in stemness indices were observed between the two ubiquitination modification patterns; particularly, the dedifferentiation phenotype was most apparent in the UPS-cluster B, whereas the differentiation phenotype was most apparent in the UPS-cluster A (Supplementary Figure 4F).

To further investigate the correlation between cancer stem cell phenotype and TME in the two ubiquitination modification patterns, weighted correlation network analysis (WGCNA) was used to identify designated phenotype-related module genes. The phenotypes were defined using the ESTIMATE score and the mRNAsi. Differentially expressed genes (DEGs) were found using a cutoff criterion of |logFC| > 1 and P< 0.05, and these DEGs were used to construct a scale-free system. The scale-free plot indicated that 3 is the best parameter for transforming the adjacency matrix into a scale-free topology (Supplementary Figure 5A–5B). Finally, 12 modules (merged dynamic) were gained by merging modules whose distance was less than 0.25 (Figure 3A). The module Eigengenes were then determined, and the correlation between the module Eigengenes and the designated phenotypes was calculated (Figure 3B). Surprisingly, the blue module correlated the most with ESTIMATE scores and the mRNAsi index. Therefore, we chose hub genes in the blue module for further investigation and discovered 18 stemness phenotype-associated genes and 88 TME phenotype-associated genes (Figure 3C, 3D). Importantly, 18 of these two kinds of genes overlap, indicating that ubiquitination modifications have cross-reactivity with TME and stem cell regulation.

Figure 3. Characteristics of ubiquitination modification-related phenotypes (CSC stemness). (A) Hierarchical clustering dendrograms of identified co-expressed genes in modules. (B) Correlations between the gene modules and clinical traits. (C, D) The hub genes in blue modules. (E) Survival analyses for the two CSC-clusters based on 523 patients with low-grade glioma from TCGA-LGG cohort. (F) GSVA enrichment analysis in the two CSC-clusters. (G) Abundance of Immune cell in the two CSC-clusters.

To explore the underlying biological behavior of the ubiquitination modification modes, we performed an unsupervised consensus clustering analysis of these 18 overlapped genes and further classified the samples into the corresponding cancer stem cell subtypes. As with the clustering results of the ubiquitination modification pattern, LGG patients were clustered into 2 distinct stem cell phenotypes (Supplementary Figure 5C–5E). Patients with CSC-cluster B had a better prognosis than patients with CSC-cluster A (Figure 3E). The expression patterns of ubiquitination modification regulators were significantly different in these two stem cell clusters, which again suggest that the two ubiquitination modification patterns have different effects on the tumor stem cell phenotype (Supplementary Figure 5F). The GSVA enrichment analysis showed immune-related biological processes such as leukocyte transendothelial migration, chemokine signaling pathways, B cell receptor signaling pathways, and cytokine-cytokine receptor interaction were significantly upregulated in the CSC-cluster A group (Figure 3F). Moreover, the CSC-cluster A group exhibited prominent enrichment of immune pathways (CD8 T effector, immune checkpoint, and antigen processing pathways) (Supplementary Figure 5G), as well as massive infiltration of immune cells (both pro- and anti-tumor immune cells) (Figure 3G), confirming again the intersectionality of ubiquitination modification regulating cancer stemness and immunity. The CSC-cluster B showed the opposite pattern (relatively few immune cells and weak enrichment of immune pathways). Importantly, the stromal-related pathways such as EMT and pan-fibroblast TGF response signaling pathways were also significantly enhanced in the CSC-cluster A, which can prevent the immune cell from infiltrating into the tumor core, resulting in a worse prognosis for the CSC-cluster A group.

Distinct patterns of ubiquitination modification associated with cellular senescence

Cellular senescence is a stress response induced by damage, which can result in the release of multiple cytokines, chemokines, and proteinases, ultimately leading to tumor immune microenvironment remodeling [17]. To further investigate the correlation between the cellular senescence phenotype and TME in the two ubiquitination modification patterns, we obtained 278 cellular senescence genes from CellAge and then 270 differentially expressed cellular senescence genes by differential analysis (p <0.05). Surprisingly, similar to the ubiquitination modification pattern clustering results, unsupervised clustering analysis based on these 270 cellular senescence-associated genes clustered LGG patients into two distinct cellular senescence phenotypes (Senescence-cluster A and B; Supplementary Figure 6A–6C). Patients with senescence cluster A had a better prognosis than patients with senescence cluster B (Figure 4A). The expression modes of ubiquitination modification regulators were greatly distinct in the two senescence clusters (Supplementary Figure 6D). Furthermore, patients in cellular senescence cluster A predominantly exhibited ubiquitination modification pattern A, whereas patients in cellular senescence cluster B predominantly exhibited ubiquitination modification pattern B, which suggests that the two ubiquitination modification patterns have different effects on the tumor cellular senescence phenotype (Figure 4B).

Figure 4. Characteristics of ubiquitination modification-related phenotypes (cellular senescence). (A) Survival analyses for the Senescence-clusters based on 523 patients with low-grade glioma from TCGA-LGG cohort. (B) The different ubiquitination modification patterns distribution between the Senescence-clusters. (C) GSVA analysis in the Senescence-clusters. (D) Enrichment score of biological process in the Senescence-clusters. (E) Abundance of Immune cell in the Senescence- clusters. (F) KEGG enrichment analysis of the differential expression genes between UPS- clusters.

The GSVA enrichment analysis showed that the immune-related biological processes such as leukocyte transendothelial migration and cytokine-cytokine receptor interaction and stromal-related pathways like ECM receptor interactions were significantly upregulated in the Senescence-cluster B group, whereas they were relatively decreased in Senescence-cluster A (Figure 4C). We further identified differences in the enrichment of stromal-related biological processes among the cellular senescence phenotypes. Stromal-related biological processes such as angiogenesis, EMT, and Pan−F−TBRS were markedly increased in the senescence-cluster B group compared to the senescence-cluster A group. Moreover, the Senescence-cluster B group exhibited prominent enrichment of immune pathways (CD8 T effector, immune checkpoint, and antigen processing pathways) as well as abundant infiltration of immune cells (both pro- and anti-tumor immune cells), confirming the intersectionality of ubiquitination modification regulating cellular senescence and immunity (Figure 4D, 4E). The Senescence-cluster A displayed the inverse trend (relatively few immune cells and weak enrichment of immune pathways).

Construction of ubiquitination modification signature

To verify the underlying biological course of the 2 ubiquitination modification modes, we recognized 216 DEGs linked to the ubiquitination modification modes using the limma package (Supplementary Figure 6E). KEGG analysis of DEGs revealed enrichment mainly in cell cycle, cellular senescence, cell adhesion molecules, ECM receptor interactions, and antigen processing and presentation pathways, confirming that ubiquitinated modifications regulate the tumor microenvironment through tumor stem cell phenotype and cellular senescence phenotype (Figure 4F). To confirm the adjustment mechanism of ubiquitination modification on the TME, we conducted the Lasso method on 216 DEGs to gain 37 signature genes of ubiquitination modification mode and then classified patients into distinct gene clusters relying on an unsupervised cluster algorithm of these genes (Supplementary Figure 7A). Aligned with the clustering results of ubiquitination modification regulators, the unsupervised clustering based on signature genes divided LGG patients into 2 classifications, which were dubbed UPS gene clusters A and B (Supplementary Figure 7B–7D). This outcome proved that LGG does have two distinct ubiquitination modification modes. To confirm the relationship of ubiquitination modification modes to the TME, we investigated the enrichment level of common TME signatures in the UPS gene clusters. Gene-cluster B had markedly increased stromal activity, as demonstrated by angiogenesis and upregulation of EMT signatures (Supplementary Figure 7E). Meanwhile, both pro- and anti-tumor immune signatures were abundant in UPS gene cluster B (Supplementary Figure 7F). Similarly, the patients belonging to gene cluster A (immune-inflamed phenotype) had a better prognosis (Supplementary Figure 7G). These outcomes were consistent with previous analysis and showed the importance of ubiquitination modifications in shaping the TME of LGG, which can categorize LGG patients as having an immune-inflamed phenotype or an immune-exclude phenotype.

In consideration of the heterogeneity and complexity of ubiquitination modifications, we needed to be able to accurately evaluate ubiquitination modification modes in individual LGG patients. We identified 7 independent prognosis-related signature genes using a multivariate Cox analysis of 37 ubiquitination modification mode signature genes (Supplementary Table 2). We then created a prognosis-related signature gene-based score system by Lasso analysis to assess ubiquitination modification patterns in individual patients with LGG, and Supplementary Table 3 shows 6 genes used to construct the score and their coefficients; the model was dubbed the UM-score (Supplementary Table 4). We discovered the UM-score of the UPS-cluster B group was considerably higher than that of the UPS-cluster A group (Figure 5A). In consistence, gene-cluster B had a considerably higher UM-score than gene-cluster A (Figure 5B). To evaluate the clinical significance of the UM-score further, we split patients into high- and low-UM-score subgroups using the media score. Patients with a low UM-score had a significant survival advantage (Figure 5C).

Figure 5. Construction of ubiquitination modification characteristic signature. (A, B) Distinction in the UM-score between ubiquitination modification modes and gene clusters in TCGA-LGG cohorts. (C) Survival analyses for the UM-score groups based on 523 patients with low-grade glioma from TCGA-LGG cohort. (D) The forecast value of UM-score in patients among the TCGA-LGG cohorts (AUC: 0.867, 0.896, and 0.895; 1, 2, and 3-years overall survival). (E) Multivariate Cox regression analysis, which included the factors of UM-score, Age, Grade, Gender, IDH1 status and patient outcomes in the TCGA-LGG cohort. (F) Survival analyses for the high- and low-UM-score groups based on 592 patients with low-grade glioma from CGGA-LGG cohort. (G) Multivariate Cox regression analysis, which included the factors of UM-score, Age, Grade, Gender, IDH1 status and patient outcomes in the CGGA-LGG cohort. (H, I) Distinction in the UM-score between CSC-clusters and Senescence-clusters in TCGA-LGG cohorts.

The AUCs of the ROC analysis for the UM-score were 0.867, 0.896, and 0.895 at 1, 2, and 3 years of OS, respectively (Figure 5D). We used the multivariate Cox method with patient clinical features such as age, gender, grade, and the IDH1 mutation to see if the UM score could serve as an independent prognostic factor. We discovered that the UM-score was a reliable and independent prognostic biomarker for assessing patient outcomes (Figure 5E; HR =0.43, 95% CI 0.27–0.67, p<0.001). The dependability of the UM-score was validated utilizing 553 CGGA-LGG patients. Aligned with above results, the low-UM-score group had a better prognosis than the high-UM-score group (Figure 5F) and multivariate Cox method also affirmed that the UM-score could be used as an independent prognosis factor (Figure 5G; HR = 0.37, 95% CI 0.25–0.55, p<0.001). These findings suggest that the UM-score can assess the ubiquitination modification modes and forecast the prognosis of LGG patients.

The role of the UM-score in assessing tumor microenvironment and immunotherapy

We further analyzed the correlation between UM-score and the two ubiquitination modification-related phenotypes. We discovered that UM-score of CSC-cluster A was markedly higher than CSC-cluster B (Figure 5H). Moreover, Senescence−cluster B had markedly higher UM-score than Senescence−cluster A (Figure 5I). Given the two ubiquitination modification-related phenotypes play a significant role in tumor immune microenvironment, we speculated the UM-score could be used to assess the immunal traits in patients with LGG. To explore the biological traits difference between UM-score groups, we used GSVA analysis on the two groups (Figure 6A). Samples in the high-UM-score group exhibited markedly enrichment of the cell cycle, nucleotide excision repair, DNA replication as well as mismatch repair pathways, whereas the samples in the low-UM-score group were on the contrary. The high-UM-score group had higher immune activity and stroma activity (Figure 6B and Supplementary Table 5), and the UM-score was considerably positively linked to immune scores (Supplementary Figure 8A; R = 0.42, p < 2.2e-16) as well as stroma scores (Supplementary Figure 8B; R = 0.46, p < 2.2e-16). Stromal-related biological processes such as angiogenesis, EMT and Pan−F−TBRS were significantly enhanced in high-UM-score group compared to low-UM-score group (Supplementary Figure 8C). Moreover, the high-UM-score group exhibited prominent enrichment of immune pathways (CD8 T effector, immune checkpoint, and antigen processing machinery) as well as substantial infiltration of immune cells (both pro- and anti-tumor immune cells; Figure 6C). While the low-UM-score subtype patients had a higher enrichment of Fibroblast Growth Factor Receptor 3 (FGFR3), the RTK/RAS pathway, and the PI3K pathway (Supplementary Figure 8D). These findings indicate that the UM-score can evaluate the tumor immune microenvironmental traits of LGG patients.

Figure 6. UM-score in the role of immunotherapy. (A) GSVA enrichment analysis in the UM-score subgroups. (B) Distinction in the immune scores and stromal scores between the UM-score subgroups in TCGA-LGG cohorts. (C) Abundance of immune cell in the UM-score subgroups. (D) Distinction in the TIDE score between the UM-score subgroups in TCGA-LGG cohorts. (E) Distinction in the MSI score between the UM-score subgroups in TCGA-LGG cohorts. (F) Distinction in the T-cell dysfunction score between the UM-score subgroups in TCGA-LGG cohorts. (G) Distinction in the T-cell exclusion score between the UM-score subgroups in TCGA-LGG cohorts. (H) Distinction in the TIDE score between the UM-score subgroups in CGGA-LGG cohorts. (I) Distinction in the MSI score between the UM-score subgroups in CGGA-LGG cohorts. (J) Distinction in the T-cell dysfunction score between the UM-score subgroups in CGGA-LGG cohorts. (K) Distinction in the T-cell exclusion score between the UM-score subgroups in CGGA-LGG cohorts.

Given that the expression of most immune checkpoints is significantly enhanced in the high-UM-score group, we speculated that the UM-score could be used to assess the immunotherapy responses in patients with LGG. TIDE scores have been proven to accurately predict the immunotherapeutic response. Patients with higher TIDE scores had a poorer immunotherapeutic response, implying immunotherapy was less likely to benefit them. In this study, TIDE scores were lower in patients with a low UM score, suggesting that ICI treatment was more effective on patients in the low UM score group (Figure 6D). To identify the source of the distinctions in immunotherapeutic response, we compared microsatellite instability (MSI), T-cell exclusion, and T-cell dysfunction scores between the UM-score groups (Supplementary Table 6). The low-UM-score group had higher scores for MSI and T-cell dysfunction, whereas the high-UM-score group had higher scores for T-cell exclusion (Figure 6E–6G). These results were verified by the CGGA cohort (Figure 6H–6K).

The utility of the UM-score across cancer types

Due to the previously mentioned close connection between UM-score and immunotherapy response, we subsequently investigated the role of the UM-score in various cancers. We discovered that the UM-score was linked to the prognosis of multiple kinds of cancer, as shown in Supplementary Figure 9A. Notably, UM-score had the highest HR on LGG prognosis, implying that ubiquitination modifications perform essential parts in the progression of low-grade glioma. MSI and TMB were meaningful biomarkers for immune checkpoint blockade therapy response. Radar plots revealed a significant correlation between UM-score and TMB in 12 of 33 cancers (Supplementary Figure 9B). Subsequently, we looked at the correlation between MSI and UM-score and discovered that ESCA had the strongest negative correlation (Supplementary Figure 9C). Supplementary Figure 9D shows that the levels of PD-L1 expression were significantly related to the UM-score. These findings confirm the scoring system’s ability to accurately predict immunotherapy response. Furthermore, we discovered a different relationship between the three biomarkers and UM-score in some cancer types. This phenomenon could be explained by the fact that cancers differ in terms of immune infiltration. We further investigated the relationship between the content of 22 immune cells and the UM-score, and discovered that the ratio of M1 to M2 macrophages correlated with the UM-score of the majority of cancer types (Supplementary Figure 9E). Furthermore, with the exception of ACC, CESC, DLBC, ESCA, KICH, OV, THYM and UCS, we discovered a link between UM-score and stemness index in 24 cancers (Supplementary Figure 9F).

All these results indicate that ubiquitination modifications play an indispensable role in shaping the tumor microenvironment.

Data collection

Low-grade glioma transcriptome data were obtained from TCGA-LGG dataset. The RNAseq transcriptome data for normal tissues were gained from GTEx. CNV data were gained from Xena. CGGA-LGG transcription data were downloaded from the CGGA website for validation purposes. 119 ubiquitination modification regulators (Supplementary Table 7) were obtained from previous studies [3, 43–46]. Based on the differential analysis of the 119 ubiquitination modification regulators, we respectively selected the top 5 ranked E1, top 5 ranked E2, top 5 ranked E3 and top 5 ranked DUB based on |Log₂FC| for further analysis.

Unsupervised consensus clustering on LGG patients

Unsupervised clustering method was performed to recognize the ubiquitination modes as well as the phenotypes associated with these patterns (CSC stemness and cellular senescence).

Gene set variation analysis

In non-parametric and unsupervised methods, GSVA was performed to assess the biological behaviour difference in distinct LGG patient subgroups (UPS-clusters, CSC-clusters, Senescence-clusters and UM-score subgroups).

Evaluation of immune infiltration features

The ssGSEA method was performed to calculate relative abundances of immune cell per sample. The genomic information used to label immune cell subtype was gained from Charoentong’s research, which annotated human immune cell subtypes.

Construction of LGG ubiquitination modification score (UM-score)

Among 37 ubiquitination modifications pattern signature genes, multivariate Cox method was performed to confirm 7 independent prognosis-related signature genes, which we selected for the development of a scoring system with the LASSO regression method. Ultimately, six genes and their coefficients were confirmed by the minimum criteria. The UM-score was calculated using the formula:

$$\text{UM-score}={\Sigma }_{\text{i}=1}^{\text{n}}{\text{Coef}}_{\text{i}}\ast {\text{x}}_{\text{i}}$$

Statistical analysis

The limma R package was applied to identify DEGs. The Spearman method was applied to conduct correlation analysis. The Wilcoxon rank sum test was applied to calculate the statistical difference between the 2 clusters. R tool was applied for all statistical tests.

Availability of data and materials

The data used in the study are described in detail in “Data collection and processing”.