LGG transcriptomic data [fragments per kilobase of transcript per million mapped reads (FPKM)], copy number variation (CNV), somatic mutation data and data on phenotypic characteristics were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). To re-annotate lncRNA-associated probes in the gene microarray, we downloaded the appropriate lncRNA genomic sequence information from the GENCODE database (GRCh38.gtf, https://www.gencodegenes.org/human/). Furthermore, the sequence information of the microarray probe was used to match with the sequence of lncRNAs to construct an lncRNA expression profile of the re-annotated microarray. Patients with LGGs with missing OS data or survival time less than 30 days were excluded to minimise statistical bias. In addition, we selected three LGG cohorts [CGGA mRNA-seq-693 (sample size: 332, available: 332 ), CGGA mRNA-seq-325 (sample size: 162, available: 332) and GSE16011 (sample size: 80, available: 80)] for external validation. It was worth noting that CGGA mRNA-seq-693 and CGGA mRNA-seq-325 datasets were created by same organization, and the above data sets have similar clinical information publicly available.

We searched the literature related to DNA methylation modifications and selected 20 DNA methylation regulators (Meng et al., 2021) to determine different modification profiles of DNA methylation. These regulators included 14 readers (UHRF1, UHRF2, ZBTB33, MBD4, MBD3, MBD2, MBD1, ZBTB4, ZBTB38, SMUG1, NTHL1, TDG, MECP2 and UNG), 3 erasers (TET3, TET2 and TET1) and 3 writers (DNMT3B, DNMT3A and DNMT1). We screened DNA methylation-related lncRNAs through gene expression correlation analysis and selected 2698 DNA methylation-related lncRNAs. The following parameters were used to screen for DNA methylation-related lncRNAs: |Pearson R| > 0.3 and p < 0.001.

To assess the association with genome instability-related lncRNAs, we combined the expression profiles of these lncRNAs and somatic mutation profiles using a bioinformatic model for tumour genomes derived from a mutation hypothesis. The model calculates the cumulative number of somatic mutations in every sample and ranks patients based on the number of somatic mutations from the highest to lowest (Bao et al., 2020). We designated 25% of patients with the lowest number of mutations and 25% of patients with the highest number of mutations as genomic stability-like (GS) and genomic instability (GU) groups, respectively. The “limma” R package was used to determine differences between groups by comparing the mean expression of DNA methylation-related lncRNAs using the Wilcoxon rank-sum test. Differentially expressed lncRNAs were selected (|log₂ Fc filter| > 0.585 and false discovery rate [FDR]-adjusted p-value < 0.05) and denoted as DNA methylation-related lncRNAs (DMlncRNAs). Subsequently, we normalised the expression of all DMlncRNAs using Z-score analysis. In addition, hierarchical clustering analysis was performed using the “limma”, “pheatmap” and “sparcl” R packages and the lncRNAs were grouped into two clusters by calculating the Euclidean distance. The cluster with low mutation levels was denoted as a GS-like cluster, whereas the one with high mutation levels was denoted as a GU-like cluster (p < 0.05, Mann–Whitney U test). Differential analysis of DMlncRNAs in genomic instability is shown in Supplementary Table S1.

The whole TCGA dataset was divided into the validation and training sets in a ratio of 3:7 to verify the findings (“caret” package). A signature of DMlncRNAs was constructed using the training set and validated using the validation and TCGA sets. The baseline characteristics of patients in the three cohorts are presented in Table 1.

The relationship between the OS of patients and expression levels of DMlncRNAs were examined through univariate Cox proportional risk regression analysis using LGG survival data retrieved from TCGA (p < 0.05). The “glmnet” R package was used for Lasso Cox regression analysis (with 1,000 iterations) to identify DMlncRNAs associated with the OS of patients with LGGs. Finally, the risk coefficients of prognostic DMlncRNAs were determined via multivariate Cox proportional risk regression analysis. A prediction model (DMlncRNA risk score) was constructed based on the expression levels of these prognostic DMlncRNAs and the coefficients evaluated via multivariate regression analysis, using the following formula: DMlncRNA Riskscore = ∑ i = 1 n coefDMlncRNA i × ExprDMlncRNA i In this formula, coef represents the coefficient in multivariate Cox regression analysis, ExprDMlncRNAi represents the expression level of lncRNAs and coefDMlncRNAi represents the coefficients of lncRNAs associated with survival. Patients with LGGs were divided into the low- and high-risk groups according to the median DMlncRNA risk score as the cut-off value. Survival curves were generated for both groups using the Kaplan–Meier method, and logarithmic tests were performed using the ‘survminer’ and ‘survival’ R packages. A p-value < 0.05 was considered statistically significant. In addition, time-dependent receiver operating characteristic (ROC) curves were generated using the ‘survivalROC’ R package to examine the prognostic significance of the risk score. The validation and TCGA cohorts were used to validate the DMlncRNA risk signature.

Clinical data from the training, test and TCGA cohorts and risk scores were used for univariate and multivariate Cox regression analyses to assess the clinical value of the risk score as an independent prognostic marker. The area under the ROC curve (AUC) was determined to evaluate the accuracy of the prognostic signature using the “pROC” R package.

Univariate and multivariate Cox regression analyses were performed for each variable in the training, test and TCGA cohorts using the “survivor” R package to determine whether the DMlncRNA risk score was an independent prognostic marker for other important clinicopathological characteristics. A p-value < 0.05 was considered statistically significant. The reliability of the DMlncRNA risk score in predicting prognosis was determined through clinical stratification analysis. Patients in the TCGA cohort were assigned to subclasses based on the following clinical characteristics: tumour stage (stage II and III), sex (female and male), history of radiation therapy, tumour type (primary and recurrent) and age (≥40 and <40 years). Patients in each clinical subgroup were divided into the low- and high-risk groups based on the median risk score. Survival differences between the high- and low-risk groups in the subgroups were compared using the log-rank test and Kaplan–Meier analysis.

A nomogram was constructed based on the clinical characteristics of DMlncRNAs and the risk score to improve the prognostic value. The “rms” R package was used to generate column line plots for independent prognostic factors and relevant clinical parameters as variables in the training and validation cohorts. Points for each variable were indicated using a horizontal line according to the different variable characteristics. The total number of points for every patient was determined by calculating the sum of points for each value, and the values were normalised to a range of 0–100. The 1-, 3- and 5-years OS of patients with LGGs were calculated by placing them between each prognostic axis and the total score axis. The ‘survcomp’ and ‘rms’ R packages were used to generate calibration plots based on the concordance index (C-index) and its 95% confidence interval (CI). Clinical decision curve analysis (DCA) was conducted using the ‘rmda’ and ‘devtools’ R packages to verify the performance of the column line plots in the validation and training cohorts.

The GSE16011 dataset retrieved from the Gene Expression Omnibus (GEO) database and the mRNAseq-325 and mRNAseq-693 datasets retrieved from the China Glioma Genome Atlas (CGGA) database were used to validate the DMlncRNA signature. The same formula (DMlncRNA risk score) was used to calculate the risk score, and a boxplot was generated to compare gene expression levels based on the following characteristics: chromosome 1p/19q combined deletion status, chemotherapy status, IDH1 mutation status, tumour type (primary and recurrent), tumour grade and age.

Immune, stromal and ESTIMATE scores were evaluated using the ESTIMATE algorithm, and the relationship between the tumour microenvironment and risk scores was examined. The correlation between the OS of patients with LGGs and algorithm scores were evaluated using Kaplan–Meier analysis. The CIBERSORT algorithm was used to determine the proportion of 22 immune cells using the gene expression data of LGG patients to further evaluate differences in immune cell infiltration between the high- and low-risk groups. Only data with a CIBERSORT p value <0.05 was filtered and reserved for the following analysis. The output was directly integrated to generate an entire matrix of immune cell fractions. The levels of immune cells with significantly different proportions in the high- and low-groups were determined using the Wilcoxon rank-sum test. The correlation between the OS of patients with LGGs and levels of the 22 infiltrating immune cells was determined through Kaplan–Meier analysis. In addition, the correlation between immune cells and risk scores was evaluated through Pearson correlation analysis.

The “clusterProfiler” R package was used to perform KEGG and GO pathway enrichment analyses. Statistical significance was indicated by p-value < 0.05.

From June 2020 to June 2021, 16 LGG and adjacent normal brain samples were obtained from eight patients who had undergone surgical dissection and pathological confirmation at the First Affiliated Hospital of Xinjiang Medical University. The present study was approved by the Medical Research Ethics Committee of the First Affiliated Hospital of Xinjiang Medical University. An RNA reagent (Wuhan Servicebio Technology Co., LTD, Wuhan, China) was used to extract total RNA. NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States) was used to determine RNA quantity. mRNA levels were quantified using a two-step reaction process, namely, reverse transcription (RT) polymerase chain reaction (PCR). Servicebio RT First Strand cDNA Synthesis Kit (Wuhan Servicebio Technology Co., LTD, Wuhan, China) was used to synthesise cDNA from RNA. The expression levels of GAPDH, colorectal neoplasia differentially expressed (CRNDE), CYTOR, MPPED2-AS1 and SNHG18 were evaluated through qRT-PCR using the SYBR Green qPCR Master Mix (High ROX) (Servicebio, Wuhan, China). The expression levels were evaluated relative to the GAPDH expression level. The following PCR primer sequences were obtained from Servicebio (Wuhan): GAPDH-F: 5′-GGA​AGC​TTG​TCA​TCA​ATG​GAA​ATC-3′, GAPDH-R: 5′-TGA​TGA​CCC​TTT​TGG​CTC​CC-3′; CRNDE-F: 5′-GGA​AAA​ATC​AAA​GTG​CTC​GAG​TG-3′, CRNDE-R: 5′-ACT​GGC​AAT​CAA​ATA​CAG​CTT​AAC​C-3′; CYTOR-F: 5′-AAA​ATC​ACG​ACT​CAG​CCC​CC-3′, CYTOR-R: 5′-AAT​GGG​AAA​CCG​ACC​AGA​CC-3′; SNHG18-F: 5-GGA​GCC​ACC​CAG​AAA​CTT​AGA​CA-3′, SNHG18-R: 5-CCC​TGG​TGG​ACT​TGA​GTG​GAA-3′ and MPPED2-AS1-F: 5′-TAG​AAA​CAC​CCC​TTC​GGA​AAC​AC-3′, MPPED2-AS1-R: 5′-CCT​TTG​GTG​ACC​TTA​TCT​AGT​TAC​TGA -3′. The expression levels of CRNDE, CYTOR, MPPED2-AS1 and SNHG18 were determined using the 2^(−ΔΔCT) method. The amplification reaction included the following steps: pre-denaturation at 95°C for 10 min, followed by denaturation at 95°C for 40 cycles of 15 s and extension at 60°C for 30 s. Fluorescence signals were recorded from 65 to 95°C at an interval of 0.3°C.

