We obtained RNA-seq data from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) for glioblastoma (GBM) and lower-grade glioma (LGG). The following criteria were used to choose participants: (1) Patients with a past pathological diagnosis of lower-grade glioma or glioblastoma; (2) Gene expression and clinical data are reported for each patient. A total of 692 patients were included in the analysis after screening. Half of the patients were randomly assigned to the training cohort, while the other half was assigned to the validation cohort.

The GENECARDS database (https://www.genecards.org/) was used to find genes relevant to ubiquitination. All ubiquitination-related genes were found by searching for “ubiquitination” in the search box. For further investigation, we extracted the top 100 most relevant genes.

Univariate COX regression was used to identify genes linked with patient survival in gliomas in order to investigate the prognostic significance of these ubiquitin-related genes. The analysis platform is R software (version 4.1.0), and the “Survival” R package is utilized for COX regression analysis.

To build a prognosis model of ubiquitination, researchers used Least Absolute Selection and Shrinkage Operator (LASSO) regression after identifying ubiquitination genes having prognostic value. Each model gene was matched to generate the relevant coefficient after achieving the optimal LAMDA value, allowing the risk score of various patients to be calculated: Risk score =∑_(I =1)^nmβ _I *(expression of ubiquitination associated gene I). Based on the median risk value as a cut-off, patients in various cohorts might be separated into high-risk and low-risk categories. The model’s prognostic usefulness was then investigated using survival analysis to measure the prognostic difference between the two groups. The prognostic model’s 1, 3, and 5-year ROC curves were also generated to assess the model’s accuracy and robustness.

To avoid bias, univariate and multivariate COX regression studies were done to further examine the model’s prognostic efficacy. To discover independent prognostic markers, risk scores and other clinical characteristics (age, sex, and Karnofsky performance score) were included in the analysis.

To analyze the single-cell data acquired from the GEO databases, we utilized the “Seurat” software (version 1.3.1). The PCA dimension reduction approach, as well as the t-Distributed Stochastic Neighbor Embedding (tSNE) method, were used to identify cell subclusters in dataset GSE162631. Using feature genes and the “SingleR” packages, the cells were re-clustered. As a result, the expression of several cells was demonstrated.

The tumor microenvironment and GBM/LGG-infiltrating immune cells were then assessed in silico. Based on bulk RNA-seq data, ESTIMATE is an algorithm for predicting the presence of invading stromal/immune cells in the tumor microenvironment. ESTIMATE was able to generate three scores based on single-sample Gene Set Enrichment Analysis (ssGSEA): stromal cell scores, immune cell scores, and ESTIMATE scores. CIBERSORT is a deconvolution technique that quantifies the proportions of distinct cell types by predicting the cellular composition of complicated tissues based on gene expression data. The link between risk score and tumor immune infiltration score was investigated using a total of seven methods.

The R package “pRRophetic” is used to predict clinical chemotherapeutic response using tumor gene expression levels. pRRophetic is capable of forecasting possibly sensitive medications that are ideal for patients based on data obtained from a vast amount of data regarding the response of various tumor cell lines to anticancer drugs. We looked for medications that may be more effective in treating high-risk patients by utilizing pRRophetic to predict the IC50 of certain anticancer agents.

Based on the results of multivariate cox regression, nomogram is a versatile approach of merging several risk factors into a single plot. We were able to visually forecast a patient’s survival probability using the nomogram created with the R packages ‘DynNom’.

American Type Culture Collection(ATCC) provided U87-MG and LN229 cells. Shanghai Institutes for Biological Sciences provided U251 and A172 cells (Shanghai, China). In four glioma cell line culture and in vitro investigations, Dulbecco’s Modified Eagle Medium (DMEM, gibco, CA, USA) with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution was utilized. Lonza provided normal human astrocytes (NHAs), which were cultivated in astrocyte growth medium containing rhEGF, insulin, ascorbic acid, GA-1000, L-glutamine, and 5% FBS. All of the cells were grown at 37°C with 5% CO2. Abcam provided antibodies against USP4, E-cadherin, and N-cadherin. Cell Signaling Technology provided the β-actin.

Ribobio (Guangzhou, China) developed and synthesized two distinct small interfering (si)RNAs against USP4. The target sequences of siRNA for USP4 were 5′- ACTGCAAAGTCGAGGTGTA-3′ (siUSP4-1) and 5′- GCAACACCTACGAGCAGTT-3′. (siUSP4-2). Genomeditech(Shanghai,China) designed and produced the USP4 over-expression plasmid. All transfections were carried out according to the manufacturer’s instructions using Lipofectamine 3000 (Invitrogen; Thermo Fisher Scientific, Inc.).

Total RNA was extracted from cell lines with TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s protocol. cDNA was synthesized with the PrimeScript RT Reagent Kit (Takara, Nanjing, China).qRT-PCR was implemented utilizing AceQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) on an ABI Stepone plus PCR system (Applied Biosystems, FosterCity, CA, USA). Primers used in this study were listed as follows: USP4(Forward): GCAGACACCATTGCAACCATC; USP4(Reverse): AACTGCTCGTAGGTGTTGCT. β-actin(Forward):GTCATTCCAAATATGAGATGCGT; β-actin(Reverse): GCATTACATAATTTACACGAAAGCA. Relative quantification was determined using the 2^-ΔΔCt method.

