Identification and validation of a novel eight mutant-derived long non-coding RNAs signature as a prognostic biomarker for genome instability in low-grade glioma

Identification of lncRNAs associated to genomic instability in Low-grade glioma patients

The flow chart or our work was shown in Figure 1. Frequency of the cumulative number of somatic mutations was calculated for each patient and then sorted in descending order to determine LncRNAs associated with genomic instability. Clinical and pathological characteristics of TCGA-LGG patients were used to explore the clinical outcome of these patients and are presented in Table 1. Patients were divided into (GU) and (GS) groups based on the cumulative number of somatic mutations. lncRNA expression profiles of 137 patients in the GS group and 133 patients in the GU group were compared, and lncRNAs with significantly different expression levels were identified. Analysis using SAM method showed that 59 lncRNAs were differentially expressed (|logFC |>1, FDR adj p-value < 0.05, Supplementary Table 1). Unsupervised hierarchical clustering analysis was performed using the 59 differentially expressed lncRNAs for the 529 samples retrieved from TCGA. The 529 samples were divided into two groups based on the levels of the 59 differentially expressed lncRNAs. Analysis showed that the two groups had significantly different somatic mutation patterns (Figure 2A). The GU group showed significantly higher cumulative somatic mutations compared with the GS group. Notably, the GU group showed significantly higher median number of cumulative somatic mutations compared with the GS group (P < 0.001) (Figure 2B). In addition, the expression level of UBQLN4 gene [20] (a genomic instability-associated gene overexpressed in aggressive tumors) was compared in the GU and GS groups (Figure 2C).

The top 10 most relevant mRNAs among the differential lncRNAs in the GS and GU groups were screened out for each lncRNA as its target genes. This dataset of lncRNA-related mRNAs was used to construct a lncRNA-mRNA co-expression network (Figure 2D). GO analysis of LncRNA-related target genes showed that the mRNAs in this network were highly involved in formation and progression of genomic instability, including mitosis, maintenance of synaptic machinery, transporter complexes and activity of various ion channels (Figure 2E and 2F). KEGG pathway analysis of LncRNA-associated target genes identified several signalling pathways associated with genomic instability, tumorigenesis, progression and treatment of low-grade glioma including B cell receptor signalling pathway, T cell receptor signalling pathway, Rap1 signalling pathway, JAK-STAT signalling pathway, PD-L1 expression and PD-1 checkpoint signalling pathway, the Wnt signalling pathway, cAMP signalling pathway and cGMP-PKG signalling pathway (Figure 2G and 2H).

Analysis showed that 59 lncRNAs are implicated in genomic instability. Variations in their expressions may destabilize the cellular genome and disrupt cellular homeostasis of regulatory networks associated with lncRNAs thus inducing various signalling pathways that promote cancer development. The 59 differentially expressed lncRNAs were then referred as genomic instability- associated lncRNAs (GILncRNAs).

Identification of a genomic instability-mutant lncRNA signature for prognosis using the training set

To further explore the predictive prognosis role of the candidate lncRNAs associated with genomic instability, 477 patients with low-grade glioma from the TCGA project were divided into a training group (n = 240) and a testing group (n = 237). Expression levels of 59 genomic instability-associated lncRNAs were analyzed using univariate Cox proportional risk regression to identify prognosis associated lncRNAs. Analysis showed 54 genomic instability-associated lncRNAs significantly correlated with prognosis of low-grade gliomas in the training set (P < 0.001). Lasso regression analysis was performed using these lncRNAs to avoid over-fitting and 16 lncRNAs associated with genomic instability in LGG were identified (Figure 3A). Notably, the optimal value of the penalty parameter was determined by performing 1000 replicates of cross-validation (Figure 3B). Further, the 16 candidate LncRNAs were screened to identify those with independent prognostic value. A total of eight genomic instability-associated LncRNAs including H19, FLG-AS1, AC091932.1, AC064875.1, AL138767.3, AC010273.2, AC131097.4, and ISX-AS1 were identified as independent prognostic risk factors using stepwise regression multi-factor Cox regression analysis identified (Table 2). A GILncSig was constructed using the coefficients of multifactorial cox analysis and expression levels of the eight independent prognosis-associated LncRNAs for assessing the prognostic risk in low-grade glioma patients. the score was calculated using the formula: GILncSig Risk score = [H19^*0.0293] + [FLG-AS1^*–0.5959] + [AC091932.1^*0.3356] + [AC064875.1^*0.3690] + [AL138767.3^*–0.6290] + [AC010273.2^*0.6908] + [AC131097.4^*0.7026] + [ ISX-AS1^*–0.5280].

