Development and validation of a two glycolysis-related LncRNAs prognostic signature for glioma and in vitro analyses

Abstract

Background

Mounting evidence suggests that there is a complex regulatory relationship between long non-coding RNAs (lncRNAs) and the glycolytic process during glioma development. This study aimed to investigate the prognostic role of glycolysis-related lncRNAs in glioma and their impact on the tumor microenvironment.

Methods

This study utilized glioma transcriptome data from public databases to construct, evaluate, and validate a prognostic signature based on differentially expressed (DE)-glycolysis-associated lncRNAs through consensus clustering, DE-lncRNA analysis, Cox regression analysis, and receiver operating characteristic (ROC) curves. The clusterProfiler package was applied to reveal the potential functions of the risk score-related differentially expressed genes (DEGs). ESTIMATE and Gene Set Enrichment Analysis (GSEA) were utilized to evaluate the relationship between prognostic signature and the immune landscape of gliomas. Furthermore, the sensitivity of patients to immune checkpoint inhibitor (ICI) treatment based on the prognostic feature was predicted with the assistance of the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. Finally, qRT-PCR was used to verify the difference in the expression of the lncRNAs in glioma cells and normal cell.

Results

By consensus clustering based on glycolytic gene expression profiles, glioma patients were divided into two clusters with significantly different overall survival (OS), from which 2 DE-lncRNAs, AL390755.1 and FLJ16779, were obtained. Subsequently, Cox regression analysis demonstrated that all of these lncRNAs were associated with OS in glioma patients and constructed a prognostic signature with a robust prognostic predictive efficacy. Functional enrichment analysis revealed that DEGs associated with risk scores were involved in immune responses, neurons, neurotransmitters, synapses and other terms. Immune landscape analysis suggested an extreme enrichment of immune cells in the high-risk group. Moreover, patients in the low-risk group were likely to benefit more from ICI treatment. qRT-PCR results showed that the expression of AL390755.1 and FLJ16779 was significantly different in glioma and normal cells.

Conclusion

We constructed a novel prognostic signature for glioma patients based on glycolysis-related lncRNAs. Besides, this project had provided a theoretical basis for the exploration of new ICI therapeutic targets for glioma patients.

Introduction

Gliomas make up virtually 80% of all lethal primary brain tumors, which seriously threatening human health and causing a heavy burden to the social economy [18]. Despite many advances in deciphering the underlying molecular mechanisms of gliomas, the efficacy of clinical comprehensive treatment options have reached a bottleneck, and as a result, the long‐term survival rate of glioma patients remains poor. The complex heterogeneous of tumor and the unique microenvironment of the brain present major challenges in treating gliomas. Therefore, there remains an urgent need to comprehensively understand the tumor microenvironment and identify a valuable biomarker for predicting the prognosis of glioma patients, which may lead to the development of new potential therapeutics for glioma patients.

Abnormal metabolism and immune evasion are two hallmarks of cancer [20]. It has been under extensive exploration in the hope of discovering new targets and effective therapies. Like other cancers, glioma presents a unique metabolic state known as the Warburg effect [21]. The exact mechanism of metabolic transformation remains unclear, previous studies have showed that the metabolic mode can switch during tumor progression. Reversal of the Warburg effect could potentially serve as a novel therapy for glioma [21, 24]. Recently, studies have identified alterations in tumor metabolism can also contribute to a potent tolerogenic immune environment [15]. With continued advancement in both of these research disciplines, the relationship between tumor metabolism and their subsequent influence on immune regulation has become increasingly recognized as an important factor contributing to tumor growth and progression. A recently study showed that there is a close correlation between the altered metabolic landscape and increased activity of infiltrated immune cells within the tumor microenvironment [17, 27]. Furthermore, CTLA-4 blockade has been found to promotes metabolic fitness and the infiltration of immune cells, especially in glycolysis-low tumours [31]. These findings suggest that the complex interdependencies exist between tumor metabolic and immune responses.