RIPA buffer with protease inhibitors (Roche) was used to lyse cellular proteins, and an equal amount of proteins was electrotransferred onto a polyvinylidene difluoride membrane (Millipore). The main and secondary antibodies were used to incubate the protein, which was then identified using enhanced chemiluminescence methods.

Cell counting kit-8 test (CCK-8) was used to assess the proliferation capabilities of GBM cells (87-MG, LN229, and A172) according to the manufacturer’s instructions. 96-well plates were used to seed the transfected cells. 10 μl of CCK-8 reagent was added to the test well at 24, 48, 72, and 96 hours after transfection and incubated for 2 hours at 37°C away from light. At a wavelength of 450 nm, the absorbance was measured.

U87-MG, LN229, and A172 cells were transfected and maintained in 6-well plates for about 12 days. The cells were then stained with 0.1 percent crystal violet for 30 minutes before being rinsed with PBS. If the colonies were larger than 1 mm in diameter, they were counted.

Cell migration and invasion were measured using transwell assays. In the upper chamber, 2×10⁴cells were cultivated in 200 μL media without serum, while in the lower chamber, 600 μL complete medium was supplied for the migration assay. According to the manufacturer’s protocols, additional Matrigel was employed for the invasion experiments (BD Biosciences, Bedford, MA, USA). Cells were fixed with 4 percent PFA and stained with 0.1 percent crystal violet solution after 24 hours of incubation at 37°C with 5% CO2.

Cells were grown in 6-well plates for 24 hours before being scratched with a sterile pipette tip (20 μL). Each wound was examined by inversion microscopy(Olympus, Japan) at 0 and 24 hours after rinsing the cells with PBS to remove cellular debris. To examine the cell migration capacity, the total wound area was analyzed using ImageJ software.

A total of 72 ubiquitin-related genes with prognostic value in glioma were identified by univariate Cox regression of 100 ubiquitin-related genes obtained from Genecards database. Through Lasso regression of the above 72 genes, we obtained a Risk Score formula consisting of 12genes: Risk Score=0.0137868801576863*UBE2D3+0.00798080087559059*UBE2D2+(-0.00918452991096172)*USP7+0.0065644500970919*GRN+0.00751875352337762*UBE2S+0.000192089285965977*UBB+(-0.00497853632444545)*UBE2G2+(-0.0968360751270892)*BTRC+0.00754166909940448*CUL1+0.104123805908217*USP4+0.0289487761185436*SIAH2+0.0271749540402458*UBE2Z ( Figures 2A, B ). Among them, UBE2D3, UBE2D2, GRN, UBE2S, UBB, CUL1, USP4, SIAH2, and UBE2Z were associated with poor prognosis of glioma(HR>1, P<0.001, Figure 2C ). USP7, UBE2G2, and BTRC were associated with better prognosis of glioma(HR<1, P<0.001, Figure 2C ). Using this formula, each patient can be calculated to obtain a risk score. Patients in different cohorts can be separated into high-risk and low-risk groups based on the median value. Figure 3A showed the survival status, score curve, and expression of model genes of the high-risk and low-risk groups of the training cohort. Figure 3B showsed the validation cohort analysis results. The dot plot of survival status, whether in the train cohort or the validation cohort, indicates that as the risk score grows, the survival time of patients gradually concentrates near the bottom, indicating a worse prognosis ( Figures 3A, B ). Survival analysis on the training cohort showed that the high-risk group has a significantly worse prognosis than the low-risk group ( Figure 3C ). The same result was found in the survival analysis of the validation cohort ( Figure 3D ). Then the results of subgroup survival analysis on train cohort showed that high risk score is associated with poor prognosis in different genders and age groups ( Figure 3E ) and the results were verified in the validation cohort ( Figure 3F ).

