We obtained transcriptome data and clinical information for low-grade glioma (LGG) from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) using the “TCGAbiolinks^” (16) R package, including data from a total of 504 LGG patients with survival information. Additional LGG validation data were sourced from the CGGA-LGG693 and CGGA-LGG325 datasets in the China Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn), where we selected patients classified as WHO grade II, resulting in 172 patients in the CGGA-LGG693 cohort and 98 patients in the CGGA-LGG325 cohort. Furthermore, single-cell transcriptomic data (GSE117891) (17) were collected from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), encompassing samples from 73 regions in 13 patients with glioma and 1 patient with brain metastases. Gene expression data from 1,152 cases of normal brain tissue were obtained from the Genotype-Tissue Expression (GTEx) project (https://commonfund.nih.gov/GTEx). For the LGG pediatric peptide vaccine immunotherapy cohort (12 patients) versus the GBM (Glioblastoma) anti-PD-1 immunotherapy cohort (34 patients), data were sourced from the TIGER database (http://tiger.canceromics.org/).

We utilized the R package “Seurat” (18) for the analysis of single-cell RNA sequencing (scRNA-seq) data from Glioma patients. Initially, we conducted quality control measures, which included filtering out low quality genes detected in fewer than three cells, removing low-quality cells with fewer than 300 identified genes, excluding cells with mitochondrial gene content greater than 15%, ribosomal gene content less than 3%, and hemoglobin gene content greater than 0.1%. We also filtered out the MALAT1 housekeeping gene and mitochondrial genes, and removed doublets. For the retained cells, we normalized the gene expression matrix using the NormalizeData function in Seurat, and applied centering and scaling using the ScaleData function. We extracted the top 2000 highly variable genes and then performed principal component analysis (PCA), utilizing the top 20 principal components for clustering. The cell populations were visualized using Uniform Manifold Approximation and Projection (UMAP). Clustering analysis was performed using the FindClusters function with a resolution of 0.3 to identify distinct cell clusters.

To identify cell types, we employed specific cell markers sourced from the official CellMarker 2.0 website (http://117.50.127.228/CellMarker/). The FindAllMarkers function in Seurat was used to determine differentially expressed genes across various cell clusters.Using the default method of the FindAllMarkers function, the Wilcoxon rank-sum test, we screened for differentially expressed genes in the microglia population, selecting those with p.adjust < 0.05 and an absolute log2FC > 0.585.

To calculate the glycosylated gene sets enrichment score within each cell cluster, we employed the irGSEA.score function from the “irGSEA” (19) package. The KEGG glycosylated gene sets (C2 classification), including “KEGG_N_GLYCAN_BIOSYNTHESIS” and “KEGG_O_GLYCAN_ BIOSYNTHESIS”, was downloaded from the MSigDB database, and several enrichment analysis algorithms were applied to normalized RNA sequencing data, including AUCell, UCell, and ssGSEA.

We computed the KEGG glycosylation pathway enrichment score for each patient in the TCGA-LGG cohort using GSVA (20). Patients were divided into high and low N-glycosylation groups based on the median N-glycosylation score, and similarly into high and low O-glycosylation groups based on the median O-glycosylation score. Differential expression analysis was performed using “Deseq2” (21), with thresholds of P.adjust < 0.05 and the absolute value of log2FC > 0.585 to identify differentially expressed genes. These genes were used to investigate the molecular characteristics associated with high and low activity in the N-glycosylation and O-glycosylation pathways.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were conducted on the differentially expressed genes using the R package “clusterProfiler” (22). Only results with p.adjust < 0.05 were considered significantly enriched.

The ESTIMATE (23) (Estimation of Stromal and Immune Cells in Malignant Tumors Using Expression Data) algorithm was utilized to assess changes in immune scores, stromal scores, ESTIMATE scores, and tumor purity in the samples. Additionally, the relative abundance of various cell types was determined using the “CIBERSORT” (24) and “GSVA” (20) R packages, employing the CIBERSORT and ssGSEA (Single Sample Gene Set Enrichment Analysis) methods, respectively.

To establish a glycosylation-based prognostic signature, we performed univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression with 10-fold cross-validation using the R package “glmnet” (25). In the LASSO regression, we selected “lambda.min” to prevent overfitting. The risk score for each LGG patient was calculated using the following formula:

Risk score = (0.382 × ASPM expression) + (0.017 × CHI3L1 expression) + (0.068 × LILRA4 expression) + (0.258 × MSN expression) + (0.116 × OCIAD2 expression) + (-0.137 × PTGER4 expression) + (0.109 × SERPING1 expression) + (0.073 × TNFRSF12A expression)

Kaplan-Meier analysis compared overall survival between high and low risk score groups, stratified by the median risk score. Receiver Operating Characteristic (ROC) curves generated with the “timeROC” (26) R package assessed the predictive accuracy of risk score for 1-, 3-, and 5-year survival.

Prognostic nomograms incorporating risk score and other clinical features were constructed using the R package “rms”. The performance of these nomograms was evaluated through calibration curves and ROC curves to assess their predictive accuracy and reliability.

We combined the expression data from the TCGA-LGG cohort with expression data from normal brain tissues in the GTEx database and performed a Wilcoxon test using the R package “rstatix”. Boxplots were then generated to visualize gene expression levels. Additionally, we calculated the proportion of each cell type in each sample from the single-cell data and determined the average gene expression within each cell type group. The gene with the highest average expression was used to compute the correlation between average gene expression and cell content in each sample, using the cor.test function from base R. The cor.test function employed Pearson’s correlation coefficient to assess the linear relationship between gene expression and cell content.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was utilized to predict the potential response of LGG patients to immune checkpoint inhibitors (ICIs) treatment (27). A higher TIDE score suggests a decreased likelihood of benefiting from immunotherapy and an increased risk of immune escape. The Dysfunction score indicates the level of immune cell dysfunction infiltrating the tumor, while the Exclusion score reflects the degree of immune cell exclusion within the tumor microenvironment.

To predict susceptibility to different drugs in LGG patients within the high and low risk score groups, drug treatment sensitivity was assessed using the R package “oncoPredict” (28). This analysis used expression matrices from the Cancer Drug Sensitivity Genomics (GDSC) database and drug treatment information as a training set. The half-maximal inhibitory concentration (IC50) served as the primary indicator for evaluating the sensitivity of LGG cancer cells drugs. Differences in IC50 between the high and low risk score groups were compared using the Wilcoxon test, with p < 0.05 considered statistically significant.

All statistical analyses were conducted using R software (version 4.3.0). The R packages employed in this study are open access and freely available. Statistical significance was determined at a threshold of p < 0.05, with levels indicated as follows: * p < 0.05; ** p < 0.01; *** p < 0.001.

Figure 1 shows the flow chart of this study. In the GSE117891 dataset, we analyzed a total of 5,368 cells and 22,559 genes. By collecting common cell marker genes for brain tissue from the CellMarker 2.0 website and combining them with classical cell markers, we ultimately annotated four cell types: oligodendrocyte (MBP, MOP, PLP1, MAG), T cell (CD3D, CD3E, CD8A), astrocyte (GFAP, AQP4, SOX9, CLU), and microglial cell (TMEM119, CX3CR1, P2RY12). Oligodendrocyte accounted for 744 cells, approximately 13.7%, T cell accounted for 379 cells, approximately 7.1%, astrocyte accounted for 2,993 cells, approximately 55.8%, and microglial cell accounted for 1,252 cells, approximately 23.4%. These cell types were visualized using UMAP and DotPlot ( Figure 2A detailed in Figure 2B and Supplementary Figure 1A ).

Based on the expression levels of genes from the “KEGG_N_GLYCAN_BIOSYNTHESIS” and “KEGG_O_GLYCAN_BIOSYNTHESIS” pathways in the KEGG pathway database, we calculated glycosylation scores for each cell type using three algorithms: AUCell, UCell, and ssGSEA. The results showed that, for both N-glycosylation and O-glycosylation, microglia exhibited the highest glycosylation scores across all three algorithms, significantly higher than other cell types ( Figure 2C ; Supplementary Figure 1B ). This suggests that glycosylation may play an important role in the function of microglia.

To understand the unique molecular characteristics of microglia, we identified 1,684 differentially expressed genes (DEGs) compared to other cell types and visualized these genes using a volcano plot ( Figure 2D ). GO and KEGG enrichment results ( Figures 2E, F ) showed that these DEGs were enriched in pathways related to cytokine regulation, immune response activation, ATP synthesis coupled with electron transport, and other related signaling pathways.

Kaplan-Meier survival analysis revealed significant survival differences based on N-glycosylation pathway activities (P < 0.001, Figure 3A ), indicating that patients with lower N-glycosylation pathway activity had better prognoses. We identified 3,019 DEGs between high and low N-glycosylation pathway activity subgroups ( Figure 3B ) and performed functional enrichment analyses using GO and KEGG, which highlighted potential alterations in immune signaling pathways ( Figures 3C, D ). Furthermore, we evaluated the immune microenvironment of each patient. Results indicated that patients with high N-glycosylation pathway activity exhibited elevated stromal scores, immune scores, and ESTIMATE scores, along with reduced tumor purity, indicative of increased tumor heterogeneity and a more complex immune microenvironment ( Figure 3E ).

There was no significant difference in prognosis between high and low O-glycosylation groups ( Figure 4A ). Similarly, enrichment analysis of 3,164 DEGs between high and low O-glycosylation pathway activity groups ( Figure 4B ) revealed that these genes were enriched in signaling pathways related to chromosome segregation, synaptic signaling regulation, and the cell cycle ( Figures 4C, D ). The high O-glycosylation pathway activity group showed higher matrix scores, immune scores, and ESTIMATE scores compared to the low-activity group, suggesting a similarly complex immune microenvironment ( Figure 4E ).These findings underscore the potential role of N-glycosylation and O-glycosylation levels in influencing immune cell types and responses within the LGG tumor microenvironment.

Recognizing the impact of glycosylation on the prognosis and immune microenvironment of patients with LGG, we further explored prognostic features to understand the potential complexity of LGG. We compared the DEGs of microglia in GSE117891 with other cell type, the DEGs of the high-low N-glycosylation subgroup in the TCGA-LGG cohort, and the DEGs of the high-low O-glycosylation subgroup, and identified 125 intersecting genes ( Figure 5A ). Through univariate Cox regression analysis, 105 of 125 genes were confirmed to be associated with prognosis. After LASSO regression narrowing, eight key genes were identified ( Supplementary Figures 2A, B ), constructing a prognostic feature model called risk score ( Figure 5B ), including ASPM, CHI3L1, LILRA4, MSN, OCIAD2, PTGER4, SERPING1, and TNFRSF12A.

We calculated the risk score for each patient in the TCGA-LGG cohort and divided the patients into high and low risk score groups based on the median risk score, showing the survival status plots of the patients. The results indicated that patients with higher risk score values exhibited increased expression levels of ASPM, CHI3L1, LILRA4, MSN, OCIAD2, PTGER4, SERPING1, and TNFRSF12A ( Figure 5C ), which were associated with poorer survival outcomes ( Figure 5D ). Kaplan-Meier analysis showed that overall survival (OS) was significantly better in the low risk score group (P < 0.001, Figure 5E ).

To assess the predictive power of our prognostic features, we generated receiver operating characteristic (ROC) curves ( Figure 5H ) for 1-, 3-, and 5-year OS, with areas under the curve (AUC) of 0.882, 0.859, and 0.785, respectively, indicating good predictive performance of our model. Additionally, we validated our prognostic features in two independent external validation sets (CGGA-LGG693, CGGA-LGG325), yielding satisfactory results with 5-year AUCs of 0.669 and 0.634, respectively ( Figures 5F, G, I, J ).

To improve clinical utility, we created nomograms combining risk score and clinical characteristics. Univariate and multivariate Cox regression analyses, which assessed the effect of clinical characteristics and risk score on LGG survival, showed that age and risk score were important predictors of survival ( Figures 6A, B ). Based on these predictors, we designed a nomogram to estimate the 1-, 3-, and 5-year survival probabilities of patients with LGG ( Figure 6C ). The ROC curve showed that the nomogram had an AUC of 0.906 at 1 year, 0.898 at 2 years, and 0.819 at 3 years ( Figure 6D ), and the calibration curve demonstrated the predictive accuracy of the nomogram ( Supplementary Figure 2H ). In addition, we validated the predictive performance of nomograms in two independent external validation sets ( Figures 6E, F ; Supplementary Figures 2I, J ).

Collectively, these results suggest that risk score-based nomograms, rooted in the concept of glycosylation, have substantial clinical predictive value for LGG.

The analysis of the TCGA-LGG cohort combined with GTEx normal brain tissue data revealed notable expression differences of the prognostic feature genes. Specifically, ASPM, LILRA4, MSN, PTGER4, and TNFRSF12A were significantly upregulated in LGG patients compared to normal tissues, while CHI3L1, OCIAD2, and SERPING1 exhibited decreased expression ( Figure 7A ). We further examined the distribution of these eight genes in the GSE117891 dataset ( Figure 7B, C ). The results indicated that ASPM, CHI3L1, and TNFRSF12A were more highly expressed in astrocyte, while LILRA4, MSN, PTGER4, and SERPING1 were more highly expressed in microglial cell. OCIAD2 was found to be more highly expressed in T cell. We calculated the proportion of each cell type in each sample and determined the average expression of genes in the cell group with the highest expression. We then analyzed the correlation between gene expression and cell proportion ( Figure 7D, E ). The results indicated that ASPM expression was positively correlated with the proportion of astrocyte (cor = 0.61, p < 0.05), OCIAD2 expression was negatively correlated with the proportion of oligodendrocyte (cor = -0.59, p < 0.05), and SERPING1 expression was positively correlated with the proportion of T cell (cor = 0.57, p < 0.05).

We also investigated the differences in the immune microenvironment between patients in different risk score groups within the TCGA-LGG cohort. The findings showed a positive correlation between glycosylation score and risk score. Patients in the high risk score group had higher stromal, immune, and ESTIMATE scores compared to those in the low risk score group, but exhibited lower tumor purity ( Figure 7F ). This indicates a more complex immune microenvironment, characterized by increased tumor heterogeneity.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm evaluates the potential for tumor immune escape by analyzing gene expression profiles of tumor samples. Our analysis revealed that in the TCGA cohort, the high risk score group exhibited higher TIDE, Exclusion, and Dysfunction scores, along with lower MSI (microsatellite instability) scores. This suggests that patients in the high risk score group are more likely to experience immune escape and may have a poorer response to immune checkpoint inhibitor (ICI) therapy ( Figure 8A ). Conversely, the low risk score group predicted a higher proportion of patients likely to respond positively to immunotherapy ( Figure 8B ). These findings indicate that patients with a low risk score may derive greater benefit from immunotherapy compared to those with a high risk score.

In the LGG pediatric peptide vaccine immunotherapy cohort, the high risk score group predicted a higher proportion of patients responding to immunotherapy ( Figure 8C ). In contrast, in the GBM anti-PD-1 immunotherapy cohort, the low risk score group indicated a higher proportion of patients responding to immunotherapy ( Figure 8D ).The ROC of risk score predicting immunotherapy response was 0.77 and 0.626 in the two cohorts, respectively ( Figures 8E, F ).

Furthermore, we examined the relationship between risk score and drug therapy effectiveness in LGG treatment. Our findings revealed that high risk score values were associated with lower IC50 of Temozolomide and Gemcitabine, and higher IC50 of Gefitinib ( Figures 8G–I , P < 0.05). Consequently, our study suggests that risk score could serve as a valuable predictor of drug therapy sensitivity in LGG patients.

In the TCGA-LGG cohort, we analyzed the survival outcomes of patients receiving medication versus those not receiving medication. At the 2-year survival time point, we observed distinct trends for short-term (< 2 years) and long-term (≥ 2 years) prognosis between the two groups ( Figures 9A–C ). Specifically, among all LGG patients, those receiving medication had a better short-term prognosis, while those not receiving medication had a better long-term prognosis.

Further analysis revealed that this phenomenon was particularly pronounced in the high risk score group. In this group, patients receiving medication had a better short-term prognosis but a worse long-term prognosis ( Figures 9D–F ). Conversely, in the low risk score group, the prognosis was not significantly affected by whether the patients received medication ( Figures 9H–J ).

The independent validation cohort CGGA-325 yielded results consistent with those of the TCGA-LGG cohort, reinforcing these observations ( Supplementary Figure 3 ). These findings suggest that the risk score not only serves as a prognostic indicator but may also help identify which patients are likely to benefit from pharmacological treatment in the short term and require careful consideration for long-term prognosis. However, consistent results were not observed in the CGGA-693 cohort, and further studies are needed to validate these findings( Supplementary Figure 4 ).