LncRNAs have been reported to be dysregulated in various types of cancers, it can regulate cancer cell proliferation, invasion, metastasis, and therapeutic resistance [11]. Tumor glycolysis could be used as a potential therapeutic target for cancer [8], the major challenge lies in the fact that metabolism is a universal cellular process [23]. Accumulating evidence has shown that lncRNAs can alter glucose metabolism either directly or indirectly [16, 29, 30, 33]. Therefore, lncRNAs are considered as a promising strategy for addressing this challenge [2, 7, 13].The relationship between glycolysis-related lncRNAs and glioma prognosis has rarely been studied. A recent study identified six glycolysis-related lncRNAs significantly related to prognosis of glioma patients through analyzed CGCA database [25, 26]. However, comprehensive analysis focusing on glycolysis-related lncRNAs exerts on glioma patients in TCGA, as well as whether these glycolysis-related lncRNAs have potential impact on immune or not is lacking. In the present study, we implemented studies with transcriptome and clinical data of glioma from The Cancer Genome Atlas (TCGA) projects to explore the glycolysis-related lncRNAs and immunity in glioma, to identify the glycolysis-related lncRNAs and relationship between the tumor glycolysis and immune cells infiltration on glioma patients.

Results

Consistent clustering of glioma and the survival rate of clusters

The study was conducted as described in the flow chart (Fig. 1). Consistent clustering analysis of glioma patients was performed based on the expression of glycolysis genes in the TCGA database. Finally, patients were decided to be divided into 2 sample clusters, with cluster 1 containing 132 samples and cluster 2 containing 512 samples (Fig. 2A). The K-M suggested a significant difference in cumulative mortality at 20 years between the two subtypes (P < 0.0001; Fig. 2B). In the glioma cohort, the best outcome was found in patients classified as cluster 2. The next question was to determine the relationship between sample clustering results and clinicopathological data (treatment, age, gender, IDH, grade, and MGMT). Matching these results to all clinicopathological characteristics, significant correlations had been noted among treatment, age, IDH, grade, MGMT and the two clusters (Table 1). In summary, these results appeared to imply that glioma could be stratified according to the expression pattern of glycolysis genes.

Fig. 1

Flow chart of the study design

Fig. 2

A Identification of consistent clustering by glycolysis genes in the TCGA database. B The “Kaplan–Meier”overall survival (OS) curve of two clusters defined by consistent expression

Table 1 The relationship between sample clusters and clinical characteristics of glioma patients

Furthermore, we also assessed differences in pathway enrichment scores between the two clusters based on the KEGG pathway using the ssGSEA algorithm (Additional file 6: Table S6). Unexpectedly, glycolytic (‘Glycolysis/Gluconeogenesis’) and glycan synthesis and metabolic (‘Amino sugar and nucleotide sugar metabolism’, ‘Glycosaminoglycan biosynthesis-keratan sulfate’, ‘Mannose type O-glycan biosynthesis’, ‘Galactose metabolism’, etc.) pathways were significantly different between the two clusters. Also, immune response- (‘cAMP signaling pathway’, ‘IL-17 signaling pathway’, ‘Cell adhesion molecules’, ‘Intestinal immune network for IgA production’, etc.) and immune cell (‘Th17 cell differentiation’, ‘Th1 and Th2 cell differentiation’, etc.)-related pathways were notable between Cluster 1 and Cluster 2. Further, PD-L1 expression and PD-1 checkpoint pathway in cancer also showed substantial variations.

Establishment of a prognostic signature based on the DE-glycolysis-related lncRNAs

We identified a total of 2 DE-lncRNAs from both subtypes, relative to cluster 2, lncRNA FLJ16779 was the up-regulated gene and lncRNA AL390755.1 was the down-regulated gene (Fig. 3A; Additional file 7: Table S7). These genes were defined as DE-glycolysis-related lncRNAs for the subsequent analysis. To verify whether these 2 DE-lncRNAs were the risk factors for glioma, we performed a univariate Cox regression analysis in the training set. Coincidentally, all DE-lncRNAs met P < 0.05 (Fig. 3B). Subsequently, a multifactorial Cox regression analysis with a STEP function indicated these 2 lncRNAs as the optimal variables for constructing a prognostic signature (Fig. 3C). Specifically, AL390755.1 (P = 8.59E−14) was possible a risk factor for glioma (HR = 1.38), whereas FLJ16779 (P = 1.82E−6) was inferred to be the oncogene (HR = 0.77).