The risk score for an individual invalid in the training set was obtained using GILncSig and LGG patients were then grouped into low- and high-risk groups using the median risk score as the threshold. Kaplan-Meier analysis showed that overall survival of patients in the low-risk group was significantly higher compared with that of patients in the high-risk group. This finding indicates that the prognostic biomarker for predicting risk scores was accurate (P < 0.001; Gehan-Breslow-Wilcoxon test: P = 4.0207E-9) (Figure 3C). Survival analysis curves showed that the three-year survival rate among LGG patients with high-risk LGG was approximately 56.8% at 95% CI confidence interval of [46.41%–69.6%], and the three-year survival rate for patients in the low-risk LGG group was approximately 94.9% at 95% CI confidence interval of [89.3%–100%]. Furthermore, the five-year survival rate for patients with high-risk LGG was approximately 39.0% at 95% CI confidence interval of [27.2%–55.9%], and the five-year survival rate for patients with low-risk LGG was approximately 78.6%, at 95% CI confidence interval [66.3%–93.3%]. The ten-year survival rate for patients with high-risk LGG was approximately 10.7%, at 95% CI confidence interval [2.25%–50.8%], whereas the low-risk LGG patients had a ten-year survival rate of approximately 60.2% at 95% CI confidence interval of [41.8%–86.6%]. ROC curve analysis of GILncSig over time for 1, 2, and 3 years showed an area under the curve of 0.921, 0.942, and 0.945, respectively (Figure 3D).

Patients in the training set were grouped based on the score and expression level of GILncSig was determined to represent the somatic mutation counts of the patients. A risk heat map, mutation scatters plot and gene expression map were generated to show the relationship between the risk score and gene expression levels of each LGG sample (Figure 3E). The heat map of the expression profiles of the 8 LncRNAs showed that lncRNA FLG-AS1, AL138767.3, and ISX-AS1 were significantly highly expressed in the low-risk group compared with the levels in the high-risk group. On the other hand, lncRNA H19, AC091932.1, AC064875.1, AC010273.2, and AC131097.4 were significantly highly expressed in the high-risk group compared with the low-risk group. Comparative analysis showed significant differences in somatic mutation patterns and UBQLN4 gene expression patterns between invalids in the low- and high-risk groups.

High-risk group patients showed a higher number of somatic mutations compared with the patients in the low-risk group (P < 0.001, Figure 3F). Interestingly, UBQLN4 expression level was significantly higher in low-risk patients compared with the level in high-risk patients (P = 0.008, Figure 3G).

Validation of a genomic instability-mutant lncRNA signature for prognosis using the test set and TCGA set

The GILncSig model was validated using a testing set of 237 patients to test its prognostic performance. Using the same GILncSig and risk thresholds as the training set, the 237 patients in the testing set were grouped into a low-risk group (n = 120) and high-risk group (n = 117). Kaplan-Meier curves showed that the high-risk group sample had significantly lower OS compared with the low-risk group (p < 0.001; Gehan-Breslow-Wilcoxon test: P = 2.0904E-9) (Figure 4A). Survival analysis curves showed that the 3- and 5-year survival rates of LGG patients in the high-risk group were approximately 53.8% and 31.1%, at 95% CI confidence intervals of [42.42%–68.3%] and [18.65%–51.8%], respectively. The 3-, 5-, and 10-year survival rates of LGG patients in the low-risk LGG group were approximately 94.1%, 81.7%, and 48.4% at 95% CI confidence intervals of [88.49%–100%], [69.97%–95.4%] and [29.00%–80.6%], respectively. Analysis of the 1-year, 2-year, and 3-year ROC curves for GILncSig in the test group over time showed AUC values of 0.920, 0.909, and 0.886, respectively (Figure 4B). Expression of GILncSig and somatic mutation counts in the test samples are presented in Figure 4C. Significant differences in somatic mutation counts were observed between patients in the low- and high-risk groups (p < 0.001, Figure 4D). The high-risk group showed significantly lower expression levels of UBQLN4 compared with the low-risk group (p < 0.001, Figure 4E).