Connectivity Map (cMap), a gene expression profiling database, was used to screen for potential drug compounds against LGGs. We uploaded differentially expressed DMlncRNAs to the database to identify potential connections and bioactive compounds. The linkage score was set between -1 and 1 to assess the closeness of the active agent associated with the query feature. A positive score implied that the drug promoted the expression of high-risk lncRNAs, whereas a negative score implied that the drug inhibited the expression of high-risk lncRNAs. A threshold of p-value < 0.05 was set to indicate significance.

R (version 4.0.3) and Perl tools were used to conduct all statistical analyses. Continuous data were analysed using the Wilcoxon test, and categorical variables were analysed using the Fisher’s exact or chi-square test. Survival differences were estimated using the KM and log-rank tests. In addition, differential expression analyses were performed for LGG (N = 518) and matched normal samples obtained from TCGA and samples obtained from Genotype-Tissue Expression (GTEx) (N = 207) through Gene Expression Profile Interaction Analysis (GEPIA) (http://gepia.cancer-pku.cn/). A p-value < 0.05 indicated statistical significance.

A schematic illustration of the construction of the DMlncRNA prognostic signature and subsequent analyses is presented in Figure 1A. We extracted the TCGA cohort matrix comprising 13,868 lncRNAs and evaluated the expression levels of 20 DNA methylation-related genes. The lncRNAs highly associated with any of the 20 DNA methylation-associated genes were identified as DMlncRNAs (|Pearson R| > 0.3 and p < 0.001). A Sankey plot was generated to visualise the DMlncRNA co-expression network (Figure 1B), and a total of 2,698 DMlncRNAs were identified. Figure 1C shows the correlation between DNA methylation-associated genes and DMlncRNAs in the TCGA set.

The total number of somatic mutations was calculated for all patients, and the patients were ranked in descending order to determine the genomic instability of DMlncRNAs. The bottom 25% (n = 137) and top 25% (n = 133) patients were divided into the GU- and GS-like groups, respectively. A total of 165 differentially expressed lncRNAs were identified as novel DMlncRNAs by comparing the significantly different lncRNA expression profiles of patients in the GS- and GU-like groups. Of the 165 lncRNAs, 85 were downregulated and 80 were upregulated (log FC filter| > 0.585, FDR-adjusted p-value < 0.05; Wilcoxon test) (Supplementary Table S1). A heat map was generated based on the top 20 downregulated and upregulated lncRNAs (Supplementary Figure S1A). A total of 529 samples in the TCGA set were used for unsupervised hierarchical cluster analysis based on the expression levels of the 165 differentially expressed DMlncRNAs, and patients were divided into the GU- and GS-like groups (Supplementary Figure S1B). The findings showed that some lncRNAs expression were significantly different between the two groups.

A total of 468 patients with LGGs in TCGA were divided into the test set (n = 139) and training (n = 329) sets to assess the predictive role of the selected DMlncRNAs in prognosis. Univariate Cox analysis was conducted to examine the association between the expression of 165 DMlncRNAs and the OS of patients in the training set to screen for prognosis-related DMlncRNAs. The results revealed that 132 DMlncRNAs were highly correlated with the prognosis of patients with LGGs (p < 0.001). Lasso regression analysis was conducted on these lncRNAs to minimise overfitting. Lasso regression is a commonly used multiple regression analysis used to fit generalised linear models while performing variable screening and complexity adjustment and enables simultaneous variable selection and regularisation. It is widely used to optimise feature selection with low correlation and prominent predictive values in high-dimensional data. Therefore, it can accurately discriminate among most predictive markers and help to identify prognostic indicators for the effective prediction of clinical outcomes. The first rank value of log λ with a minimum segment likelihood deviation was represented as a dashed vertical line. Subsequently, we identified seven lncRNAs involved in DNA methylation in LGGs (Figure 2A) and determined the optimum value of the threshold via 1,000 rounds of iterations (Figure 2B). Furthermore, stepwise multivariate Cox regression analysis of these seven candidate lncRNAs identified five DMlncRNAs as prognostic risk factors (Figure 2C and Table 2), including CRNDE, AC010273.2, MPPED2-AS1, SNHG18 and CYTOR. We constructed a risk score based on these DMLncRNAs (DMlncRNA risk score) to assess the risk of prognosis of patients with LGGs based on the expression levels of these five independent prognosis-related DMlncRNAs and multifactorial Cox coefficients. The formula for calculating the risk score is as follows: DMlncRNA Riskscore = [CRNDE ∗ 0.2160] + [AC010273.2 ∗ −0.5662] + [MPPED2-AS1 ∗ 0.4249] + [SNHG18 ∗ 0.3898] + [CYTOR ∗ 0.3979].

We first calculated the risk score of all patients in the training cohort, and patients were divided into the low- and high-risk groups using the median risk score as the threshold. Furthermore, we examined the relationship between the expression of the five DMlncRNAs and OS status and risk score in the training and test cohorts (Figures 2D–F), and the heat map indicated that the five DMlncRNAs were significantly upregulated in the high-risk group. Kaplan–Meier analysis revealed significantly better survival outcomes (OS) for low-risk patients than for high-risk patients, indicating that the prognostic risk score was reliable (p < 0.001) (Figure 2G). The survival analysis curves revealed that the 3-, 5- and 10-years survival rates of high-risk patients were approximately 53.6% (95% CI: 44.57–64.5%), 34.62% (95% CI: 24.76–48.4%) and 5.55% (95% CI: 0.98–31.4%), respectively, whereas those of low-risk patients were approximately 91.0% (95% CI: 85.0–97.4%), 77.9% (95% CI: 67.5–90.0%) and 52.3% (95% CI: 35.1–78.0%), respectively. The AUC values of the ROC curves of the DMlncRNA signature for predicting 1-, 2- and 3-years OS were 0.893, 0.919 and 0.866, respectively, (Figure 2H).

A test set including 139 patients was used to validate the prognostic value of the DMlncRNA signature. Using the threshold risk score of the training cohort, the 139 patients in the test set were divided into the high- and low-risk groups. We generated risk curves, scatter plots and heat maps to demonstrate the relationship between risk scores and the OS of patients with LGGs in the test cohort (Figures 3A–C). Similar to the OS of patients in the training set, the OS of high-risk patients in the test set was significantly poorer than that of low-risk patients as demonstrated by Kaplan–Meier curves (p < 0.001) (Figure 3D). Furthermore, the 3-, 5- and 10-years survival rates of high-risk patients were approximately 69.6% (95% CI: 56.66–85.5%), 47.1% (95% CI: 29.31–75.8%), and 23.6% (95% CI: 5.45–98.6%), respectively, whereas the 3-, 5- and 10-years survival rates of low-risk patients were approximately 96.3% (95% CI: 69.21–100%), 84.7% (95% CI: 89.43–100%) and 53.4% (95% CI: 29.69–96.0%), respectively. The 1-, 2- and 3-years ROC curves for the DMlncRNA model showed AUC values of 0.931, 0.892 and 0.930 (Figure 3E). The prognostic performance of DMlncRNAs in TCGA cohort was consistent with that in the test cohort. Patients in TCGA cohort were divided into the low- and high-risk groups, and the results of risk curves, scatter plots and heat maps for the association between risk scores and the survival status of patients with LGGs in TCGA cohort were consistent with those in the training and validation sets (Figures 3F–H). Similarly, the survival curves revealed that the OS of the low-risk group was significantly higher than that of the high-risk group (p < 0.001) (Figure 3I). The 3-, 5- and 10-years survival rates of high-risk patients were approximately 59.9% (95% CI: 52.39–68.4%), 41.0% (95% CI: 32.02–52.4%) and 10.1% (3.27–31.4%), respectively, whereas those of low-risk patients were approximately 93.51% (95% CI: 89.12–98.1%), 79.86% (95% CI: 70.98–89.9%) and 57.28% (95% CI: 43.17–76.0%), respectively. In addition, the AUC values for predicting 1-, 2- and 3-years OS were 0.902, 0.906 and 0.884, respectively, (Figure 3J).

Kaplan–Meier curves were generated to predict the OS of LGGs based on the five DMlncRNAs in the training cohort to evaluate their prognostic value (Figure 4A). The five lncRNAs, namely, AC010273.2, CRNDE, CYTOR, MPPED2-AS1 and SNHG18, were associated with OS (p < 0.001) and identified as risk factors because their upregulation was correlated with a poor prognosis. Similar results were obtained in the validation set (Figure 4B). Subsequently, we constructed a nomogram based on these DMlncRNAs in the training set and included the following factors: age, ATRX mutation status, sex, PTEN mutation status, tumour grade, diagnosis type, tumour staging, IDH1 R132 mutation status, EGFR mutation status, TP53 mutation status and the risk scores of the five DMlncRNAs (Figure 4C). The nomogram was used to predict the 1-, 3- and 5-years survival rates of patients. A C-index of 0.8439 (95% CI: 0.8055–0.8823) was obtained, which verified the predictive value of the nomogram. In addition, a calibration plot (Figure 4D) and clinical decision curve (Figure 4E) demonstrated that the column line plot had an excellent prediction value. The prediction value of the column line plot was higher than that of the risk score. The findings of the validation cohort were consistent with those of the training cohort (Figures 4F–H).

Tumour type (primary and recurrent), tumour grade, sex, age, LGG diagnosis type, TP53 mutation status, IDH1 R132 mutation status, EGFR mutation status, ATRX mutation status and PTEN mutation status were used to evaluate the independence of the DMlncRNA signature in predicting prognosis. In addition, multivariate and univariate Cox regression analyses were conducted. The findings of multivariate analysis indicated that the DMlncRNA risk score was highly correlated with the OS of patients in the training, test and TCGA cohorts (p < 0.05), and LGG diagnosis type and age were significantly correlated with OS (p < 0.05) as well as the DMlncRNA risk score (Table 3 and Figures 5A–F). Furthermore, ROC curves were generated to compare the DMlncRNA risk score with the clinical characteristics of patients. The findings revealed that the risk score was better than tumour type (primary and recurrent), tumour grade, sex, age, LGG diagnosis type, PTEN mutation status, IDH1 R132 mutation status, EGFR mutation status, ATRX mutation status and TP53 mutation status (Figures 5G–I). These findings suggested that the performance of the DMlncRNA risk score in predicting the survival of patients with LGGs was significantly higher.

In addition, we performed a stratified analysis of the DMlncRNA risk score to examine whether its predictive value was dependent on tumour type, tumour grade, radiotherapy, sex and age. Patients in TCGA cohort were divided into the following two groups based on the median age of patients (40 years): a group with patients aged ≥40 years (n = 253) and a group with patients aged <40 years (n = 215). Based on the DMlncRNA risk score, patients in the two groups were further divided into the low- and high-risk groups. OS was significantly different between the low- and high-risk groups in the group with patients aged <40 years (p < 0.001) (Figure 6A) and in the group with patients aged ≥40 years (p < 0.05) (Figure 6B). Furthermore, patients were grouped based on sex, and 212 patients were women and 256 patients were men in TCGA cohort. These patients were further divided into the low- and high-risk groups based on their risk scores. OS was significantly higher in the low-risk group than in the high-risk group (p < 0.001, Figures 6C,D). Furthermore, patients were grouped based on tumour grade, and 241 patients were included in the Grade III group and 226 patients were included in the Grade II group in TCGA cohort. Patients in the Grade II group were divided into the low- and high-risk groups based on the DMlncRNA risk scores, and OS was significantly different in both groups (p < 0.001) (Figure 6E). Similarly, the DMlncRNA risk scores were used to divide patients in the Grade III group into the low- and high-risk groups, and OS was significantly different in both groups (p < 0.001) (Figure 6F). Finally, stratification analyses were performed based on tumour type (primary and recurrent) and whether patients received radiotherapy. The findings indicated that OS was significantly lower in the high-risk group than in the low-risk group, irrespective of whether patients received radiotherapy and tumour type (p < 0.05) (Figures 6G–J). These findings suggested that the DMlncRNA risk score was an independent factor for predicting the OS of patients with LGG.

The CGGA mRNA-seq-693 (n = ) and CGGA mRNA-seq-325 (n = ) datasets were used to validate the prognostic performance of DMlncRNAs. First, we investigated the relationship between the lncRNAs CRNDE, CYTOR and SNHG18 and DNA methylation in the CGGA mRNA-seq-693 cohort. The results indicated that the expression levels CRNDE and CYTOR were significantly correlated (p < 0.05) with tumour grade, age (<40 and ≥40 years), PRS type, 1p19q chromosome codeletion, IDH mutation status and chemotherapy status (Figures 7A–L). Furthermore, SNHG18 was significantly correlated (p < 0.05) with IDH mutation status and age (<40 and ≥40 years), and SNHG18 expression levels were significantly different in the 1p19q chromosome codeletion subgroups (Figures 7M–O). In addition, CRNDE exhibited a good correlation (p < 0.05) with IDH1 (R132) mutation status (Figure 7P) in the GSE16011 dataset.

Similar results were observed in the CGGA mRNA-seq-325 dataset for the lncRNAs CRNDE, CYTOR and SNHG18. the expression levels of these lncRNAs were strongly correlated with chromosome 1p19q codeletion, IDH mutation status, tumour grade, chemotherapy status and age (<40 and ≥40 years) (p < 0.05) (Supplementary Figures 2A–O).

We evaluated immune scores and immune cell infiltration levels in the high- and low-risk groups to examine the role of DMlncRNAs in immune cell infiltration in LGGs. The ESTIMATE (Figure 8A, p = 4.5e-08), immune (Figure 8B, p = 9.6e-07) and stromal (Figure 8C, p = 9.4e-10) scores of the high-risk group were higher than those of the low-risk group. The ESTIMATE, immune and stromal scores of patients with LGGs were significantly correlated with prognosis (Figures 8D–F) (p < 0.05). The survival rate of the high-score group was lower than that of the low-score group, implying that TME features are highly associated with the onset of LGG progression. Furthermore, the association between the expression of DMlncRNAs and immune components was determined via the CIBERSORT algorithm using 22 types of immune cell profiles in the low- and high-risk groups. A total of 149 tissues were found to be eligible for the analysis of CIBERSORT (p < 0.05). In addition, the proportion of tumour-infiltrating immune subtypes was evaluated, and heat maps were generated (Figures 8G,H). Interestingly, the infiltration of M2-type macrophages and monocytes was mainly in all LGG samples. Among the high-risk patients had a higher infiltration of M0-type macrophages. In contrast, low-risk patients seemed to have a higher infiltration of naive B cells. Correlation and differential analyses revealed six TIC types that were strongly correlated with the expression of DMlncRNAs (Figure 8I) (p < 0.05). Finally, immune cells were divided into the high- and low-expression groups based on the median. The results revealed that OS was better in the high-expression group with activated mast cells and activated NK cells than in the low-expression group (Figures 8J,L) (p < 0.05), whereas OS was better in the low-expression group with resting mast cells than in the high-expression group (Figure 8K) (p < 0.05). These results indicated that the expression of DMlncRNAs significantly affected the immune activity in TME.

To further evaluate the biological role of DMlncRNAs in LGGs, we screened for differentially expressed mRNAs in the low- and high-risk groups using the most enriched KEGG pathways and GO terms related to molecular functions (MFs), cellular components (CCs) and biological processes (BPs). The GO terms were mainly associated with extracellular matrix organisation, structural organisation, structural components of the extracellular matrix, sister chromatid separation and chromosome segregation involved in DNA methylation associated with genomic instability (Figures 9A,B). In addition, KEGG analysis revealed that DMlncRNAs were mainly enriched in the Hippo signalling pathway, MAPK signalling pathway, PI3K–Akt signalling axis, p53 signalling pathway, Ras signalling pathway, transcriptional dysregulation in cancer, antigen processing and presentation and autoimmune thyroid disease (Figures 9C,D). These signalling pathways are associated with core biological oncogenic processes, most of which involve regulation of the oncogenic activation pathways and immune checkpoint expression. Therefore, the effects of these pathways on DNA methylation modifications during immunotherapy should be investigated further.

For exploring the potential regulatory relationship of lncRNAs in risk signature, we combined the results of multiple databases to construct a ceRNA network. Firstly, We used starbase online tool to explore the lncRNA-miRNA axis; however, unfortunately, only CRNDE and CYTOR had prediction results for downstream miRNAs. Next, we explored target-mRNAs of miRNAs in three databases (miRDB, miRtarBase, and TargetScan). The above target-mRNAs were overlapped with differentially expressed genes to obtain the final target genes that may be involved in genomic instability. It's worth noting a total of 1,040 differentially expressed mRNAs were identified in the GS- and GU-like groups (log FC filter| > 0.585, FDR-adjusted p-value < 0.05; Wilcoxon test). Finally, we constructed a ceRNA network (316 edges) based on 2 lncRNAs, 112 miRNAs, and 58 mRNAs (Supplementary Figure S3). Interestingly, CRDNE may be closely associated with a higher number of genomic instability-associated mRNAs. In microRNAs, mir-29 family (mir-29a-3p, mir-29b-3p, and mir-29c-3p) had the largest number of edges, which may indicate the important regulatory in ceRNA network.

To verify the expression of DMlncRNAs, three lncRNAs were first analysed using the GEPIA and GTEx databases. The findings indicated that the expression level of lncRNA CRNDE was significantly higher in LGG samples than in normal brain specimens (p < 0.05; Figure 10A). However, the expression levels of CYTOR and SNHG18 were lower in LGG tissues than in normal brain tissues (Figures 10B,C). Furthermore, survival analysis was performed using GEPIA to examine the association between survival rates and the lncRNAs CRNDE, CYTOR and SNHG18. Low expression levels of DMlncRNAs were significantly (p < 0.05, Figures 10D–F) associated with the DFS of patients with LGGs. Finally, to better characterise the expression levels of DMlncRNAs in normal and LGG tissues, eight normal brain tissue samples and eight LGG samples were collected. The expression levels (qRT-PCR) of CRNDE and MPPED2-AS1 were significantly higher in LGG samples than in normal brain tissue samples (p < 0.05, Figures 10G–I, Supplementary Figure S4). However, the expression levels of SNHG18 and CYTOR were higher in normal brain samples than in LGG samples (p < 0.05), which was consistent with the results of GEPIA.

Potential drug compounds significantly associated with the differentially expressed DMlncRNAs were identified using the cMap database. We found that 10 small molecule drugs were negatively associated with LGG, such as calmidazolium (specificity = 0.123), etacrynic acid (specificity = 0.006), monobenzone (specificity = 0.0405), parthenolide amiprilose (specificity = 0.1724), ciclopirox (specificity = 0.0952) and gelsemine (specificity = 0.0409), which may inhibit the development of LGG (Figure 11, Supplementary Table S2). It is worth mentioning that we found that calmidazolium could cross the blood-brain barrier (Lee and Hait, 1985).