Fig. 3

A Volcano plot and heatmap of differentially glycolysis-related lncRNAs in both cluster 1 and cluster 2 subtypes. B Univariate cox regression analysis in the training set. C Multifactorial cox regression analysis in the training set

Assessment and validation of the effectiveness of the prognostic feature

Risk scores were calculated separately for each sample in the training and testing sets according to the aforementioned formula, and patients were categorized into high- and low-risk groups based on the median risk score of each set. In both the training and testing sets, with increasing risk scores, the number of patient deaths climbed sharply (Fig. 4A and D). Subsequently, K-M-survival analysis was performed on this prognostic feature, while AUC values were calculated by time-dependent ROC analysis for the accuracy of which in predicting survival, using the outcome variable. In the training set, K–M curves could effectively distinguish between high-risk and low-risk groups (Fig. 4B), and the risk scores all reached an AUC of 0.8 or more at 1 to 5 years (Fig. 4C). Similar results were also reproduced in the testing set (Fig. 4E and F). Based on the expression heatmap of the prognostic lncRNAs, which was observed in the high-risk group associated with poorer OS, AL390755.1 was overexpressed. In contrast, FLJ16779 was overexpressed in the low-risk group possessing a longer OS characteristic. Moreover, we found that the vast majority of Cluster 2 patients were in the low-risk group. These results suggested that the prognostic signature was highly sensitive and specific and that these signature lncRNAs could be utilized as prognostic biomarkers in clinically. Furthermore, the statistical tables of clinical information in the training and testing sets were displayed in Table 2.

Fig. 4

Construction of the glycolysis-related LncRNAs signature for survival prediction. A, D Risk score distribution, survival status for patients in low- and high-risk groups by the LncRNA signature in training and testing sets. B, E Kaplan–Meier curve based on the identified survival-related glycolysis LncRNAs in training and testing sets. C, F ROC curves of the signature for predicting1, 3, 5, 7- year survival of glioma in training and testing sets

Table 2 Clinical information of glioma patients in the training and testing sets

Independent prognostic analysis

Here, risk scores and clinicopathological characteristics (age, gender, treatment type, cluster, MGMT, grade, and IDH) of the full Gliomas sample in the TCGA database were included in the Cox analysis to explore their potential for independent prognosis. Ultimately, age, IDH, MGMT, grade, and risk score were identified as independent prognostic factors for glioma through univariate and multivariate Cox analyses (Fig. 5A and B). These independent prognostic factors were then incorporated into the Nomogram to explore the prediction of the Nomogram model for patient survival at 1, 3, and 5 years, which with a c-index of 0.8705766 (Fig. 5C). The status of each variable corresponded to a score, with a higher total score for a patient indicating poorer survival for that patient. Besides, the calibration curves also demonstrated that the Nomogram model was effective in predicting the survival of patients (Fig. 5D).

Fig. 5

A, B Univariate (A) and multivariate (B) Cox analyses of clinicopathological characteristics and risk score. C Nomogram to predict the 1-, 3-, 5-year OS. D The calibration curve of the nomogram model prediction

Functional enrichment analysis of risk score-related DEGs

