Mutation-derived, genomic instability-associated lncRNAs are prognostic markers in gliomas

Data collection

We collected the clinical characteristics, transcription group data, and somatic cell mutation data from GBM and low-grade glioma (LGG) patients. The datasets generated during this study are available in The Cancer Genome Atlas repository (https://portal.gdc.cancer.gov/). These data were matched using sample names. Patients who had no corresponding information about their survival, or who had less than 30 days of treatment were excluded to eliminate interference from non-cancer causes. We differentiated mRNA and lncRNAs using human genome profiles; mRNAs and lncRNAs were annotated using the HUGO Gene Nomenclature Committee (HGNC2) database (https://www.genenames.org/). A total of 629 samples retained paired lncRNA and mRNA expression profiles, and survival information, somatic mutation information, and common clinicopathological features were obtained for further study. All patients with glioma were randomly divided into two groups: training and test sets. A total of 316 patients in the training set were used to identify the prognostic features of lncRNAs and establish a risk model for the outcome. The test set included 313 patients and was used to independently validate the performance of the prognostic risk model. The GSE43378 dataset generated and analyzed during the current study is available in the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE43378) for the external validation of the model.

Technical route

The study process is shown in Fig. 1. We collected and then analyzed the data from somatic cell mutations and transcription groups to obtain genomic instability-related lncRNAs (GIrlncRNAs). The relationship between GIrlncRNAs and mRNA was then analyzed using co-expression analysis. Next, we randomly divided the patient cohort into training and testing sets. Cox and Lasso regression analyses of GIrlncRNAs were conducted to construct a prognostic signature. The signature was evaluated using mutation correlation analysis, model comparison, independent prognostic value analysis, clinical stratification, examination of the external dataset and cell line validation.

Identification of GIrlncRNAs

We used a method from the Mutator Hypothesis (Bao et al., 2020) to identify GIrlncRNAs. Patients with the highest cumulative gene mutation count and who were at the lowest 25% cumulative gene mutation count were assigned to genomic unstable-like (GU) and genome stable-like (GS) groups, respectively. The average expression of lncRNAs between the two groups was compared using the Wilcoxon rank-sum test in the limma package of the R software. The cut-off thresholds were intended to be |fold change| > 2.0 and false discovery rate (FDR) < 0.05.

Construction of GIrLncSig

The “survival” software package of R was used to conduct univariate Cox regression analysis on the training set to assess the relationship between the expression level of GIrlncRNAs and the overall survival time of patients. The least absolute shrinkage and selection operator (LASSO) regression algorithm was used to further screen candidate GIrlncRNAs to construct the GIrlncRNAs prognostic signature (GIrLncSig). The following formula, based on a combination of the Cox coefficient and gene expression, was used to calculate the signature risk score:

$GIrLncSigscore=\sum _{i=1}^{n}coefi\times Ei$

GIrLncSig is the prognostic risk score for patients with glioma. Ei represents the expression level of lncRNAi in patients and coefi represents the coefficient of lncRNAi. The median GIrLncSig score was used as the risk cut-off point to divide glioma patients into low-risk and high-risk groups. The survival curves of the two groups were plotted by Kaplan-Meier method using “Survminer” and “Survival” packages in R language, and a log-rank sum test that obtained P < 0.05 was considered significant.

Real time-PCR validation

The U87, U251, LN229, and U343 cell lines of human glioblastoma and the immortalized cell line SVGp12 were used for the cellular level validation of lncRNA in the model. All cells were cultured in DMEM supplemented with 10% FBS, 100 U/ml penicillin, and 100 U/ml streptomycin. The cultures were maintained at 37 °C in a humidified environment containing 5% carbon dioxide and were confirmed to be mycoplasma-free prior to experimental use.

A total of 1 ml Trizol (Invitrogen, Waltham, MA, USA) was added to the collected cells to extract total cellular RNA, and the absorbance value of RNA at 260 nm was measured using a Nanodrop 2000 UV spectrophotometer. Then, 1 ug of RNA was removed and used to synthesize cDNA by reverse transcription (New England Biolabs, Ipswich, MA, USA), the SYBR Green (Applied Biosystems, Foster City, CA, USA) method and CFX96 real-time PCR system (Bio-Rad, Hercules, CA, USA) were used for real-time polymerase chain reaction (RT-PCR), and actin was used as an internal control. Amplification was performed at 95 °C/120 s followed by 39 cycles at 95 °C/5 s and 60 °C/30 s. Relative expression of RNA was calculated using the 2^−ΔΔCt method. The primers were generated by Sangon Biotechnology (Shanghai, China). The primers for lncRNAs in GIrLncSig are listed in Table 1.

Identified genomic instability-related lncRNAs (GIrlncRNAs) in glioma samples

To identify genomic instability-related lncRNAs, we sorted glioma patients with the number of gene somatic mutation sites and defined the top 25% (n = 227) as the genome unstable (GU) group and the last 25% (n = 110) as the genome stable (GS) group. Next, we compared the expression profiles of lncRNAs in these two groups to identify the lncRNAs that were significantly different. After screening, we identified 91 differentially expressed lncRNAs (Table S1) and selected the top forty significantly different lncRNAs between the GS and GU groups for heatmap plotting (Fig. 2A). These 91 differentially expressed lncRNAs were significantly associated with genomic instability; thus, they were defined as genome instability-related lncRNAs (GIrlncRNAs). Consensus cluster analysis was then performed on 896 samples from The Cancer Genome Atlas (TCGA) collection (including 390 GBMs and 506 LGGs), and all samples were divided into two groups based on the differential expression of GIrlncRNAs and the median number of accumulated somatic mutations (Fig. 2B). The group with the higher number of accumulated somatic mutations was defined as the GU-like group and the group with a lower number of accumulated somatic mutations was defined as the GS-like group. The two clustered groups had significantly different somatic mutation patterns (Fig. 2C). We analyzed the expression correlation between lncRNAs and their target mRNAs. We also selected the top nine strongly correlated mRNAs as target genes since lncRNAs cannot perform direct biological functions but can regulate mRNAs. We then constructed a co-expression network from our results (Fig. 2D).

Construction of GIrlncRNAs signature in the training set

A total of 629 glioma patients from the TCGA project were randomly divided into a training dataset (n = 316) and a test dataset (n = 313). Patients’ clinicopathological characteristics are shown in Table 2. Eighty prognostic related lncRNAs were identified based on univariate Cox proportional risk regression in the training set (Table S2 and Fig. S2). Lasso regression analysis and stepwise multifactor Cox proportional risk regression analysis were performed since lncRNAs have distinct biological functions and many lncRNAs interact with each other. These steps were taken because there were large amounts of data and the potential for inaccurate model construction. A total of 17 lncRNAs were screened as independent prognostic factors to build the prognostic model (Figs. 3A, 3B and Table 3). A prediction model was finally obtained: genomic instability-related lncRNA signature (GIrLncSig) score = (0.0441217761138621 × LINC01579) +(−0.197814767909076×AL022344.1) + (0.0387926268692573 × AC025171.5) + (0.000726291473736049 × LINC01116) + (0.140438109893274 × MIR155HG) + (0.205877895933553 × AC131097.3) + (−0.0369769918481854 × LINC00906) + (0.0588040388194378 × CYTOR) + (−0.0336307964944069 × AC015540.1) + (−0.108512558356276 × SLC25A21.AS1) + (0.00613294332460668 × H19) + (0.0819381463964954 × AL133415.1) + (0.00495065493579931 × SNHG18) + (0.0180933234213193 × FOXD3.AS1) + (−0.0263113038864061 × LINC02593) + (0.0103411968758757 × AL354919.2) + (0.0460814562382004 × CRNDE). A positive coefficient for lncRNAs suggested that high expressions were associated with long survival times as a protective factor in glioma. In contrast, a negative coefficient for lncRNAs indicates that it is a risk factor for glioma. The risk score of each patient in the training set was obtained using GIrLncSig, and then these patients were divided into high- and low-risk groups based on the median risk score (0.0103411968758757). Survival analysis revealed that the low-risk group had significantly better survival rates than the high-risk group (Fig. 3C). The 5-year survival rates were 4.43% in the high-risk group and 14.56% in the low-risk group. The ROC curve showed that the AUC was 0.934, 0.898, and 0.904 at 1, 3, and 5 years, respectively (Fig. 3D). We sorted patients in the training set according to their risk score and there were different expression levels of GIrlncRNAs in the two groups (Fig. 3E). Patients with high-risk scores exhibited an increased expression of risk genes and decreased expression of protective genes, whereas those with low-risk scores exhibited the opposite trend. We also found a significant difference in somatic mutation patterns between patients in the high-risk and low-risk groups, implying that the model could accurately reflect the somatic mutation situation in gliomas (Fig. 3F).

Independent validation of GIrLncSig in the glioma cerebri data set of transcription data and external validation in the GEO data set

An independent TCGA test set of 313 patients was used to examine the performance of GIrLncSig. The test set used the same GIrLncSig and risk thresholds as the training set and the 313 patients in the test set were divided into high-risk (n = 173) and low-risk groups (n = 140). There was a significant difference in overall survival (Fig. 4A). The overall survival rate was significantly lower in the high-risk group than in the low-risk group, which is consistent with the results of the training group. The 5-year survival rate was 7.51% in the high-risk group, which is lower than the 16.4% in the low-risk group (Fig. 4A). ROC curve analysis showed that the AUC of GIrLncSig was 0.875 for 1 year and 0.773 for 3 years (Fig. 4B). The expression levels of GIrLncSig, patient death distribution counts and model lncRNA expression are shown in Fig. 4C. A significant difference was observed in the somatic mutation pattern between patients in the high-risk and low-risk groups (Fig. 4D), and this result was similar in the training group. To further validate the prognostic significance of GIrLncSig, a cross-platform was performed in other independent datasets from different platforms. GSE43378, a dataset from the Gene Expression Matrix data set, was downloaded for further analysis because of the large sample size and complete clinicopathological features. We investigated the relationship between glioma and genomic instability in this independent dataset and found that four lncRNAs (Linc01116, CRNDE, Linc00906, and SNHG18) in GIrLncSig were included in GSE43378. The expressions of Linc01116 and CRNDE were positively correlated with tumor grade (Figs. 4E and 4F) and the survival time was significantly different between their high- and low-expression subgroups (Figs. 4G and 4H). These results were consistent with those observed in the training and test sets.

Evaluation of independent prognostic significance of GIrLncSig and clinical stratification analysis

In the TCGA dataset, the prognosis of GIrLncSig was analyzed by adjusting for clinical stratification, including age (>65 and ≤65 years), gender (male and female), tumor classification (WHO grade II, III and IV), and other clinical factors. In all clinical subgroups, the survival rate in the low-risk group was higher than that in the high-risk group (Fig. 5). This demonstrated that GIrLncSig exhibited a significant independent prognostic prediction value on the overall survival of glioma patients under different clinical stratification conditions.

GIrLncSig was associated with IDH1 mutation status and comparison between GIrLncSig prediction and other lncRNA model predictions

We observed a difference in the IDH1 status in a cohort of glioma patients. IDH1 is a security mutant in the nervous region, as reported, therefore we determined that it was important to evaluate the connection between GIrLncSig and IDH1 mutation status. We compared the differences between the high-risk and low-risk groups on the training set and test set and the results revealed that the proportion of IDH1 mutation was significantly higher in the low-risk group than in the high-risk group, regardless of whether the patient was affected by glioma cerebri (GC) (Fig. 6A) or low-grade glioma (Fig. 6C). These results imply that patients with IDH1 mutations were at lower risk. We then divided the patients into four subgroups (IDH1 mutation/GS-like, IDH1 mutation/GU-like, IDH1 wild-type/GS-like and IDH1 wild-type/GU-like) to consider both the GIrLncSig and IDH1 mutation status. As shown in Figs. 6B and 6D, there were significant differences in the survival time among the four groups. This suggests that IDH1 mutation status might be related to genomic instability and combining them together is better for predicting prognosis in glioma patients.

We then compared the predictability of GIrLncSig and two recently reported lncRNA signatures for survival prediction using the same TCGA glioma patients. Compared to the 10-lncRNA prognostic signature reported by Pan et al. (2020) (PanlncSig), 4-lncRNA prediction signature reported by Li et al. (2019) (LiDlncSig), 5-lncRNA prediction signature reported by Li et al. (2021) (LiXlncSig), 10-lncRNA prediction signature reported by Luan et al. (2019) (LuanlncSig), and 9-lncRNA prediction signature reported by Tao et al. (2021) (TaolncSig) as displayed in Fig. 6E, our GIrLncSig had an AUC value of 0.889 in 1-year OS. These results were more effective than those of PanLncSig (AUC value = 0.852), LiDlncSig (AUC value = 0.835), LiXlncSig (AUC value = 0.790), LuanlncSig (AUC value = 0.842), and TaolncSig (AUC value = 0.862) in predicting patient survival. In summary, the above results confirmed the reliability and effectiveness of GIrLncSig in predicting GC patients.

Validation of GirLncSig

To further validate the prognostic significance of GirLncSig in gliomas, we verified eight lncRNAs in GirLncSig in four glioblastoma cell lines (U87, U251, LN229, and U343) and one normal cell line (SVGp12). The expression levels of these eight lncRNAs were significantly higher in glioblastoma cell lines than in normal cells (Fig. 7A). Similarly, the results of clinical samples also proved that the expression of GirLncSig in normal brain tissue (N1–N3) was lower than that in glioma tumor tissue (T1–T6), further validating the prognostic significance of GirLncSig.

We then examined the relationship between GirLncSig and cell proliferation and found that in the glioma samples from the TCGA database, the expression of almost all model-related sub-risk lncRNAs was positively correlated with the expression of cell proliferation markers Ki67 and PCNA (Figs. 8A and 8B). Similarly, we also detected CDH2 and VIM, which are markers of cell migration, invasion, and mesenchymal transition. The eight lncRNAs were found to be positively correlated (Figs. 8C and 8D). These findings tentatively show that GIrLncSig can promote tumors by affecting cell proliferation and EMT.