GILncSig's prognosis in the TCGA group was similar to the findings of the test group. Patients in the TCGA group were divided into low- and high-risk groups. Overall survival analysis showed that the overall survival time of patients in the high-risk group was significantly less compared with that of the low-risk group (p < 0.001, Figure 4F). The survival rates at three, five and ten years for patients with high-risk LGG were approximately 56.65%, 36.74%, and 11.64%, at 95% CI confidence intervals of [48.85%–65.7%], [27.57%–49.0%], and [4.21%–32.2%], respectively. The survival rates at 3, 5, and 10 years for patients with low-risk LGG were approximately 95.5%, at 95% CI confidence intervals of [91.94%–99.2%], [73.98%–91.2%], and [39.74%–73.8%], respectively. The 1-year, 2-year, and 3-year ROC curve analysis of the TCGA group over time showed AUC values at 0.919, 0.920, and 0.913, respectively (Figure 4G). GILncSig expression levels and somatic mutation count in TCGA samples are presented in Figure 4H. Analysis showed a significant difference in somatic mutation pattern between patients in the high- and low-risk groups (p < 0.001, Figure 4I). The high-risk group showed significantly lower expression levels of UBQLN4 compared with those for the low-risk group (p < 0.001, Figure 4J).

Construction of 8 genomic instability-mutant lncRNAs signature and correlation with clinical features for Lower-grade glioma

Correlation between overall survival and the eight genomic instability-associated lncRNA models was determined using Kaplan-Meier curves using the training set (Figure 5A). Analysis showed that lncRNA H19, AC091932.1, AC064875.1, AC010273.2, and AC131097.4 were negatively correlated with overall survival, indicating that they are risk factors for low-grade glioma, therefore, high expression levels of the lncRNAs are associated with poor prognosis. On the contrary, lncRNA FLG-AS1, AL138767.3, and ISX-AS1 were positively correlated with overall survival (p < 0.001). These findings show that lncRNA FLG-AS1, AL138767.3, and ISX-AS1 play a protective role in low-grade glioma. Similar findings were obtained through validation using the validation set (Figure 5B). A nomogram was constructed based on tumor grade, gender, age and the eight genomic instability-related lncRNA risk scores in the training group using multivariate Cox regression results (Figure 5C). Analysis of the nomogram, showed that the prognosis-related factors can be used to accurately predict 1-, 3- and 5-year survival of patients. The consistency index (C-index) was used for validation to evaluate the predictive value of the nomogram. Analysis showed that the C-index of the nomogram was 0.8693 and the 95% CI confidence interval of the C-index was [0.8356–0.9030]. The calibration curve (Figure 5D) and clinical decision curve (Figure 5E) showed that the nomogram had a good predictive effect.

Predictive performance of the nomogram was significantly higher compared to that of the risk score model. Similar results were obtained using the validation set (Figure 5F–5H).

The 8 genomic instability-associated lncRNAs were correlated with infiltration of immune cell subtypes in LGG

To explore the relationship between the expression levels of eight genomic instability-associated lncRNA models and the tumor microenvironment, the abundance of tumor-infiltrating immune subpopulations was determined using the CiberSort algorithm. Furthermore, the correlation between the risk scores of the eight lncRNA-constructed models and infiltration of immune cell subtypes in LGG was analyzed using the training set. The correlation coefficients between macrophages M0, macrophages M1, memory resting CD4 T cells, and CD8 T cells and risk scores were 0.23, 0.26, 0.27, and 0.22, respectively (P < 0.001, Figure 6A). These findings indicate that infiltration of tumor immune cells is positively correlated with the risk scores of the eight genomic instability-associated lncRNAs. Activated mast cells (R = –0.21, p < 0.05) and monocytes (R = –0.29, p < 0.001) were negatively correlated with risk scores of the eight lncRNAs. In addition, verification using the validation set showed that macrophages M1 (R = 0.31, p < 0.001) and memory resting CD4 T cells (R = 0.29, p < 0.001) were positively correlated with the risk scores of the eight genomic instability-associated lncRNAs. On the other hand, analysis using the validation set showed that activated mast cells (R = –0.27, p < 0.001) and monocytes (R = –0.2, p < 0.05) immune cell subtypes were negatively correlated with risk scores of the eight genomic instability-associated lncRNAs (Figure 6B). The above correlations between GILncSig signature and the infiltrating immune cells were confirmed with xCell platform in both the training set and validation set (Figure 7A–7B).