To further reveal the potential function of prognostic lncRNAs, we first screened 680 risk score-related DEGs between high- and low-risk groups (Additional file 8: Table S8) and performed GSEA on them to reveal potential functions. Figure 6A and B illustrated the top10 GO and KEGG results, respectively. GO analysis revealed that these DEGs were mainly enriched in terms related to immune responses (‘ADAPTIVE IMMUNE RESPONSE’, ‘ADAPTIVE IMMUNE RESPONSE BASED ON SOMATIC RECOMBINATION OF IMMUNE RECEPTORS BUILT FROM IMMUNOGLOBULIN SUPERFAMILY DOMAINS’, ‘HUMORAL IMMUNE RESPONSE’, ‘NEGATIVE REGULATION OF IMMUNE SYSTEM PROCESS’, etc.) and immune cell regulation (‘LEUKOCYTE PROLIFERATION’, ‘REGULATION OF LYMPHOCYTE ACTIVATION’, ‘REGULATION OF T CELL ACTIVATION’, etc.). Notably, a multitude of biological processes associated with neurons (‘NEURON PROJECTION TERMINUS’, ‘NEURON TO NEURON SYNAPSE’, ‘REGULATION OF NEURON PROJECTION DEVELOPMENT’, etc.), neurotransmitters (‘NEUROTRANSMITTER RECEPTOR ACTIVITY’, ‘NEUROTRANSMITTER SECRETION’, ‘NEUROTRANSMITTER TRANSPORT’, etc.), and synapses (‘POSTSYNAPTIC MEMBRANE’, ‘POSTSYNAPTIC NEUROTRANSMITTER RECEPTOR ACTIVITY’, ‘POSTSYNAPTIC SPECIALIZATION MEMBRANE’, ‘PRESYNAPTIC MEMBRANE’, etc.) were significantly enriched. KEGG results indicated that risk score-related DEGs were significantly associated with immune disorders such as ‘AUTOIMMUNE THYROID DISEASE’ and ‘SYSTEMIC LUPUS ERYTHEMATOSUS’. Consistently, immune response-related pathways (‘CYTOKINE CYTOKINE RECEPTOR INTERACTION’, ‘PRIMARY IMMUNODEFICIENCY’, ‘NATURAL KILLER CELL MEDIATED CYTOTOXICITY’, etc.) were also significantly enriched. Moreover, the ‘AMINO SUGAR AND NUCLEOTIDE SUGAR METABOLISM’ pathway was inevitably featured in the results. More detailed results of GO and KEGG analysis were reported in Additional file 9: Table S9 and Additional file 10: Table S10.

Fig. 6

A, B Top10 GO (A) and KEGG (B) enrichment results by GSEA enrichment analysis. C The lncRNA-mRNA co-expression network by pearson correlation analysis. Two lncRNAs (pink rhombus) are connected to their respective co-expressed (|cor|≥ 0.6 and P < 0.05) mRNAs (blue circles) with lines (pink: postive correlation; purple: negative correlation)

Moreover, we illustrated the network of interactions of 150 risk score-related DEGs with prognostic lncRNAs for |cor|≥ 0.6 and P < 0.05 in Fig. 6C. Among them, 35 genes were common DEGs for both prognostic lncRNAs (Table 3). Elaborately, ARPP21, JPH3, NTNG2, KCNIP3, and ELFN2 were positively correlated with lncRNAs FLJ16779, but negatively associated with AL390755.1. The remaining 30 genes such as EMP3, TNFRSF12A, MCUB, and ANXA2 were positively associated with lncRNA AL390755.1, but negatively related to FLJ16779.

Table 3 List of common risk score-related DEGs for prognostic lncRNAs (|cor|≥ 0.6 and P < 0.05)

The impact of risk scores on the immune landscape of glioma patients

Inspired by the results of the ssGSEA, we hypothesized that prognostic lncRNAs may operate in the patient’s immune microenvironment to influence patient outcomes. ESTIMATE analysis indicated that the immune, stromal, and ESTIMATE scores were significantly higher in the high-risk group than in the low-risk group, suggesting that the high-risk group had more components of the immune microenvironment (Fig. 7A). Subsequent immune cell enrichment analysis demonstrated that more immune cells, such as Natural killer T cells, Myeloid-derived suppressor cell (MDSC), and Type 1T helper cell, were presented in the high-risk group (Fig. 7B). This evidence suggested that prognostic lncRNAs were involved in the altered immune microenvironment of patients.

Fig. 7

A Comparison of immune, stromal, and ESTIMATE scores between the high-risk group and low-risk group. B Volcano plots for the enrichment of immune cell types based on the normalized NES score between low-risk and high-risk patients from the GSEA. C The relationship between risk score and the expression of nine ICIs. D Comparison of TIDE scores between the high-risk group and the low-risk group

The relationship between risk scores and patient response to ICI therapy