To determine the independent prognostic usefulness of risk score, we used univariate and multivariate Cox regressions. First, univariate Cox regression in the training cohort revealed that age and risk score are independent prognostic predictors of glioma ( Figure 4A ). Age and risk score were also revealed to be independent prognostic predictors of glioma in validation cohort analysis using univariate Cox regression ( Figure 4B ). After that, the multivariate Cox regression was run. Age and risk score were determined to be independent prognostic predictors of glioma in the training cohort using multivariate Cox regression ( Figure 4C ). Gender, age, and risk score were confirmed to be independent prognostic predictors of glioma in a second validation cohort analysis using multivariate Cox regression ( Figure 4D ). We then created ROC curves for this signature in both the train and validation cohorts to assess its accuracy. The area under the curve (AUC) of 1, 3, and 5 years was 0.869, 0.925, and 0.868, respectively, according to the ROC curve of the training cohort ( Figure 4E ). The AUC of 1, 3, and 5 years was 0.854, 0.867, and 0.796, respectively, according to ROC curves for the validation cohort, demonstrating that the signature can accurately determine the prognosis of patients with glioma ( Figure 4F ).

Tumor formation and progression are influenced by the immune microenvironment. Understanding the effects of the immunological microenvironment on tumor prognosis and treatment is beneficial. As a result, we looked at how immune infiltration differed between the high-risk and low-risk groups. To begin, we used multiple algorithms to create an immunological infiltraion heat map for high-risk and low-risk groups, with red representing high invasion levels and blue representing low invasion levels ( Figure 5A ). Following that, we conducted a correlation analysis between immune cells and risk ratings, finding that many immune cells were substantially connected with risk scores ( Figures 5B–I ).

Tumor immunotherapy is a promising treatment option. Immune checkpoint-related gene expression and microsatellite instability are crucial indications for evaluating immunotherapy’s effectiveness. Between high-risk and low-risk groups, we looked at differences in immune checkpoint-related genes and microsatellite instability. The high-risk group had higher levels of expression of immune checkpoint genes, according to the findings ( Figure 6A ). The high-risk group had less microsatellite instability ( Figure 6B ). Microsatellite instability decreased as the risk score grew, according to correlation analyses ( Figure 6C ).

We used the Seurat package to process the single-cell transcriptome data. Raw data from GSE162631 were downloaded. A total of 4 GBM samples were downloaded for further analysis. The data with a mitochondrial RNA percentage larger than 0.10 were filtered, and we eventually acquired 51,449 cells that meet the standard. PCA reduction plot showed no significant differences in cell cycles ( Figure 7A ). Meanwhile, we selected the top 3000 variable features, which were labeled in red( Figure 7B ). And the top 10 variable features were labeled. In Figures 7C , the principle component analysis showed the distribution of the samples, and the results showed no significant batch effects. After dimension reduction and the tSNE clustering, the immune cells were identified using their feature genes ( Figure 8A ). It could be clearly recognized that the cells were divided into approximately 6 categories of cell types, namely endothelial cells(EC), neutrophils, T cells/B cells, mural cells, tumor-associated macrophages(TAMs), didenric cells and microglias ( Figure 8B ). The genes involved in the signature were identified in the single-cell tumor microenvironment atlas( Figure 8C ). In a fleeting glimpse, we could see the expression of the gene UBE2D3, UBE2D2, GRN, UBB was ubiquitous in immune-microenvironment in the GBM patients, and UBE2G2, CUL1, and USP4 showed moderate expression in immune cells. Besides, almost all ubiquitination associated genes were invariably expressed in TAMs and Dendritic Cells, indicating the strong correlation of ubiquitination within those immune cells.