Verification of the genomic instability-mutant lncRNA signature model genes using two external independent LGG datasets

Cross-platform validation of genomic instability-associated lncRNA models using independent datasets from different platforms showed that only four of the eight lncRNAs in GILncSig (AC064875.1, AC131097.4, FLG-AS1, H19) retrieved from CGGA mRNA-seq-693 represented a large sample size and were matched clinicopathological features. Therefore, we explored the relationship between lncRNA AC064875.1, AC131097.4, FLG-AS1, H19 and low-grade glioma and genomic instability using the CGGA mRNA-seq-693 dataset. Analysis showed that expression levels of AC064875.1 and H19 were significantly correlated with age (< = 41 and >41 years), tumor grade, IDH1 mutation status and chromosome 1p19q joint deletion (p < 0.05), however, expression levels of these lncRNAs were not correlated with gender and MGMT methylation status (Figure 8A and 8B). AC131097.4 expression levels were significantly different between age subgroups (p < 0.05, Figure 8C). Moreover, FLG-AS1 was significantly correlated with gender and tumor grade (especially between grade II vs IV, III vs IV) (p < 0.05, Figure 8D).

Similar findings were obtained using the GSE16011 dataset for lncRNA H19. Expression level of H19 was significantly correlated with tumor grade and IDH1(R132) mutation status (p < 0.05, Figure 8E). Survival analysis was performed using lncRNA in the model obtained from CGGA mRNA-seq-693 and GSE16011 datasets to determine their correlation with prognosis. In the CGGA mRNA-seq-693 datasets, the OS of patients with high expression of AC064875.1 and H19 was significantly lower compared with that of patients with low expression levels of AC064875.1 and H19 (p < 0.001, Supplementary Figure 1A–1B). On the contrary, high expression level of FLG-AS1 was correlated with a higher overall survival compared with low expression levels of FLG-AS1 (p < 0.05, Supplementary Figure 1C). In GSE16011 dataset, patients with low H19 expression level showed significantly higher overall survival compared with the OS of patients with high expression level of H19 (p < 0.001, Supplementary Figure 1D). Notably, these findings were consistent with the findings obtained using the TCGA training group and test group.

Genomic instability-associated lncRNA signature performance compared with existing lncRNA-related signatures in survival prediction