Over the past decade, ICIs have proven to be promising agents for many solid tumor malignancies, adapting this treatment strategy to an increasingly important role in Neuro-oncology [1, 34]. Therefore, we investigated the potential relationship between risk score and the expression of nine ICIs. Seven ICIs, excluded LAG3 and TIGIT were detected to have higher expression levels in the high-risk group, which may be related to the higher tumor grade of patients in which group (Fig. 7C). The TIDE algorithm then predicted the sensitivity of patients in the high- and low-risk groups to ICI treatment. TIDE scores were significantly higher in the high-risk group than in the low-risk group, suggesting that patients in the low-risk group may benefit from ICI treatment (Fig. 7D).

Expression of AL390755.1 and FLJ16779 mRNA in NHA, U87 and A172 cells

By qRT-PCR validation of the expression of the two lncRNAs, we found that the expression of AL390755.1 was significantly lower in glioma cells than that in NHA cells, while the expression of FLJ16779 showed the opposite results (Fig. 8).

Fig. 8

Expression of AL390755.1 and FLJ16779 in NHA, U87 and A172 Cells, ^∗p < 0.5, ^∗∗p < 0.01, ∗∗∗p < 0.001

Discussion

Glioma is one of the most common, aggressive and challenging malignant tumors in neurosurgery. In this study, we identified an efficient prognostic model consisting of two glycolysis-related lncRNAs using TCGA databases. This model demonstrated strong predictive capacity for glioma patient prognosis as patients in the high-risk group had poorer overall survival compared to those in the low-risk group. Additionally, the nomogram model showed superior performance in predicting prognosis. GO analysis revealed that DEGs associated with risk scores were mainly involved in immune responses and immune cell regulation. KEGG results indicated that risk score-related DEGs were significantly associated with immune disorders and immune response-related pathways. Immune landscape analysis revealed extreme enrichment of immune cells in the high-risk group. Furthermore, patients in the low-risk group appeared to benefit more from ICI treatment. These findings suggest that glycolysis-related lncRNAs could serve as important biomarkers and potenial therapeutic targets for glioma.

Recently, the fields of tumor metabolism and immune oncology have received significant attention. Gliomas cells display higher levels of glycolytic activities compared with normal brain tissue, particularly in glioblastomas (GBM) [5], which are widely acknowledged to be a hallmark of immune cell activation [3]. Zehang Jiang et al. found that there is a positive correlation between glycolytic activity and immune score across all 14 cancer types, with GBM and LGG having the highest correlation [14]. Previous cancer imaging studies have revealed the complex relationship between tumor glucose metabolism rate and immune microenvironment [9, 19]. Reversing this metabolic program could provide adequate glucose for immune cells to stimulate antitumor response. Although checkpoint inhibitors and CAR-T cell therapy are currently the most common immunotherapies for cancer patients, only 20–40% of patients respond to immunotherapy [12, 22], the impact of glucose metabolism shift on immune cell function may be another explanation for this issue. Modulating glioma metabolism represents a logical therapeutic approach. Our study showed that these glycolysis-associated lncRNAs were mainly involved in immune responses based on functional enrichment analysis. Additionally, immune landscape analysis suggested an extreme enrichment of immune cells in the high-risk group, patients in the low-risk group were more likely to benefit from ICI treatment. This observation is consistent with recent findings that increased tumor glycolysis suppressed antitumor immunity [4], and immunotherapy was more favorable in tumors with low glycolysis [17, 27, 31]. Therefore, adjusting cell metabolism along with immunotherapy might be an effective treatment regime. However, large cohorts are required to verify the association between tumor glycolysis and immunotherapy response.

Several studies have indicated that specific metabolism-related genes may serve as prognostic indicators for patients with glioma [6, 10, 25, 26, 32]. In contrast to our results, Wang Jia and colleagues identified six glycolysis-related lncRNAs in glioma patients based on the CGCA database, the best results they obtained had ROC-AUC 0.875 and 0.879 for the training and validation set, respectively [25, 26]. Other models displayed best ROC-AUC ranging from 0.847 to 0.873 [10, 32]. In contrast, our prediction model achieved ROC-AUC 0.943 and 0.909 for the training and validation set, respectively, while the corrected C-index was 0.8705766, suggesting superior accuracy in predicting patient outcomes. Our identification of these two glycolysis-related lncRNAs signatures has substantial clinical significance due to their high sensitivity and specificity as prognostic biomarkers.