According to the findings, high-risk patients have a worse prognosis. Therefore, in order to conduct precise intervention in high-risk patients, we performed drug sensitivity analysis to identify drugs that might be effective. Results showed Salubrinal, Rapamycin, Paclitaxel, Midostaurin, JW.7.52.1, Dasatinib, Cyclopamine, Bryostatin.1, Bexarotene, Bortezomib, MG.13 2., and Parthenolide had a lower semi-inhibitory concentration (IC50) in the high-risk group, meaning that the high-risk group was more sensitive to these drugs ( Figure 9 ).

By merging the clinical parameters of glioma patients, a nomogram was created to evaluate the prognosis of glioma patients at 1, 3, and 5 years ( Figure 10A ).

The role of USP4 in glioblastoma was confirmed since its HR value was the highest in the signature. First, a survival analysis using the GEPIA database revealed that increased USP4 expression in glioma patients was linked to poor outcomes ( Figure 10B ). USP4 was found to be linked to a variety of immune cells in immunological investigation ( Figure 10C ). GSVA analysis of USP4 and ubiquitination related pathways is presented in Supplemental Figure S1 . In vitro tests were then carried out to confirm USP4’s function. To begin, qRT-PCR revealed that USP4 expression was up-regulated in all four glioma cell lines when compared to normal control NHAs cell lines ( Figure 11A ), with the highest expression in U87-MG and LN229 cell lines, so gene knockdown was performed in these two cell lines ( Figure 11B , *P<0.05, **P<0.01). In the U87-MG and LN229 cell lines, both siRNAs drastically reduced USP4 expression ( Figure 11C ). USP4 knockdown significantly reduced the activity of U87-MG and LN229 cell lines in the CCK-8 experiment ( Figure 11D , **P<0.01). The ability of the U87-MG and LN229 cell lines to form colonies was dramatically reduced following USP4 knockdown ( Figure 11E , **P<0.01). After knocking out USP4, the migratory and invasion capacities of the U87-MG and LN229 cell lines were dramatically reduced ( Figure 11F , **P<0.01). The ability of U87-MG and LN229 cell lines to heal was dramatically reduced following USP4 knockdown ( Figure 11G , **P<0.01). Following that, plasmid overexpressed USP4 in the A172 cell line ( Figure 12A ). The colony forming ability of the A172 cell line was greatly improved after overexpression of USP4 ( Figure 12B , **P<0.01). The viability of the A172 cell line was dramatically improved following USP4 overexpression in the CCK-8 experiment ( Figure 12C , **P<0.01). In a transwell experiment, USP4 overexpression greatly improved the A172 cell line’s migration and invasion capacity ( Figure 12D , **P<0.01). In wound healing assays, USP4 overexpression greatly improved the migratory ability of the A172 cell line ( Figure 12E , **P<0.01). The knockdown and overexpression efficiencies of USP4 were confirmed by Western-blotting, and the link between USP4 and EMT-related proteins N-cadherin and E-cadherin was investigated. A statistically significant association was established between USP4 and the EMT proteins E-cadherin and N-cadherin ( Figure 12F ). N-cadherin expression was dramatically reduced in the siUSP4-1 and SiUSP4-2 groups when the USP4 gene was knocked down in the U87-MG cell line, whereas E-cadherin expression was significantly raised. N-cadherin expression was greatly reduced when the USP4 gene was knocked down in the LN229 cell line, whereas E-cadherin expression was significantly increased. After USP4 was overexpressed in A172, N-cadherin expression was drastically increased, but E-cadherin expression was significantly decreased.