The predictive performance of the genomic instability-associated lncRNA model was compared with three glioma-associated lncRNA signatures reported in previous studies (Table 3). The three signatures included a 6-lncRNA signature reported by Lin's study [21] (hereafter referred to as LinlncSig) and the 8-lncRNA signature from Li's study [22] (hereafter referred to as LilncSig), which used the same cohort of TCGA patients. The third study was a 3-lncRNA signature reported by Qiu's study [23] (hereafter referred to as QiulncSig) which used 167 TCGA-LGG patients with radiological response information. The AUC of GILncSig for one-year OS was 0.919, which was significantly higher compared with the AUC for LinlncSig (AUC = 0.854), LilncSig (AUC = 0.796) and QiulncSig (AUC = 0.833) (Figure 9A). The AUC of GILncSig for three years OS was 0.913 which was significantly higher compared with the AUC for LinlncSig (AUC = 0.775), LilncSig (AUC = 0.769) and QiulncSig (AUC = 0.760) (Figure 9B). Further, the AUC of GILncSig for five years OS was 0.851 which was significantly higher compared with the AUC for LinlncSig (AUC = 0.699), LilncSig (AUC = 0.644) and QiulncSig (AUC = 0.757) (Figure 9C). We turned the GILncSig, LinLncSig, LiLncSig and QiuLncSig signature into the dichotomous variables in order to combine the signatures, which was named as CombinedSig signature. As the ROC curves showed (Supplementary Figure 2A–2C), we found the CombinedSig signature could increase the prognostic power for 1 year and 3 year survival rate compared to other four independent signatures. Although the combined signature could not increase the predictive power for 5 year rate, the AUC values of the GILncSig signature we constructed in this study were still higher than the other three recently published lncRNA prognostic signature as well as CombinedSig signature, indicating that the GILncSig signature had better prognostic performance in predicting survival in LGG patients. We also compared the GILncSig signature with the other clinical factors including Age, Gender, IDH status and Grade. The results showed that whether the prediction of 1 year, 3 year and 5 year, the power of GILncSig signature for predicting prognosis was better than the Age, Gender, IDH status and Grade factors (Supplementary Figure 2D–2F). These findings indicate that GILncSig had significantly higher prognostic performance for survival prediction compared with the three recently published lncRNA markers. The sensitivity and specificity of each GILncSig cut-off value for predicting 5-year survival were shown in Supplementary Table 2.

Identification of the independence of genomic instability-associated lncRNAs signature in prognostic prediction

A multivariate Cox regression analysis for age, gender, tumor grade, tumor type (primary and recurrent), LGG diagnostic types and GILncSig risk score was performed to see if the GILncSig signature was an independent prognostic predictor. Multivariate analysis showed that GILncSig was significantly correlated with OS in each group in the training set, validation set and TCGA set (p < 0.001). In the multivariate analysis, two other clinical factors, namely age and diagnostic types were significant (p < 0.05) in addition to GILncSig (Table 4 and Figure 10A–10D).

Therefore, we performed a stratified analysis to determine whether GILncSig prognostic value was independent of age and tumor grade. Patients in the TCGA set were grouped as a young < = 41 years group (n = 246) and an age > 41 years group (n = 231) based on the median age (41 years). Each age group of patients was further classified as a low- or high-risk group based on GILncSig. The OS of low- and high-risk groups in the young < = 41 years group was significantly different (p < 0.001, Figure 11A). In addition, the OS in the age > 41 years group was statistically significant for the low- and high-risk groups (p < 0.001) (Figure 11B). All patients with low-grade glioma were then stratified by tumor grade, and patients in the TCGA set were stratified into Grade II (n = 231) and Grade III groups (n = 246). Patients with pathological grade of Grade II were further grouped into high- and low- groups based on GILncSig score. Analysis showed significant differences in OS between the high- and low- groups (p < 0.001, Figure 11C). Furthermore, GILncSig was used to classify pathologically graded Class III patients into high- and low-risk groups, and analysis showed significant difference in OS between the two groups (p < 0.001, Figure 11D). These findings indicate that GILncSig is an independent prognostic factor in predicting the overall survival in patients with low-grade gliomas.

Identification of the relationship between genomic instability-associated lncRNAs signature and IDH1 somatic mutations

Chi-square test analysis showed that a significantly higher proportion of patients in the low-risk group possessed IDH1 mutant phenotype compared with the proportion in the high-risk group of the training set, testing set, and TCGA set. In the training set, 90% of patients in the low-risk group had an IDH1 mutation, which was significantly higher compared with the proportion in the high-risk group (63%) (p < 0.001, Figure 12A). In the testing set, 92% of patients in the low-risk group had the IDH1 mutation, which was significantly higher compared with the proportion in the high-risk group (63%) (p < 0.001, Figure 12B). In the TCGA group, a significantly lower number of patients in the high-risk group (63%) had the IDH1 mutation compared with the number of patients with the IDH1 mutation in the low-risk group (91%) (p < 0.001, Figure 12C). These findings indicate that GILncSig is associated with IDH1 mutation status and can be used as a mutational marker for the IDH1 gene. The IDH1 gene mutant is associated with high genomic stability, whereas an increase in the number of wild-type IDH1 genes induces genomic instability, impairs non-homologous end-joining DNA repair, and increases susceptibility to DNA damage [24, 25]. This implies that IDH1 wild type effect on genomic instability can be used to improve the clinical prognosis of low-grade gliomas.