Specifically, we have reported a significant finding that FLJ16779 and AL390755.1 were confirmed for the first time to be correlated with the prognosis prediction of glioma. The study of Yongqiang Wang et al. complements these findings by establishing a robust three-lncRNA model for predicting the OS of gastric cancer patients, including OVAAL, FLJ16779 and FAM230D [25, 26], suggesting the crucial role of FLJ16779 in predicting the prognosis of various types of cancers. Little information is available about these glycolysis-related lncRNAs in glioma and other tumors, which prompted us to evaluate the potential function by using GSEA. GO and KEGG analysis revealed that these DEGs were mainly enriched in terms related to immune responses and immune response-related pathways. Consequently, we postulate that the interaction of these two glycolysis-related genes with immune system responses may play a critical role in the carcinogenesis and progression of glioma. Thus, further experimental trials are warranted to validate these preliminary results and address any outstanding questions.

In this study, we have identified a two glycolysis-related lncRNA signature for predicting the prognosis of glioma patients based on TCGA database. This finding provides a promising avenue for exploring new immune checkpoint inhibitor (ICI) therapeutic targets for glioma patients. However, there were certain limitations in this study. It was only a preliminary study exploring the relationship between glycolysis level and immune cell infiltration, and the exact mechanism of action and regulatory relationship need to be further studied.

Materials and methods

Data source

A total of 675 samples from The Cancer Genome Atlas (TCGA)-lower grade gliomas (LGG) and -Glioblastoma (GBM) databases were utilized in this study, of which 5 were normal samples and 670 were glioma samples. Of the glioma samples, 635 samples had survival information and 580 samples had complete clinical information. Glycolysis genes (Additional file 1: Table S1) were downloaded from the Molecular Signatures Database (MSigDB) for all genes in the glycolytic pathway, including BIOCARTA GLYCOLYSIS PATHWAY, GO GLYCOLYTIC PROCESS, HALLMARK GLYCOLYSIS, KEGG GLYCOLYSIS GLUCONEOGENESIS, and REACTOME GLYCOLYSIS [28].

Consistent clustering

The consensus clustering was accomplished in R by using the ConsensusClusterPlus package. Selection of glioma subtypes based on the expression profile of the glycolysis genes employing a k-means clustering approach. The optimal number of clusters was decided by the cumulative distribution function (CDF) curve of the consensus scores and the consensus matrix heatmap. The Kaplan–Meier (K–M) curve was adopted to evaluate the prognosis of the different subtypes and P < 0.05 was considered significant.

Pathway analysis of different subtypes

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway identification was performed by single-sample gene set enrichment analysis (ssGSEA) to assess pathway variations between subtypes. Briefly, ssGSEA calculated an enrichment score (ES) for each KEGG pathway in each sample within each cluster via the GSVA software package. Each ssGSEA ES represented the extent to which the KEGG pathway was up- or down-regulated in the sample (Additional file 2: Table S2).

Differential analysis

The differential expression analysis would be performed in R using the limma package to assess the expression distribution of lncRNAs in different subtypes. The lncRNAs satisfying |log₂ fold change (FC)|> 2 and P < 0.05 were considered as DE-lncRNAs and included in the subsequent analysis. Furthermore, genes satisfying |log₂ FC|> 1 and P < 0.05 were considered to be the risk score-related DEGs between the high- and low-risk groups.

Construction and assessment of the prognostic feature

Here, we worked with a sample of 635 TCGA-glioma subjects with survival information. First, we randomly divided the sample of 635 cases into a training set (n = 445) and an internal validation set (n = 190) in a ratio of 7:3. The training set was used to locate prognostic lncRNAs and the assessment of the predictive validity of the prognostic feature, and the internal validation set was only designed to validate the predictive validity of the prognostic feature. Subsequently, DE-lncRNAs screened by univariate Cox regression analysis (P < 0.05) were further included in multivariate Cox regression analysis with a step function to determine the best variables to use to construct a prognostic feature. The risk score for each patient in the training and internal validation sets was calculated based on the coefficient (coef) value (Additional file 3: Table S3) and expression of each prognostic lncRNA as shown in the formula below:

