The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) provided us with RNA sequencing, corresponding clinicopathology, and genomic mutation data for glioma patients. 1152 normal brain tissue sample information was collected from the genotype-tissue expression (GTEx) database (https://www.gtexportal.org/home/index.html). Similarly, the RNA-sequencing expression data and matched clinical information were downloaded from the GSE16011 and GSE109857 datasets of the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/) and from the CGGA_301, CGGA_325, and CGGA_693 datasets of the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn/). A total of 3310 genes associated with palmitoylation were retrieved from GeneCards database (www.genecards.org). The raw count data were normalized using the “limma” package in the R software (Version 4.2.0) (15). The flowchart of this study design is shown in Figure 1.

The “limma” R package was utilized to identify differentially expressed genes (DEGs), applying the threshold of “|log2FC | > 2 and FDR < 0.05” (16). We identified the overlapping DEGs associated with palmitoylation in datasets GSE16011 and GSE109857, presenting them through a Venn diagram. GeneMANIA (http://genemania.org/) (17) was used to build the gene-gene interaction network. Functional enrichment analysis, including Kyoto encyclopedia of genes and genomes (KEGG) and Gene Ontology (GO), were performed using the “ClusterProfiler” package (18). Gene set enrichment analysis (GSEA) was also used to evaluate the biological activity and pathways in high-risk and low-risk glioma patients. The gene sets corresponding to various pathways were obtained from the Molecular Signatures Database (MSigDB). To quantify the enrichment level and statistical reliability, the normalized enrichment score (NES) and false discovery rate (FDR) were utilized, respectively (19). Simultaneously, Metascape database (https://metascape.org/) integrates functional enrichment, interactome analysis, gene annotation, and membership search to offer experimental biologists a comprehensive means of analyzing and interpreting gene expression profiles. Additionally, NetworkAnalyst 3.0 (https://www.networkanalyst.ca/) was utilized to construct the coregulatory networks of transcription factor (TF)-microRNA (miRNA) for PRGs.

Utilizing the “ConsensusClusterPlus” package in R (20), unsupervised hierarchical clustering was performed on glioma patients, classifying them into the most appropriate number of clusters based on the palmitoylation-associated intersection genes. The “PCA” package was utilized to conduct principal component analysis (PCA) across various subgroups.

The stromal score, immune score, and estimate score were calculated by the “ESTIMATE” R package for samples in the TCGA-GBMLGG, CGGA_301, CGGA_325, and CGGA_693 cohorts. Tumor-infiltrating immune cells (TIICs) were assessed in the TCGA-GBMLGG, CGGA_301, CGGA_325, and CGGA_693 cohorts using the CIBERSORT online tool (21). The two-sample Wilcoxon test was employed to compare immune infiltration and functions among various groups. We further compared the gene expression differences of common immune checkpoint inhibitors (ICIs), major histocompatibility complex (MHC), and T-cell stimulators between different clusters. A method for assessing the effectiveness of immune checkpoint blockade (ICB) is the Tumor Immune Dysfunction and Exclusion (TIDE) score, which is available on the TIDE website (http://tide.dfci.harvard.edu/) (22). A range of predictive indicators for the ICIs response, encompassing the TIDE score, microsatellite instability (MSI) score, T-cell exclusion score, and T-cell dysfunction score, were utilized to assess the correlation between various clusters and the efficacy of immunotherapy, as well as the relationship between different PRRS groups and the outcomes of immunotherapy.

The waterfall plots were utilized, produced by the “maftools” R package (23), to evaluate the somatic point mutation counts in each sample within the TCGA-GBMLGG cohort, allowing for comparisons of genetic variations among distinct clusters. Kaplan-Meier survival analysis was performed incorporating both TMB and subgroup classification, categorizing patients into four prognostic groups to assess the survival disparities between those with high versus low TMB values within distinct subgroups.

Based on the intersection of PRGs and DEGs, univariate Cox regression analysis was used to evaluate the prognostic significance of these genes, and genes with P values less than 0.05 were selected for further analysis. To further refine the selection of candidate prognostic-related genes, we employed the LASSO Cox regression algorithm via the R package “glmnet” to discern variations in the regression coefficients of prognostic genes (24). 9 PRGs with nonzero coefficients were selected through 10-fold cross-validation. We designated the TCGA-GBMLGG cohort as the training set, with CGGA_301, CGGA_325, and CGGA_693 serving as independent test sets. The PRRS formula is as follows: Risk score = Σ (Coefi × Exp) Coefi indicates the risk coefficient, and Exp indicates the expression level. We calculated the risk score (RS) for each patient and divided them into high-risk and low-risk groups based on the median RS. A risk curve was produced to explore the differences in survival status between patients in different risk groups. Moreover, based on both training and test datasets, we performed Kaplan-Meier survival analysis on risk scores and generated time-dependent receiver-operating characteristic curves to appraise the predictive ability of the risk groups identified above. The R package “pheatmap” was used to create heatmaps that showed the variation in PRGs expression as well as the distribution pattern across clinicopathological traits and risk score groups.

To assess whether the risk score serves as an independent prognostic factor influencing the survival of glioma patients, we performed univariate and multivariate analysis on the TCGA-GBMLGG dataset. Nomograms are extensively utilized in cancer prognosis research, serving as statistical predictive models that integrate multiple risk factors to estimate the individual survival probabilities (25). Utilizing the TCGA-GBMLGG dataset, we constructed a prognostic nomogram model to predict the 1-, 3-, and 5-year overall survival (OS) probabilities employing the “rms” package in R. Based on the training and test datasets, calibration curves were plotted for graphical assessment. The concordance index (C-index) and the receiver operating characteristic (ROC) analysis were employed to assess the predictive accuracy of the nomogram, while the decision curve analysis (DCA) was utilized to evaluate the clinical utility of the nomogram.

The protein expression levels corresponding to PRGs screened by the LASSO Cox regression algorithm in glioma were explored through the UALCAN data analysis portal (https://ualcan.path.uab.edu/). In the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/), the protein expression differences of PRGs in glioma of different malignant degrees and normal tissues were compared based on immunohistochemical images.

The TIMER database (http://timer.cistrome.org) (26) was utilized to evaluate the correlations between the expression levels of PRGs and the infiltration rates of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells. The infiltration enrichment of 24 common immune cells was presented using the ssGSEA method from the R package “GSVA” (27). In addition, the relationships between PRGs and immunomodulator expression levels in GBM and LGG were examined using the TISIDB database (http://cis.hku.hk/TISIDB/index.php). Ultimately, the expression status of PRGs across various cell types within the tumor microenvironment (TME) was investigated based on RNA sequencing datasets from 10x Genomics (GSE131928 and GSE163108) and Smart-seq2 (GSE89567 and GSE84465), retrieved through the Tumor Immune Single-cell Hub 2 (TISCH2) database (http://tisch.compbio.cn/) (28).

The cancerSEA database was employed to ascertain the mean correlation between PRGs and various functional states in glioma, encompassing angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, epithelial-mesenchymal transition (EMT), hypoxia, inflammation, invasion, metastasis, proliferation, quiescence, and stemness.

The mutational landscape of PRGs in glioma was investigated utilizing the cBioPortal database (http://cbioportal.org) (29). We conducted an in-depth analysis of the heterozygous and homozygous copy number variations (CNV) of PRGs, as well as the correlation between the CNV and mRNA expression levels of the respective genes via the GSCALite platform (http://bioinfo.life.hust.edu.cn/web/GSCALite/). Furthermore, the platform was employed to investigate the correlation between PRGs methylation and their expression profiles in glioma, as well as to assess the influence of varying methylation states on survival outcomes. MethSurv (https://biit.cs.ut.ee/methsurv/) is an online program that performs variable survival analysis using DNA methylation data. Hypomethylation at certain CpG sites in glioma was explored for its prognostic significance.

Based on the GSCALite website, we studied the correlations between PRGs expression levels and drug response. Considering the absence of biomarkers capable of precisely forecasting chemotherapeutic drug sensitivity in glioma patients, we utilized the “pRRophetic” and “ggplot2” packages to perform drug susceptibility analysis, comparing the half-maximal inhibitory concentration (IC50) values of various chemotherapeutic agents against glioma between high- and low-risk groups via the Wilcoxon signed-rank test.

Utilizing Autodock (https://ccsb.scripps.edu/mgltools/downloads/) for molecular docking analysis allowed for the investigation of the interactions between the identified PRGs and potential therapeutic compounds. The small molecular drugs and protein structures corresponding to the PRGs were sourced from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/) and PDB database (https://www.rcsb.org/), respectively. Ultimately, biological macromolecules and small molecular drugs were automatically docked, adhering to the standardized docking procedure. The results were visualized using the PyMol. A 100 ns molecular dynamics (MD) simulation was conducted using Gromacs software to further substantiate the plausibility and dependability of the docking outcomes (30). In order to evaluate the binding stability of the receptor-ligand complex, MD simulation was carried out after achieving equilibrium in a human-like environment using the TIP3P water model, which represents a 0.145 mol/L neutral sodium chloride solution. The physical conditions were set to constant pressure (101 kPa), constant temperature (310 K), and periodic boundary conditions (31, 32). Conformations were recorded and computed every 10 ps. The root mean square deviation (RMSD) and root mean square fluctuation (RMSF) of the MD simulation results analysis and visualization were performed using a GROMACS inbuilt program with the (Equation 1) RMSD and (Equation 2) RMSF formula as follows.

In the above equations, Rt−Rref denotes the displacement of the t-th atom at a specific frame from its position in the reference conformation (ref) (i.e., positional offset), N stands for the total number of atoms, and T indicates each sampling time point.

The LN229 and U251 cell lines were obtained from BeNa Culture Collection (BNCC, Henan, China). Cells were cultured in high-glucose DMEM (HyClone, USA) supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and 1% penicillin-streptomycin (PS, Gibco, USA), and the cultures were maintained at 37 °C in a humidified 5% CO₂ atmosphere. The small-molecule compound was dissolved in dimethyl sulfoxide (DMSO) to prepare a 10 mM stock solution. This stock solution was subsequently subjected to serial twofold dilutions using cell culture medium, thereby establishing a series of concentration gradients. In none of the groups did the final concentration of DMSO surpass 0.1%.

Cell viability was evaluated using the Cell Counting Kit-8 (CCK-8) assay (EpiZyme, Cat. No. CX001M). Cells were seeded into 96-well microplates at a density of 5 × 10³ cells per well and incubated overnight to allow for complete cell adherence. The cells were then cultured in culture medium with varying drug concentrations for 48 hours. Following treatment, 10 µL of CCK-8 reagent was added to each well, followed by incubation at 37 °C for 2 hours. Using a microplate reader (Bio-Rad), the absorbance was measured at 450 nm, cell viability was calculated relative to the control group, and the IC50 value was determined with GraphPad Prism 9.0.

Cells were inoculated into 6-well culture plates and incubated until they attained 95% confluence. Subsequently, a consistent linear scratch was generated using a sterile 200 µL pipette tip. After washing with phosphate-buffered saline (PBS), the cells were incubated in a serum-free culture medium containing a drug concentration equivalent to 0.1 times the IC50 value. Wound area images were captured at 0, 24, and 48 h under a microscope. The scratch area was measured using ImageJ software.

Cell apoptosis was assessed using the Annexin V-FITC/PI Apoptosis Detection Kit (Beyotime, Cat. No. C1062S-2), strictly adhering to the guidelines provided by the manufacturer. After being cultured for 24 hours in a medium with a drug concentration equivalent to 0.1 times the IC50 value, the cells were harvested and then washed with PBS. The cells were resuspended in binding buffer, followed by incubation with 5 µL of Annexin V-FITC and 10 µL of propidium iodide (PI) staining solution at room temperature for 15 minutes in the dark. The CytoFLEX flow cytometer was used to analyze apoptosis as soon as staining was completed.

R software (version 4.1.0, https://www.r-project.org/) was used to perform all statistical analyses in this study. The independent Student’s t-test was used to compare continuous variables between two groups, the Wilcoxon rank-sum test was used to compare non-normally distributed variables, and the chi-squared test was used to compare data on categorical variables between two groups. Survival curves were generated employing the Kaplan–Meier method and assessed for comparison through the log-rank test. P value< 0.05 was considered statistically significant.

The datasets GSE16011 and GSE109857 have been collected and utilized in order to screen for DEGs between glioma and normal brain tissue. 225 genes and 534 genes are significantly upregulated and downregulated in the glioma tissues, respectively, according to the analytical results of the GSE16011 dataset (Figure 2A). 76 genes and 332 genes are significantly upregulated and downregulated in the glioma tissues, respectively, according to the analytical results of the GSE109857 dataset (Figure 2B). Then, the intersection of DEGs from the two datasets and a set of 3310 PRGs was determined, ultimately identifying 46 common DEGs associated with palmitoylation (Figure 2C). Figure 2D displays the heatmap of 46 intersecting genes. Furthermore, the 46 DEGs were uploaded to the GeneMANIA database to establish a gene co-expression network, which demonstrated their associated functionalities. As depicted in Figure 2E, these genes exhibit close correlations across diverse biological functions, encompassing signal release from synapse, GABA receptor activity, neurotransmitter transport, regulation of neurotransmitter levels, ligand-gated anion channel activity, and presynapse. Subsequently, GO and KEGG enrichment analyses were conducted to investigate the underlying biological roles of the 46 DEGs. Based on KEGG pathway enrichment analysis, the results showed that these genes were primarily involved in morphine addiction, insulin secretion, synaptic vesicle cycle, retrograde endocannabinoid signaling, and GABAergic synapse (Figure 2F). Furthermore, GO analysis demonstrated that these genes were mainly connected with signal release from synapse, neurotransmitter secretion, and synaptic vesicle cycle (Figure 2G).

Given the heterogeneity inherent within tumors, TCGA-GBMLGG patients were divided into two subgroups with different molecular and clinical characteristics by unsupervised consensus analysis, based on the expression levels of 46 palmitoylation-related DEGs, including 502 samples in cluster 1 and 164 samples in cluster 2 (Figures 3A, B). PCA revealed that the 46 genes linked to palmitoylation could clearly discriminate between the two clusters (Figure 3C). Following that, a heatmap using a consensus clustering solution (k = 2) displays the two clusters (Figure 3D). Figure 3E shows the distribution of clinical pathological characteristics in clusters 1 and 2, as well as the transcriptomic characteristics of PRGs differentially expressed in the two subtypes. Except for gender, there were statistically significant differences in survival rates, age, WHO grade, isocitrate dehydrogenase (IDH) mutation status, 1p/19q codeletion status, and O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation status between the two clusters (Figures 3G-M). In addition, we found that the OS of patients in cluster 2 was significantly better than that of patients in cluster 1 (Figure 3F). These suggest that PRGs may exert an influence on the progression of glioma via certain underlying mechanisms. 1223 genes were screened out by DEGs analysis between the two clusters, of which 301 genes were upregulated and 922 genes were downregulated (Supplementary Figure 1A). The heatmap was used to display the expression of DEGs based on hierarchical clustering analysis (Supplementary Figure 1B). Then, GO enrichment and KEGG pathway analysis were performed on the up- and down-regulated DEGs, respectively (Supplementary Figures 1C–F).

Considering the close relationship between palmitoylation and the tumor microenvironment, the differences in immune cells and immune infiltration between the two clusters were assessed utilizing the ESTIMATE and CIBERSORT algorithms (33). We found that the stromal score, immune score, and estimate score in cluster 1 were significantly higher than those in cluster 2 (Figure 4A). Cluster 1 was distinguished by a prominent infiltration of memory B cells, activated CD4+ memory T cells, regulatory T cells (Tregs), resting natural killer (NK) cells, M0 macrophages, M1 macrophages, M2 macrophages, resting dendritic cells (DCs), activated dendritic cells, and resting mast cells. In contrast, Cluster 2 was characterized by a significant infiltration of naive B cells, plasma cells, naive CD4+ T cells, resting CD4+ memory T cells, follicular helper T cells (Tfh), gamma delta T cells (γδ T cells), activated NK cells, monocytes, activated mast cells, and neutrophils (Figure 4B). These results indicated that the aforementioned two clusters exhibited distinctly disparate infiltration characteristics within the TME. Then, we studied the relationship between the two clusters and ICIs, MHC, as well as T-cell stimulators. Immune checkpoint genes (HAVCR2, CTLA4, PDCD1, CD28, CD80, LAG3, CD96, CD86, and PDCD1LG2), MHC, and T-cell stimulators (CD2, CD226, CD27, CD28, CD40LG, ICOS, TNFRSF14, TNFRSF4, TNFRSF8, and TNFRSF9) were highly expressed in cluster 1 (Figures 4C–E). The TIDE algorithm was employed to assess the risk of tumor immune escape. The findings revealed that, when compared with cluster 2, cluster 1 demonstrated a suboptimal response to immunotherapy and a higher TIDE score, suggesting an increased probability of immune escape (Figures 4F, G). In addition, we found that the T-cell exclusion score was higher in cluster 1, while the MSI and T-cell dysfunction scores were higher in cluster 2 (Figures 4H–J).

TMB exhibits a close correlation with the immune cell infiltration and therapeutic efficacy of immunotherapy (34, 35). The somatic mutation data for patients with TCGA-GBMLGG was acquired, and the TMB scores were subsequently calculated. The results showed that the TMB values in cluster 1 were slightly higher than those in cluster 2 (Figure 4K). The mutation differences of the top 20 genes across various glioma clusters were displayed using waterfall plots (Figures 4L, M). In glioma patients of clusters 1 and 2, missense mutations were the most prevalent type of mutation (Supplementary Figures 2Aa, Ba), single nucleotide polymorphism (SNP) held an absolute position in contrast to insertion (INS) or deletion (DEL) (Supplementary Figures 2Ab, Bb), and the most common mutation type identified was C>T (Supplementary Figures 2Ac, Bc). Supplementary Figures 2Ad, Bd displayed the number of mutations per sample. A type of mutation was represented by the box diagram for each color in Supplementary Figures 2Ae, Be. The horizontal histograms illustrated genes with a higher mutation frequency, such as in cluster 1 (IDH1 (54%), TP53 (45%), ATRX (27%), TTN (16%), CIC (14%)) and cluster 2 (IDH1 (79%), TP53 (35%), ATRX (26%), CIC (23%), NOTCH1 (13%)) (Supplementary Figures 2Ae, Be). Given the potential correlation between palmitoylation-related clusters and TMB, we conducted a stratified analysis and found that the integration of clusters with TMB scores offers a more accurate prediction of the prognosis for glioma patients (Figure 4N). Furthermore, given the pivotal roles that N6-methyladenosine (m⁶A) methylation and ferroptosis play within the immune microenvironment, we conducted an investigation into the expression profiles of m⁶A-associated genes and ferroptosis-related genes across distinct clusters. The results showed that most m⁶A-related genes and ferroptosis-related genes were highly expressed in cluster 1 (Supplementary Figures 2C, D).

Based on the TCGA-GBMLGG cohort, a univariate Cox regression analysis was conducted on 46 DEGs, resulting in the identification of 35 PRGs that were significantly associated with prognosis (Figure 5A). The LASSO Cox regression analysis was employed to filter the best prognostic indicators from 35 PRGs. After incorporating variables into the LASSO Cox regression model corresponding to the minimum lambda value, 9 PRGs were selected for the construction of the prognostic PRRS model (Figures 5B, C). The PRRS was computed using the following formula: risk score = (0.237684096 * APOC1) + (0.296416314 * FXYD1) + (0.033378972 * ZCCHC12) + (0.031624620 * F2R) + (0.043403341 * PTBP1) + (0.308436072 * NCAPG) + (-0.317672694 * BMP2) + (-0.108996577 * PDE2A) + (0.342006851 * IFI30). The risk curve analysis conducted within the TCGA-GBMLGG, CGGA_325, and CGGA_301 cohorts demonstrated a correlation between elevated risk scores and an increased risk of mortality (Figures 5D–F). The outcomes of the Kaplan-Meier survival analysis have confirmed that the survival prospects for patients in the high-risk group are markedly inferior to those in the low-risk group. The p-values obtained from both the training cohort (TCGA-GBMLGG) and the validation cohorts (CGGA_325 and CGGA_301) were all statistically significant, being less than 0.001 (Figures 5G–I). A more in-depth time-dependent ROC curve analysis was subsequently performed to ascertain the robustness of the model. The results demonstrated its robust predictive capability for 1-year, 3-year, and 5-year survival rates, which was further corroborated in the validation cohort (Figures 5J–L). These findings underscore the potential efficacy of the model in clinical risk assessment scenarios.

We observed a significant correlation between the risk score and factors including age, WHO grade, IDH mutation status, 1p/19q codeletion status, and MGMT promoter methylation status within the TCGA-GBMLGG cohort (Figure 6A). Specifically, glioma patients older than 60 years have a higher risk score. A higher WHO grade is associated with an increased RS among patients. Individuals with wild-type IDH demonstrate higher risk scores. Patients presenting with 1p/19q non-coding deletions show elevated risk scores. Furthermore, patients with non-methylated MGMT promoters exhibit higher risk scores (Figure 6B). Similarly, these observations were corroborated in the external datasets, namely CGGA_325, CGGA_301, and CGGA_693 (Figures 6C–E).

Due to the significant correlation between PRRS and high malignancy in glioma, it has been demonstrated that PRRS was an independent predictor of survival outcomes through univariate and multivariate Cox regression analysis based on the TCGA-GBMLGG cohort (Figures 7A, B). To evaluate the prognostic utility of combining PRRS with clinicopathological parameters, we constructed a nomogram based on WHO grade, age, 1p/19q codeletion status, MGMT promoter methylation status, and PRRS to predict the 1-, 3-, and 5-year survival rates of glioma patients (Figure 7C). The results of the proportional hazards (PH) (36) assumption testing for the multivariate Cox regression model showed that both the nomogram and the included variables conformed to the PH assumption (Supplementary Table 1). The C-index of the nomogram was 0.884 in the TCGA-GBMLGG cohort, 0.749 in the CGGA_325 cohort, 0.761 in the CGGA_301 cohort, and 0.731 in the CGGA_693 cohort, respectively. The calibration curves derived from multiple cohorts demonstrated excellent concordance between the survival probabilities predicted by the nomogram and those observed in reality (Figures 7D–G). The area under the ROC curve (AUC) values were greater than 0.75 at multiple time points, indicating that the model had satisfactory discriminative capacity (Figure 7H). DCA further demonstrated that the overall net benefit derived from the nomogram exceeded that of any individual clinical characteristic (Figure 7I). The above results suggest that the nomogram we established has good prognostic value for glioma patients and may have practical significance for clinical decision-making.

Kaplan-Meier survival analysis revealed a statistically significant correlation between the 9 genes incorporated in the PRRS model and overall survival in the TCGA-GBMLGG cohort, highlighting their crucial prognostic significance (Supplementary Figure 3A). The differential expression analysis results showed that the mRNA expression levels of FXYD1, ZCCHC12, and PDE2A were highly expressed in normal brain tissues, while the other six genes were significantly upregulated in glioma samples (Supplementary Figure 3B). We further assessed the protein expression levels corresponding to nine PRGs utilizing the UALCAN and HPA databases. Regrettably, data pertaining to protein expression and immunohistochemical outcomes for certain PRGs were deficient in these two databases. The results of the UALCAN database showed that in GBM, the protein expression of APOC1, F2R, PTBP1, NCAPG, and IFI30 significantly increased, while the protein expression of PDE2A significantly decreased (Supplementary Figure 3C). The immunohistochemical findings within the HPA database further elucidated the differences in protein expression levels of PRGs between normal brain tissue and LGG and HGG samples (Supplementary Figure 4). These findings lend credence to the possible involvement of these PRGs in the development of glioma as well as their value in forecasting patient outcomes.

For a more in-depth analysis of the differences in gene functions and pathways involved among the subgroups classified by PRRS, 605 and 1105 DEGs were determined between the low-risk and high-risk groups in the GSE16011 dataset and GSE109857 dataset, respectively, with 297 genes being common to both datasets (Figures 8A–C). Following that, KEGG pathway analysis and GO enrichment analysis were carried out on these DEGs in the TCGA-GBMLGG dataset. The results of enrichment analysis indicated that these DEGs were predominantly enriched in pathways such as external encapsulating structure organization, extracellular structure organization, and extracellular matrix organization within the biological process category, as well as in pathways associated with complement and coagulation cascades and extracellular matrix (ECM)-receptor interaction within the KEGG category (Figures 8D, E). The Metascape enrichment results were colored by cluster ID, where nodes sharing the same cluster ID were usually close to each other. The enrichment results encompassed NABA CORE MATRISOME, PID INTEGRIN1 PATHWAY, complement and coagulation cascades, cellular response to growth factor stimulus, and inflammatory response (Figure 8F). Additional GSEA analysis was undertaken within the GSE16011 and GSE109857 datasets to supplement and verify the functional annotations derived from KEGG and GO pathways. The most enriched pathways in the GSE16011 dataset were the integrin1 pathway, TYROBP causal network in microglia, complement system, complement and coagulation cascades. The most enriched pathways in the GSE109857 dataset were the cell cycle checkpoints, PLK1 pathway, systemic lupus erythematosus, and TYROBP causal network in microglia (Figures 8G, H). Additionally, the coregulatory network of TF-miRNA for the 9 PRGs was established utilizing NetworkAnalyst (Figure 8I).

We evaluated the differences in immune infiltration between different risk groups in the training cohort (TCGA-GBMLGG) and the validation cohorts (CGGA_325, CGGA_301, and CGGA_693), respectively. The findings revealed that both the immune and stromal scores were higher in the high-risk group compared to those in the low-risk group (Figures 9A–D). High immune score corresponded to low survival rate (Figure 9E). Then, we compared the differences in TIICs, ICIs, MHC, and T-cell stimulators between different risk groups in multiple cohorts. The results of the training and validation cohorts showed that the abundance of M2 macrophages in the high-risk group was significantly higher than that in the low-risk group (Figures 9G–J). The mRNA expression levels of immune checkpoint genes, including CD274, HAVCR2, CTLA4, PDCD1, TIGIT, CD80, LAG3, CD96, and PDCD1LG2, were notably elevated in the high-risk group relative to the low-risk group, whereas the mRNA level of RGMB was conspicuously decreased (Supplementary Figures 5A-D). Moreover, with the exception of HLA-DOB and HLA-DQA1, the mRNA levels of MHC were significantly higher in the high-risk group compared to the low-risk group (Supplementary Figures 5E-H). Lastly, the mRNA expression levels of T-cell stimulators such as CD2, CD226, CD27, ICOS, TNFRSF14, TNFRSF18, TNFRSF4, and TNFRSF8 were notably increased in the high-risk group when compared to the low-risk group (Supplementary Figures 5I-L). In addition, 9 PRGs were correlated with the infiltration of various TIICs (Figure 9F).

Accumulating evidence indicates that the infiltration of immune cells within the TME plays a pivotal role in tumor progression (37). However, it remains unclear whether elevated expression of PRGs influences the recruitment of immune cells in glioma. An analysis conducted utilizing the TIMER database revealed a positive correlation between the expression levels of APOC1 and IFI30 in GBM and the abundance of infiltrating B cells, CD8+ T cells, macrophages, and neutrophils. Similarly, in LGG, the expression levels of APOC1, PTBP1, NCAPG, and IFI30 exhibited a positive correlation with the infiltration of CD4+ T cells, macrophages, neutrophils, and dendritic cells (Supplementary Figure 6). Subsequently, we constructed lollipop plots to illustrate the associations between 9 PRGs and the infiltration levels of 24 immune cell types in glioma. The expression levels of APOC1 and IFI30 were positively correlated with the infiltration rates of most immune cells (Supplementary Figure 7). TIICs residing within the microenvironment are capable of secreting a diverse array of cytokines, thereby making them essential to the progression of glioma. We visualized the relationship between the expression levels of 9 PRGs and various cytokines associated with infiltrating immune cells in the form of heatmaps. The findings demonstrated that the majority of immunostimulators, immunoinhibitors, MHC molecules, chemokines, and chemokine receptors were positively related to the expression levels of APOC1 and IFI30 in GBM and LGG (Supplementary Figures 8A-E). The expression of important genes is changed by DNA methylation, which impacts the development of glioma (38). Therefore, we investigated the DNA methylation associated with PRGs, and the results showed that the methylation of APOC1 and IFI30 was negatively correlated with most cytokines (Supplementary Figure 8F-J). These results unveil new glioma treatment targets and prognostic indicators.

Given the prognostic significance and abnormal distribution of PRGs in glioma tissues, we conducted an investigation into the cell types in which these genes are enriched based on single-cell RNA sequencing (scRNA-seq) data. Analysis of the scRNA-seq data in the GSE131928 dataset revealed the identification of 27 cell clusters and 8 cell types in glioma tissues. We observed significant enrichment of APOC1 and IFI30 in monocytes/macrophages (Figure 10A). 19 cell clusters and 6 cell types were identified in glioma tissues by scRNA-seq data analysis of the GSE163108 dataset. It was observed that APOC1 and IFI30 were markedly enriched within monocytes/macrophages, notably in the C17 cluster. Similarly, NCAPG exhibited significant enrichment in Tprolif cells, predominantly in the C8 cluster. Furthermore, notable enrichments of F2R and PTBP1 were identified within conventional CD4+ T cells, CD8+ T cells, Tprolif cells, and regulatory T cells (Figure 10B). To substantiate our discoveries, we conducted an analysis of the scRNA-seq datasets GSE89567 and GSE84465. Consistently, our results demonstrated a significant enrichment of APOC1 and IFI30 within monocytes/macrophages (Supplementary Figures 9A, B).

To further explore the functional significance of PRGs in glioma, single-cell analysis was conducted using the CancerSEA database, revealing a significant correlation between PRGs expression and pivotal cellular functions at the single-cell level. The research results found that in HGG, PTBP1 was positively correlated with the cell cycle, DNA damage, DNA repair, proliferation, and stemness, while negatively correlated with angiogenesis, hypoxia, inflammation, and quiescence. FXYD1, ZCCHC12, and PDE2A exhibited a negative correlation with DNA damage, DNA repair, epithelial-mesenchymal transition (EMT), hypoxia, and invasion in GBM. NCAPG was positively correlated with cell cycle, DNA damage, DNA repair, EMT, invasion, metastasis, and proliferation in glioma (Supplementary Figure 10). The above results suggest that PRGs may be involved in regulating multiple biological functions of glioma cells.

The genomic variations of 9 PRGs in glioma were analyzed using the cBioPortal database. In the database, 8.6% (95/1099) of glioma patients were found to carry PRGs gene variations, with PTBP1 accounting for the highest proportion at 2.8% (Supplementary Figure 11A). Mutation and amplification constituted the most common genetic alterations observed in PRGs among patients with GBM and LGG (Supplementary Figure 11B). Furthermore, we investigated the correlation between the mRNA expression levels and CNV of PRGs (Supplementary Figure 11C). Next, we investigated the CNV status of these genes in GBM and LGG. PRGs showed a high sensitivity to mutation, with APOC1, NCAPG, FXYD1, and F2R mainly exhibiting heterozygous deletion in GBM, whereas PDE2A, BMP2, IFI30, and PTBP1 exhibited heterozygous amplification as the main CNV type. In contrast, ZCCHC12 and PDE2A primarily demonstrated heterozygous deletion in LGG, while APOC1, FXYD1, BMP2, IFI30, and PTBP1 showed heterozygous amplification as the main CNV pattern (Figure 11A). The susceptibility of these genes to heterozygous deletion and amplification was validated by the CNV analysis in Figure 11B. Figure 11C presented a notable correlation between CNV and the expression of PTBP1, PDE2A, APOC1, and FXYD1 in LGG, alongside a significant association between CNV and the expression of PTBP1 specifically in GBM.

We found that the expression levels of PRGs were negatively correlated with methylation levels. Notably, APOC1, ZCCHC12, NCAPG, IFI30, and FXYD1 showed significant statistical significance, indicating that they were dysregulated (Figure 11D). We further analyzed the survival differences between high and low methylation of PRGs, and the results showed that high methylation of APOC1, ZCCHC12, NCAPG, IFI30, and FXYD1 was associated with a lower risk of survival in LGG, while high methylation of BMP2 indicated a higher risk of survival. Additionally, high methylation of IFI30 and ZCCHC12 exhibited a correlation with a lower survival risk in GBM (Figure 11E). The influence of methylation levels at various cytosine-phosphate-guanine (CpG) sites on the prognosis in LGG and GBM was further investigated utilizing MethSurv. The results indicated that patients with LGG exhibiting higher methylation levels of BMP2 demonstrated inferior OS compared to those with lower methylation levels (Figure 11F). Conversely, patients with LGG who exhibited higher methylation levels of APOC1, ZCCHC12, NCAPG, IFI30, and FXYD1 displayed superior OS compared to those with lower methylation levels (Supplementary Figure 12A-E). Additionally, patients with GBM who harbored higher methylation levels of IFI30 and ZCCHC12 exhibited better OS compared to those with lower methylation levels (Supplementary Figures 13A, B). In conclusion, there was a negative correlation between the methylation levels of APOC1, NCAPG, and IFI30 and their expression levels. It can be inferred that the reduced methylation of APOC1, NCAPG, and IFI30 contributed to their elevated expression in glioma, which eventually led to a poor prognosis for glioma patients.

To assess the potential of PRRS as a biomarker for predicting the clinical efficacy of immunotherapy, we compared the distribution of TIDE score, MSI score, and TMB score in different risk groups. According to our findings, glioma patients who showed a high efficacy of immunotherapy responses had significantly lower PRRS than those who were immunotherapy resistant (Figure 12A). The TIDE and TMB scores were higher in the high-risk group, whereas the MSI score was the opposite (Figures 12B–D). Furthermore, the Kaplan-Meier survival analysis revealed that patients with elevated TIDE scores exhibited decreased OS rates compared to those with lower TIDE scores (Figure 12E). However, patients with high MSI scores demonstrated enhanced OS rates in comparison to those with low MSI scores (Figure 12F). Conversely, patients with high TMB scores displayed diminished OS rates when compared to those with low TMB scores (Figure 12G). In addition, by integrating the TIDE score with PRRS, we stratified glioma patients into four distinct subgroups. Similarly, using the same methodology, we combined the MSI and TMB scores with PRRS, respectively, uncovering marked disparities in survival outcomes (Figures 12H–J). TIDE and PRRS appear to exert a joint and common influence on the survival outcomes of glioma patients, with TMB similarly sharing an effect with PRRS in this regard. These findings indicate that glioma patients with high PRRS may be more prone to cancer immune evasion. To delve into the role of PRRGs in the immune responses and mechanisms within TME, we evaluated the correlation between PRRGs and TMB as well as MSI. The results demonstrated a positive correlation between APOC1, F2R, PTBP1, NCAPG, and IFI30 (P< 0.001) with TMB, whereas BMP2 and PDE2A (P< 0.001) exhibited a negative correlation with TMB (Supplementary Figure 14A). On the contrary, APOC1, F2R, PTBP1, NCAPG, and IFI30 (P< 0.001) were negatively correlated with MSI, while BMP2 and PDE2A (P< 0.001) were positively correlated with MSI (Supplementary Figure 14B).

To fully look into the potential worth of PRGs as new therapeutic targets, we first selected a number of drugs from the CTRP and GDSC datasets, which presented a significant association between drug sensitivity and the PRRS model (Figures 13A, B). In the high-risk group, the IC50 values of AT-7519, BIX02189, CUDC-101, PIK-93, THZ-2-102-1, Trametinib, 17-AAG, and PD-0325901 were significantly lower than those in the low-risk group, indicating that the sensitivity of glioma patients in the high-risk group to these drugs was higher than that in the low-risk group (Figure 13C). Subsequently, the lollipop plots demonstrated a statistically significant correlation between the IC50 values of these drugs and nine PRGs (Supplementary Figure 15). These findings suggested that these drugs could potentially serve as therapeutic options for glioma by modulating the expression of PRGs products. Utilizing both AutoDock and AutoDock Vina, we performed molecular docking between the proteins corresponding to the 9 PRGs and the drugs mentioned, with the results presented in Supplementary Table 2. We selected the best small molecules that bind to each protein and then used the PyMol process for visual analysis. The results indicated the following optimal docking binding energy: -5.4 kcal/mol for the binding energy between APOC1 and THZ-2-101-1, -6.1 kcal/mol for FXYD1 and AT-7519, -8.0 kcal/mol for ZCCHC12 and THZ-2-101-1, -11.0 kcal/mol for F2R and THZ-2-101-1, -8.4 kcal/mol for PTBP1 and BIX02189, -8.2 kcal/mol for NCAPG and THZ-2-101-1, -7.5 kcal/mol for BMP2 and THZ-2-101-1, -10.8 kcal/mol for PDE2A and THZ-2-101-1, and -7.3 kcal/mol for IFI30 and THZ-2-101-1 (Figures 14A–I).

The outcomes of MD simulations, encompassing RMSD and RMSF, constitute crucial evidence for assessing the stability of protein-ligand complexes. A lower RMSD value indicates a higher degree of stability (39). Moreover, the computation of RMSF offers deep insights into the fluctuations of protein residues throughout the simulation, influencing protein functionality, with elevated values signifying enhanced residue flexibility and diminished values indicating limited motion (40). We performed MD simulations on these 9 molecular docking systems and visualized their output data (Supplementary Figure 16). The RMSD results obtained from MD simulations revealed suboptimal overall performance for the NCAPG and THZ-2-101–1 complex system. Significant RMSD fluctuations suggested that the system failed to demonstrate effective inhibition of glioma (Supplementary Figure 16F). The RMSD of the FXYD1-AT-7519 complex exhibits initial fluctuations but stabilizes after 15 ns, whereas the RMSD of the PTBP1-BIX02189 complex stabilizes after 10 ns (Supplementary Figures 16B, E). It was also observed that APOC1, ZCCHC12, F2R, BMP2, PDE2A, and IFI30 could form stable complex molecular dynamics systems with THZ-2-101-1, respectively (Supplementary Figures 16A, C, D, G-I). Supplementary Figures 16J-R displays the fluctuations of residues pertaining to nine palmitoylation-associated proteins, respectively. These results suggested that APOC1, FXYD1, ZCCHC12, F2R, PTBP1, BMP2, PDE2A, and IFI30 might be new therapeutic targets for glioma.

To validate the anti-glioma effects of the aforementioned compounds, we conducted a series of in vitro experiments employing LN229 and U251 glioma cells. Firstly, we employed the CCK-8 assay to assess the cytotoxicity. In LN229 cells, a concentration-dependent decline in cell viability was observed after 48 hours of treatment with AT-7519, BIX02189, and THZ-2-101-1, with IC50 values of 406.3, 89.5, and 240.6 μm/mL, respectively. Similarly, in U251 cells, the IC50 values for AT-7519, BIX02189, and THZ-2-101–1 were determined to be 380.1, 127.7, and 50.7 μm/mL, respectively (Figure 15A). Wound healing assays were employed to assess the effects of the three compounds on cell migration. Following 24 and 48 hours of treatment, a notable decrease in wound closure was evident in the drug-treated groups relative to the control, thereby underscoring the capacity of these compounds to suppress glioma cell motility (Figure 15B). Furthermore, following 24 hours of drug treatment, we employed flow cytometry to assess cell apoptosis. When compared to the cells in the control group, those treated with AT-7519, BIX02189, and THZ-2-101–1 demonstrated a notably elevated apoptotic index. (Figure 15C). Collectively, these findings confirm that AT-7519, BIX02189, and THZ-2-101–1 can inhibit the migration of glioma cells in vitro and also promote their apoptosis, thereby supporting their potential as effective anti-glioma agents.