Previous studies report that IDH1 mutations are linked to more prolonged survival and may serve as an independent prognosis biomarker for low-grade gliomas [26, 27]. Therefore, we compared the prognosis performance of GILncSig on patient outcome with prediction of outcome using the IDH1 mutation status. The log-ranch test was used to group the patients into four groups, namely IDH1 Mutation/GS group, IDH1 Mutation/GU group, IDH1 Wild/GS group, and IDH1 Wild/GU group (p < 0.001). Use of GILncSig for patients with IDH1 wild-type and IDH mutation showed longer overall survival in the IDH1 Mutation/GU group compared with that for the IDH1 Wild/GU group (Figure 12D). Interestingly, the IDH1 Wild/GS group showed a higher OS rate and had a better outcome compared with that of the IDH1 Mutation/GS group. These findings indicate that GILncSig combined with IDH1 mutation status has greater prognostic significance compared with use of IDH1 mutation status alone.

Acquisition and processing of low-grade glioma data

Low-grade glioma gene expression profile data, clinical features and somatic mutation information were obtained from TCGA database (https://portal.gdc.cancer.gov), a cancer genome mapping database, and VarScan2, a genomic variant-based investigation platform (GenomeVIP) [46, 47]. VarScan2 is used for analysis of somatic mutations and copy number alterations (CNA) in exome data using tumor-normal tissue pairs. TCGA-LGG cohort comprised 529 low-grade glioma samples, which we identified using the Ensembl database and GENCODE Release 29 (GRCh38.p12) (https://www.gencodegenes.org/human/) annotations were used to discriminate the pattern of low-grade glioma lncRNA copy number alterations and expression [48]. We postulated that most GENCODE lncRNAs (97%) are located in the covered regions. All low-grade glioma samples obtained were randomly divided into training and testing groups. A total of 477 patients were included after excluding patients with incomplete clinical information and survival time less than 30 days. The training set comprised 240 patients, which were used for analysis of lncRNA signatures and generation of a prognostic risk model. A test set of 237 cases was used to independent validation of the performance of the prognostic risk model. Two additional sets for independent verifications of glioma, GSE16011 and CGGA mRNA-seq-693 were acquired from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16011) and Chinese Glioma Genome Atlas (http://www.cgga.org.cn/download.jsp), respectively, characterized by large sample sizes and multiple clinicopathological features. All the data we used were quoted from the TCGA and GEO databases, and had obtained the consent of the data submitter. And the patient data in this work were also acquired from the publicly available datasets whose informed consent of patients were complete. All the authors declared that there was no conflict of interest in this study.

Identification of lncRNAs associated with genome instability in LGG

lncRNA somatic mutation profiles and expression profiles in low-grade glioma genomes were retrieved for identification of mutation-derived binding genomic instability-associated lncRNAs. The proportions of somatic mutations for each LGG individual were calculated then patients were sorted in descending order based on the number of somatic mutation frequencies using Wilcoxon rank-sum test. Further, gene mutation Unsupervised clustering was performed to define the genomic unstable group (GU) which comprised the first 25% of patients and the genomic stable group (GS) was defined as the last 25% of patients. The lncRNAs expression profiles of the GU group and GS group were compared by SAM method (|logFC|>1 and FDR-adjusted with p < 0.05) and low-expressed genes were excluded (The mean value of lncRNA expression in all samples was less than 0.5) to determine genomic instability-related lncRNAs.

Identification of genome instability-associated lncRNAs prognostic signature for low-grade glioma

A univariate Cox proportional risk regression profiling was used to determine the correlation among expression levels of genomic instability-associated lncRNAs and overall survival (p < 0.001) of patients with low-grade gliomas. The cox regression analysis was downscaled using the Lasso analysis to further select the best lncRNAs for construction of risk models. A total of 1000 cross-validations were performed to prevent overfitting before identifying the genomic instability-related lncRNAs associated with OS. Multivariate cox regression analysis was then performed to determine the risk coefficients of prognostic markers for lncRNAs associated with genomic instability.

Genomic instability-associated lncRNAs signature (GILncSig) was constructed based on the multiple regression analysis of coefficients and expression levels of lncRNAs associated with genomic instability for prognosis prediction using the following equation:

$$\text{GILncSig Riskscore}={\displaystyle \sum _{i=1}^n{\text{coefGILncRNA}}_i\times \text{Exp}}{\text{ GILncRNA}}_i$$

where, GILncSig Risk score represents the prognostic risk score of patients with low-grade glioma. Coef represents multivariate cox regression coefficient, GILncRNA_i represents the ith genomic instability-associated lncRNA, whereas Exp GILncRNA_i represents the expression level of genomic instability-associated lncRNAs. Further, univariate and multivariate Cox regression analysis were used to determine whether genomic instability-associated lncRNAs are independent predictors for low-grade gliomas.

Identification and verification of genome instability-associated lncRNAs survival analysis for Low-grade glioma

Patients were grouped into GILncSig low- and high- risk groups using the median scores of LGG patients in training set as risk cutoff values. Kaplan-Meier survival curve analysis was performed to compare differences in overall survival (OS) between the low- and high- groups. Hazard ratios (HRs) and 95% confidence interval (CI) were used to analyze whether the prognoses of the two groups were significantly different. A time-dependent subject receiver operating characteristic curve (ROC) analysis was utilized to compare the specificity and sensitivity of GILncSig risk score for prognosis of patients with low-grade glioma. The above results were then validated using a validation set.

Construction of a nomogram based on genome instability-associated lncRNAs

Genomic instability-associated lncRNA risk scores and clinical indices were incorporated into the construction of the model to optimize its predictive capability. Nomogram plots were generated to predict patient survival at 1-, 3- and 5- years based on the multivariate Cox regression results. Calibration curves were used to verify its predictive value. The C index and 95% CI confidence interval were evaluated to determine the prediction performance of the nomogram using clinical decision curve analysis (DCA). These results were then validated using the validation set.

Association of genome instability-associated lncRNAs signature with the immune infiltration

The relationship between lncRNA models and immune cell infiltration was evaluated using expression profiles using ESTIMATE and CIBERSORT algorithms to perform correlation with levels of immune cell subpopulations in the samples [49, 50]. Unsupervised hierarchical clustering of the samples was performed. The infiltrative expression of 22 immune cells in the LGG training set was determined using CiberSort R package. Further, the correlation among immune cell infiltration and risk score was determined using Spearman test. The results were then validated with a validation set to explore the potential mechanisms of different risk groups reflected by GILncSig.

Functional enrichment analysis

Pearson correlation analysis was performed to establish the correlation between lncRNA and mRNA paired expression, using the top 10 mRNAs that co-expressed lncRNA-associated chaperone proteins. Further, to explore the potential function of lncRNAs, functional enrichment analysis of lncRNA-related mRNA chaperone proteins was performed. The significantly enriched gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were then identified.

External validation of genes in the genome instability-associated lncRNAs signature model

Validation of the signature model genes was performed using CGGA mRNAseq-693 tool in the CGGA database and the GSE16011 dataset from GEO database. The equation used to calculate risk scores was used for validation. Kaplan-Meier curve analysis was performed, and expression of genes of different ages, genders, grades and IDH1 mutation status were presented in boxplots.

Statistical analysis

All statistical analyzes and plots were performed or generated using the R software version 4.0.3. We used the Pearson correlation analysis to determine the relationship between the two. "Limma" R package was used for differential gene screening. The "glmnet" package was used for LASSO Cox regression model analysis. The "survivor" and "survminer" packages were used for perform survival data analysis and for generation of Kaplan-Meier plots. A two-tailed log-rank test was used to determine statistical significance of survival curves. ROC analysis was performed using the "TimeROC" package. The "rms" package was used to construct the Nomogram and for generation of calibration curves. The "clusterProfler" package was used for GO and KEGG enrichment analysis. P < 0.05 was statistically significant. We used non-parametric tests to deal with non-normal distributions, and used parametric tests to deal with normal distributions.