$${\text{Risk score}} = \frac{e^{{\text{sum(each gene's expression levels}} \times {\text{corresponding coefficient}}})}{e^{{\text{sum (each gene's mean expression levels}}\,\times {\text{corresponding coefficient)}}}}$$

Patients were categorized into high- and low-risk groups based on the median risk score. The K-M analysis for detecting the difference in OS between the high- and low-risk groups was used. To assess and validate the prognostic predictive validity of the prognostic feature, ROC curves were performed in both the training set and the internal validation set.

Creation of the Nomogram

Univariate and multivariate Cox analyses were conducted to explore the independent prognostic factors in patients with glioma. Variables initially included were risk score, age, gender, type of treatment (type), IDH mutation type (IDH), grade, MGMT status (MGMT), and sub-clustering (cluster). Ultimately, independent prognostic factors identified by multivariate Cox analysis (P < 0.05) would be available for the construction of the Nomogram. Additionally, calibration curves were obtained to assess the predictive accuracy of the Nomogram model for patients with 1, 3, and 5-year OS.

GSEA in risk score-related DEGs

Gene Ontology (GO) and KEGG analysis of risk score-related DEGs was implemented based on R software using GSEA in the clusterProfiler package. Here, as functional enrichment studies of lncRNAs are not yet straightforward, we focused on potential functions with risk score-related DEGs to understand the biological processes and pathways that lncRNAs may be involved in. GO analysis was used to reveal the potential biological functions of risk score-related DEGs, while KEGG was responsible for the enrichment of the pathways.

lncRNA-mRNA network

Due to the large size of risk score-related DEGs, we filtered out risk score-related DEGs that were strongly correlated with prognostic lncRNAs (|coefficient (cor)|≥ 0.6 and P < 0.05) using Pearson correlation analysis (Additional file 4: Table S4 and Additional file 5: Table S5). The lncRNA-mRNA network was then mapped and embellished using Cytoscape software.

Immune landscape analysis

The stromal score, immune score, and ESTIMATE score for each sample (n = 635) were calculated by applying the estimate R package. Further Wilcoxon tests were performed to assess the relationship between these scores and risk score (n_high-risk = 317, n_low-risk = 318).

In addition, we performed a GSEA using the clusterProfiler package to identify immune cell types that differed between the high- and low-risk groups. Briefly, genes between the high- and low-risk groups were sorted by |log₂ FC| and then GSEA of immune cells was performed. The normalized ES (NES) was used to identify differences in immune cell types between the two groups.

Cell culture and qRT-PCR

NHA were purchased from Lonza and cultured using an AGM Bullet Kit™ (Lonza, Walkersville, MD) as recommended by the manufacturer. U87 and A172 cells were purchased from ATCC and were cultured in Dulbecco’s modified.

Eagle’s medium (Gibco, Carlsbad, CA) containing 10% fetal bovine serum.

(Gibco, Carlsbad, CA) according to standard protocols. Then, place the cells in a 37 °C, 5% CO2 incubator for culture. Change the medium once a day. The RNA is extracted when the cells grow to 80% confluent. Use TRIzol (ThermoFisher Scientific, USA) to extract total cell RNA. Follow the steps of PrimeScrip reverse transcription kit (Takara, Japan) to reverse transcription into cDNA. Configure the PCR reaction system and analyze it according to the SYBR Premix Ex Taq (Takara, Japan) instruction. Human actin was used as an endogenous control, the relative gene expression was calculated by the 2 ^−ΔΔCt method. Primer sequence utilized for the qRT-PCR analysis is listed in Table 4. Repeat the experiment 3 times independently for each sample.

Table 4 qRT-PCR primer sequence

Statistical analysis

Survival analysis was performed in the survival package. The ggplot2 package was used to plot volcanoes. pROC was used for ROC curve generation. Nomograms were obtained using the rms package. Sensitivity prediction for ICI treatment was achieved through the TIDE online website (http://tide.dfci.harvard.edu/). R software was employed for statistical analysis and P < 0.05 was considered statistically significant.