N6-Methylandenosine-Related lncRNAs Are Potential Biomarkers for Predicting the Overall Survival of Lower-Grade Glioma Patients

The prognostic value of N6-methylandenosine-related long non-coding RNAs (m6A-related lncRNAs) was investigated in 646 lower-grade glioma (LGG) samples from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) datasets. We implemented Pearson correlation analysis to explore the m6A-related lncRNAs, and then univariate Cox regression analysis was performed to screen their prognostic roles in LGG patients. Twenty-four prognostic m6A-related lncRNAs were identified as prognostic lncRNAs and they were inputted in a least absolute shrinkage and selection operator (LASSO) Cox regression to establish a m6A-related lncRNA prognostic signature (m6A-LPS, including 9 m6A-related prognostic lncRNAs) in the TCGA dataset. Corresponding risk scores of patients were calculated and divided LGG patients into low- and high-risk subgroups by the median value of risk scores in each dataset. The m6A-LPS was validated in the CGGA dataset and it showed a robust prognostic ability in the stratification analysis. Principal component analysis showed that the low- and high-risk subgroups had distinct m6A status. Enrichment analysis indicated that malignancy-associated biological processes, pathways and hallmarks were more common in the high-risk subgroup. Moreover, we constructed a nomogram (based on m6A-LPS, age and World Health Organization grade) that had a strong ability to forecast the overall survival (OS) of the LGG patients in both datasets. We also establish a competing endogenous RNA (ceRNA) network based on seven of the twenty-four m6A-related lncRNAs. Besides, we also detected five m6A-related lncRNA expression levels in 22 clinical samples using quantitative real-time polymerase chain reaction assay.


INTRODUCTION
Gliomas are obstinate intracranial tumors with very poor mortality and recurrence rates (Jiang et al., 2016). They are extremely difficult to remove by neurosurgical resection due to their invasiveness. Lower-grade gliomas (LGGs) can evolve to higher-grade gliomas and develop resistance to chemotherapy. These factors result in gliomas having a persistently high mortality rate. Therefore, searching for therapeutic targets for treating gliomas is urgent (The Cancer Genome Atlas (TCGA) Research Network, 2008;Brat et al., 2015;Jiang et al., 2016).
N6-methylandenosine (m6A) modification, the most abundant epigenetic methylated modification of messenger RNAs (mRNAs) and non-coding RNAs (ncRNAs), plays a vital role in RNA splicing, export, stability and translation Dai et al., 2018). m6A modification is a invertible and dynamical RNA epigenetic process that is regulated by m6A regulators, including "writers" (methyltransferases), "readers" (signal transducers) and "erasers" (demethylases) (Zaccara et al., 2019). Recent research had revealed that m6A modification can regulate oncogenesis and tumor progression in several kinds of cancers, including glioma. For example, METTL14 could enhances leukemia stem cells self-renewing and acts as a cancer promoter in the development and maintenance of acute myeloid leukemia (Weng et al., 2018); knockdown of FTO limits glioblastoma stem cells self-renewing and the proliferation and invasion of lung squamous cell carcinoma (SCC) cells (Cui et al., 2017;Liu et al., 2018); and YTHDF2 restrains cell proliferation by reducing the mRNA stability of EGFR in liver cancer (Zhong et al., 2019). Recently, bioinformatic research has revealed that dysregulation of m6A regulators involved in the malignant development of gliomas (Chai et al., 2019).
Aberrant lncRNA expression is also strongly related to tumor malignancy, and dysregulation of long non-coding RNAs had been confirmed to play a critical role in the pathogenicity of gliomas (Yang G. et al., 2014;Bhan et al., 2017). For instance, HOXA11-AS is reported to be a cell cycle-related lncRNA as well as a biomarker of glioma patients (Wang et al., 2016;Se et al., 2017); MALAT1 was a suppressor of glioma by downregulating MMP2 and devitalizing ERK/MAPK signaling (Han et al., 2016); and knockdown of DANCR could inhibit glioma cell growth and invasion through downregulating miR-135a-5p/BMI1 axis (Feng et al., 2020). However, the full role of m6A regulators in the dysregulation of lncRNAs in cancers remains unclear, and few studies have been performed to explore the mechanisms underlying how m6A modifications contribute to lncRNA-dependent glioma occurrence and development. Thus, understanding how m6A modifications of lncRNAs are involved in glioma progression may help to identify biomarkers that can act as useful therapeutic targets.
Here, based on The Cancer Genome Atlas (TCGA) dataset (n = 476) and the Chinese Glioma Genome Atlas (CGGA) dataset (n = 170), we identified the prognostic significance of m6Arelated lncRNAs by bioinformatic and statistical analysis of data from patients with LGG. Our results showed that 24 m6A-related lncRNAs had prognostic value in both TCGA and CGGA LGG patients. Furthermore, we constructed an m6A-related lncRNA prognostic signature (m6A-LPS) based on the ability of nine m6A-related lncRNAs to predict the OS of LGG patients. In the meanwhile, LGG patients in low-and high-risk subgroups (categorized based on the m6A-LPS) had different prognosis and tumor hallmarks were more common in the high-risk subgroup. Furthermore, an accurate nomogram was constructed to predict OS in patients with LGG and a ceRNA network was built to search the target miRNAs and mRNAs of these m6A-related prognostic lncRNAs.

Datasets and m6A-Related Genes
For TCGA training set, mRNA expression files [Fragments Per Kilobase of transcript per Million mapped reads (FPKM) normalized] were acquired from the Genomic Data Commons Data Portal 1 and the corresponding clinicopathological data was obtained from the cBioPortal website 2 . To obtain a CGGA validation set, RNA-seq data and related clinicopathological data were downloaded from the CGGA website 3 . Gliomas of World Health Organization (WHO) grade II and III were defined as lower-grade gliomas (LGGs), and LGG patients with missing OS values or OS < 30 days were excluded in order to reduce statistical bias in our analysis. Ultimately, we obtained a TCGA training dataset involving 476 patients and a CGGA validation dataset involving 170 patients. In addition, based on previous publications, expression matrixes of 21 m6A-related genes were extracted from the TCGA and CGGA databases, respectively, including expression data on writers (METTL3, METTL14, METTL16, WTAP, VIRMA [KIA1499], RBM15, RBM15B, and ZC3H13), erasers (FTO and ALKBH5) and readers (YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1, and RBMX). The clinicopathological and molecular features for the samples including in current study were collected in Supplementary Table S1.

Annotation of lncRNAs
The long non-coding RNA annotation file of Genome Reference Consortium Human Build 38 (GRCh38) was acquired from the GENCODE website 4 for annotation of the lncRNAs in the TCGA dataset. Based on recognizing the Ensemble IDs of the genes, 14,247 lncRNAs were identified in the TCGA dataset, and 4304 lncRNAs were identified based on the gene symbols in the CGGA dataset. In this part, lncRNAs defined in our analysis included eight types of transcripts (lincRNA, antisense, processed transcript, sense intronic, 3prime overlapping ncRNA, sense overlapping, and macro lncRNA).

Bioinformatic Analysis
Pearson correlation analysis was first implemented to mining m6A-related lncRNAs (with the | Pearson R| > 0.5 and p < 0.001) in each dataset for finding m6A-related lncRNAs. Then univariate Cox regression analysis was implemented to filtrate the prognostic m6A-related lncRNAs in the two datasets, respectively. The m6A-related prognostic lncRNAs screened from the two databases are intersected to obtain the 24 shared m6A-related prognostic lncRNAs. Thereafter, using the R package "glmnet" (Friedman et al., 2010) to conduct least absolute shrinkage and selection operator (LASSO) Cox regression (with the penalty parameter estimated by 10-fold cross-validation) (Tao et al., 2020), we developed a m6A-related lncRNA prognostic signature (m6A-LPS) for the LGG patients involving 9 m6Arelated lncRNAs. The risk score calculating formula is: where Coef i means the coefficients, x i is the FPKM value of each m6A-related lncRNAs.
Risk scores were computed for all patients including in our study. For both two datasets, principal component analysis (PCA) was performed using R programming language (verison3.6.1) (see footnote 5) and scatter diagrams were plotted using the R package "ggplot2." Using the TCGA cohort, differentially expressed genes (DEGs) in the high-risk subgroup in contrast to the low-risk subgroup were identified based on the standards of | log2(Fold change)| > 1 and p < 0.05 using the R package "limma" (Ritchie et al., 2015). Perl programming language was used to perform the prediction analysis of the target miRNAs of the 7 m6A-related lncRNAs in the miRcode database 5 and 59 shared target mRNAs of these miRNAs was found in the miRTarBase 6 , miRDB 7 , and TargetScan 8 databases. Respectively, the 2571 differential expression genes (DEGs) between lowand high-subgroups and the 59 target mRNAs in the ceRNA network were then inputted into the "Metascape" website (Zhou et al., 2019) 9 for functional and pathway enrichment analysis, which involved Canonical Pathways, Reactome Gene Sets, Gene Ontology (GO) Biological Processes and Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG pathway). Additionally, we used GSEA software 10 to investigate the tumor hallmarks that were more common in the high-risk subgroup compared with the low-risk subgroup. The ceRNA network was plotted using the software of "Cytoscape" (Shannon et al., 2003).

Samples and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
We totally collected 22 non-neoplastic and neoplastic samples from patients who underwent surgical treatments in the Neurosurgical Department of The Second Affiliated Hospital of Nanchang University from 2015 to 2019. In these 22 clinical samples, 16 glioma samples were classified into 8 WHO grade II gliomas and 8 WHO grade III gliomas, and 6 non-neoplastic brain tissues (NBTs) were obtained from intractable epilepsy patients. Fresh tumor and non-neoplastic brain tissues were frozen and stored at -80 • C. This research was approved by the Medical Ethics Committee of The Second Affiliated Hospital of Nanchang University and the sample acquisition and usage was performed in accordance with the approved guidelines. Informed consent was acquired from each involved patient.

Statistical Analyses
Kaplan-Meier curves and the log-rank test were used to compare the OS between various subgroups, comprising the low-and high-risk subgroups and additional subgroups based on the expression of each of the nine m6A-related lncRNAs. The student's t-test was used to compare the risk scores (based on the m6A-LPS) between pairs of subgroups based on the following clinicopathological features: isocitrate dehydrogenase (IDH) mutation status (mutant or wildtype IDH), 1p/19q co-deletion status (co-deletion or non-co-deletion), O(6)-methylguanine-DNA methyltransferase (MGMT) promoter methylation status (methylated or unmethylated), age (≤ 40 or > 40 years), WHO grade (II or III) and gender (male or female). Chi-square tests were used to compare the distribution of age, gender, WHO grade, IDH status, and 1p/19q status between low-and high-risk subgroups (partitioned by the median value of risk scores) in LGGs of both datasets (Supplementary Table S2). Univariate and multivariate Cox regression analyses were utilized to evaluate the independent prognostic value of the m6A-LPS regarding OS.
We performed multivariate Cox regression to establish a nomogram, the calibration plots showedthe prognostic predictive accuracy of the nomogram and the C-index was calculated for the nomogram models in both two datasets. These analyses were performed using the R package "rms." The prognostic ability of the nomogram and other predictors (risk score, age, and WHO grade) for 1/3/5-year OS was evaluated by receiver operating characteristic (ROC) curves (R package "timeROC") and the area under the curve (AUC) values (Blanche et al., 2013). The statistical analysiscarried out in this study was using R programming language (version 3.6.1) 11 and SPSS Statistics 11 https://www.r-project.org/ 25 software 12 . The related R codes were uploaded to the git-hub repository 13 .

Identification of m6A-Related lncRNAs in LGG Patients
Firstly, using the downloaded file from the "GENCODE" website, we identified 14247 lncRNAs in the TCGA dataset and 4304 lncRNAs in the CGGA dataset for the following analysis. We then extracted the expression matrixes of 21 m6A-related genes from the TCGA and the CGGA datasets, respectively. A lncRNA whose expression value was correlated with one or more of the 21 m6A-related genes (| Pearson R| > 0.5 and p < 0.001) was defined as a m6A-related lncRNA. Pearson correlation analysis was performed to search m6A-related lncRNAs in each dataset, and we obtained 75 lncRNAs which were significantly correlated with m6A-related genes in both two datasets. Combined with the prognostic information, univariate Cox regression was then implemented to screen m6A-related prognostic lncRNAs from the 75 m6A-related lncRNAs in both the TCGA and CGGA datasets (p < 0.05), respectively. Finally, we found that 24 m6Arelated lncRNAs were significantly correlated with the OS of LGG patients in both two datasets. The work flow was shown in Figure 1A and the correlations between the 24 lncRNAs and the m6A-related genes in the TCGA dataset are shown in Figure 1B. 12 https://www.ibm.com/products/software 13 https://github.com/tzw2019/m6A-LPS-R The results of univariate Cox analysis of the twenty-four m6Arelated lncRNAs are shown in Table 1.

Construction of the m6a-LPS in the TCGA Dataset
To build the m6A-LPS for forecasting the OS of LGG patients, we performed a LASSO Cox analysis on the basis of the 24 m6A-related prognostic lncRNAs in the TCGA cohort and it generated the m6A-LPS which contains 9 m6A-related lncRNAs and coefficient of each (Figures 2A,B). The m6A-LPS involved nine lncRNAs and, for each patient in the TCGA dataset, a risk score was calculated based on the coefficient for each lncRNA ( Figure 2C). Patients in the TCGA cohort were divided into lowand high-risk subgroups based on the median value of risk scores. Kaplan-Meier survival curves depicted that LGG patients with higher risk scores had worse clinical outcomes (lower OS rates and a shorter OS time) (Figure 2D). Risk score and survival status distributions are plotted in Figure 2E. And the ROC curves demonstrated that m6A-LPS harbored a promising ability to predict OS in the TCGA cohort (1-year AUC = 0.898, 3-year AUC = 0.810, 5-year OS = 0.781; Figure 2F).

Validation of the m6a-LPS in the CGGA Dataset
To validate the prognostic ability of m6A-LPS, we calculated risk scores for patients in the CGGA cohort using the same formula.
LGG patients in the CGGA dataset were assigned to low-and high-risk groups based on the median risk score. The results were consistent with the findings in the TCGA dataset: LGG patients with higher risk scores had lower OS rates and a shorter OS time  Color shaded lncRNAs were protective lncRNAs and others were risky lncRNAs.
in the CGGA dataset ( Figure 2G). Risk score and survival status distributions are shown in Figure 2H and it showed that patients with higher risk scores had shorter overall survival time and dead status. The ROC analysis also indicated that m6A-LPS had a strong prognostic value for LGG patients in the CGGA dataset (1-year AUC = 0.808, 3-year AUC = 0.833, 5-year AUC = 0.817; Figure 2I). These results showed that the m6A-LPS had a robust and stable OS-predictive ability.

Stratification Analysis of the m6a-LPS
We attempted to identify whether clinicopathological features were associated with the risk score. The results revealed that LGG patients with wildtype IDH, 1p/19q non-co-deletion status, unmethylated MGMT, older age and WHO grade III (Figures 4A-E) had higher risk scores, while the risk score was not associated with gender ( Figure 4F). To better assess the prognostic ability of the m6A-LPS, we performed a stratification analysis to confirm whether it retains its ability to predict OS in various subgroups. In contrast with patients with lower risk, higher risk LGG patients had worse OS in both the WHO grade II and III subgroups (Figures 4G,H). Likewise, we confirmed that m6A-LPS retained its ability to predict OS for patients aged ≤ 40 or > 40 years (Figures 4I,J) and patients with mutant or wildtype IDH (Figures 4K,L). These data indicated that it could be a potential predictor for LGG patients. Besides, expression levels of the nine m6A-related lncRNAs between primary and recurrent LGGs were compared in the CGGA dataset (all patients in the TCGA cohort were primary LGGs). The result showed that C6orf3 was downregulated while LINC00152 was upregulated in current LGGs and other m6A-related lncRNAs were not differentially expressed between primary and current LGGs (Supplementary Figure S1A).

Principal Component Analysis
The nine prognostic lncRNAs and their correlated m6Arelated genes (m6A regulators) are shown in an alluvial diagram (Supplementary Figure S1B). Most of the correlated m6A regulators were m6A "readers", and no "erasers" were statistically significant associated with the nine prognostic lncRNAs. Based on the expression value of the 21 m6A-related genes, principal component analysis (PCA) was performed to assess the differences between the low-and high-risk subgroups (Supplementary  Figures S1C,D). The results showed that the lowand high-risk patients in both the TCGA and CGGA  datasets were distributed in distinct directions. These results may suggest that differential m6A statuses exist in different risk subgroups.

Pathway and Process Enrichment Analysis and Gene Set Enrichment Analysis (GSEA)
For investigating the potential biological process and pathway involving in the molecular heterogeneity between the low-and high-risk subgroups, we identified 2571 differential expression genes (DEGs) [|log2 (fold change)| > 1 and p < 0.05] between the low-and high-risk subgroups in the TCGA cohort. These DEGs were primarily enriched in these terms: NABA core matrisome, NABA matrisome-associated (Canonical Pathways); GPCR ligand binding, PD-1 signaling, ECM proteoglycans, elastic fiber formation (Reactome Gene Sets); cytokine-mediated signaling pathway, regulation of cell adhesion, blood vessel development and adaptive immune response (GO Biological Processes) (Figures 5A-C). Gene  set enrichment analysis revealed that several tumor hallmarks were enriched in the high-risk subgroup, such as epithelialmesenchymal transition, the reactive oxygen species pathway, P13K-AKT-MTOR signaling, inflammatory response, KRAS signaling, complement, IL2-STAT5 signaling, glycolysis and MTORC1 signaling ( Figure 5D) and so on. These results may give us some insights into the cellular biological effects related to the m6A-LPS.

m6A-LPS Was an Independent Prognostic Factor for LGG Patients
We used univariate and multivariate Cox analyses to assess whether the m6A-LPS was an independent prognostic factor for patients with LGG. Based on the data of LGG patients in the TCGA dataset, univariate Cox analysis indicated that m6A-LPS was remarkably associated with OS [hazard ratio (HR): 3.367, 95% CI: 2.263-5.010, p < 0.001; Figure 6A] and multivariate Cox analysis further showed that m6A-LPS was an independent predictor of OS (HR: 1.720, 95% CI: 1.069-2.770, p = 0.026; Figure 6A). The conclusion was validated in the CGGA dataset, which confirmed that m6A-LPS was an independent predictor of OS for LGG patients in the CGGA validation dataset (univariate: HR: 3.007, 95% CI: 1.920-4.709, p < 0.001; multivariate: HR: 2.056, 95% CI: 1.212-3.487, p < 0.001; Figure 6B). These results indicated that our m6A-LPS, as an independent prognostic indicator, might be useful for clinical prognosis evaluation.

Construction and Validation of the m6A-LPS-Based Nomogram
To create a clinically applicable quantitative tool to predict the OS of LGG patients, we established a nomogram using the risk status (based on m6A-LPS), WHO grade and age in the TCGA dataset and it was also tested in the CGGA dataset ( Figure 6C). Calibration plots showed that the observed vs. predicted rates of 1-, 3-and 5-year OS showed perfect concordance in the TCGA (Figures 6D-F) and CGGA cohorts (Supplementary Figure S2A). Then time-dependent ROC curves were used to assess the prognostic predictive ability of the nomogram and other predictors (risk score, age and WHO grade) in the TCGA (Figures 6G-I) and CGGA datasets (Supplementary Figure S2B and the results revealed that, compared with the other predictors, the nomogram had excellent accuracy regarding 1-, 3-and 5year OS (AUC = 0.899, 0.860, and 0.806, respectively). C-index was also calculated for assessing the predictive ability of the nomogram in both two datasets as and results showed a stable and robust predictive power (C-index for the TCGA dataset: 0.817 and CGGA dataset: 0.642). These data indicated that the nomogram has a robust and stable ability to predictive the OS for LGG patients.

Construction of the ceRNA Network and Functional Enrichment Analysis
To further understand how the m6A-related lncRNAs regulate mRNA expression by sponging miRNAs in LGG, we constructed a ceRNA network based on the m6A-related lncRNAs. Seven of twenty-four lncRNAs were extracted from the miRcode database and 351 pairs of interaction between the seven lncRNAs and twenty-four miRNAs were identified. Then we used three databases (miRTarBase, miRDB, and TargetScan) to search target mRNAs based on the twenty-four miRNAs and totally fifty-nine mRNAs were identified in all the three databases. Ultimately, seven lncRNAs, twenty-four miRNAs and fifty-nine mRNAs were included in our ceRNA network ( Figure 7A). Furthermore, the 59 target mRNAs were used to implemented functional enrichment analysis in the Metascape online tool and we found that these genes were enriched in vasculature development, fibroblast growth factor receptor signaling pathway, negative regulation of growth, mitotic sister chromatid segregation (GO Biological Processes); PID E2F pathway, PID P53 downstream pathway, NABA matrisome associated (Canonical Pathways); pathways in cancer, transcriptional misregulation in cancer, MAPK signaling pathway(KEGG Pathway) (Figures 7B-D). These data may provide us some clues for finding the potential functions of these m6A-related lncRNAs in LGGs.

Validation of the Expression Levels of Five of the m6A-Related lncRNA in Glioma Samples
For validating the expression levels of the m6A-related prognostic lncRNAs in glioma samples, we detected five m6A-related prognostic lncRNAs expression levels in our collected 22 nonneoplastic and neoplastic samples including 16 glioma samples (8 WHO grade II gliomas and 8 WHO grade III gliomas) and 6 non-neoplastic brain tissues (NBTs) by using RT-qPCR assay.
Our results showed that C6orf3, GDNF-AS1, LINC00237, and LINC00925 was downregulated in LGG samples compared with normal brain tissues and LINC00265 was upregulated in WHO grade III gliomas (Supplementary Figure S2C).

DISCUSSION
A total of 646 LGG patients from the TCGA and CGGA datasets were included in our study to exploit the prognostic significance of m6A-related lncRNAs. Twenty-four m6A-related lncRNAs were confirmed to have prognostic value in both the TCGA and CGGA datasets, and nine of them were used to establish an m6A-LPS for predicting the OS of LGG patients. Based on the median risk score, LGG patients were divided into the lowand high-risk subgroups, and the high-risk group had worse clinical outcomes and enrichment of tumor hallmarks and certain malignant related pathways. Multivariate Cox regression analysis showed that m6A-LPS was an independent risk factor for OS. Furthermore, combining m6A-LPS with age and WHO grade, we created a nomogram, and it had a robust ability to predict OS in LGG patients in the TCGA and CGGA datasets. A ceRNA network consist of seven m6A-related lncRNAs, twenty-four miRNAs and fifty-nine mRNAs were established for viewing the potential functions of these m6A-related lncRNAs. Finally, RT-qPCR was used to detected five m6A-related prognostic lncRNA expression levels in total 22 gliomas and normal brain tissues.
Multiple studies have suggested that m6A modification might function as a regulator in cancer pathogenesis, but how it acts in a lncRNA-dependent manner during glioma progression is still unknown. m6A regulators can maintain the malignancy of several kinds of tumors by modifying specific lncRNAs. Negatively controlled by YTHDF3, the lncRNA GAS5 could restrain colorectal cancer progression by causing YAP (Yes1 associated transcriptional regulator) phosphorylation and degradation (Ni et al., 2019). KIAA1429 [also named as Vir like M6A methyltransferase associated (VIRMA)] promotes liver cancer progression via m6A modification of the lncRNA GATA3 (Lan et al., 2019). METTL14 limits malignant progression of colon cancer by inhibiting the lncRNA XIST (Yang X. et al., 2020). Furthermore, lncRNA FOXM1-AS could accelerate the mutual effect of ALKBH5 and FOXM1 nascent transcripts and conduce to the maintenance of glioblastoma stem cells (Zhang et al., 2017). Study had revealed that m6A modification of lncRNAs could influence cancer initiation and progression, and lncRNAs might serve as competing endogenous RNAs, targeting m6A regulators and thereby influencing tumor aggressive progression. Taking the evidence together, we believe that m6A modification is targeted at lncRNAs, and we ought to pay more attention to the interactions and functions of lncRNAs and m6A modifications so as to identify potential prognostic markers or therapeutic targets of cancers.
We identified 24 m6A-related prognostic lncRNAs from 646 LGG patients, and nine of them were included in the m6A-LPS. LINC00152 is highly expressed and acts as an oncogene in colorectal cancer (Ou et al., 2020), and it may enhance malignancy and may be a promising prognostic biomarker of cancers (Seo et al., 2019;Wang et al., 2020). LINC00265 helps ZMIZ2 to stabilize β-catenin by acting as a sponge (and thereby inhibiting several miRNAs) in colorectal cancer (Zhu et al., 2019), and it can be used to predict poor prognosis in acute myeloid leukemia patients (Ma et al., 2018). It has been reported that LINC00665 increases the malignancy of hepatocellular carcinoma by activating the protein kinase R (PKR)/nuclear factor (NF)-κB pathway (Ding et al., 2020) and induce gastric cancer progression by activating the Wnt signaling pathway (Yang G. et al., 2014). Several of the nine lncRNAs were reported to be associated with cancer progression, but there have been few reports regarding glioma, and reports on how the lncRNAs interact with m6A-related genes have been even rarer. Thus, we hope that our results help to identify the prognostic lncRNAs that m6A regulators might target, thereby providing insights into their potential roles in LGG tumorigenesis and progression.
This study included two glioma datasets, the TCGA and CGGA datasets, and our results were derived and validated using them, but there were several limitations in our study. More independent glioma cohorts should be used to validate the identified prognostic m6A-related lncRNAs. Additionally, the roles of the lncRNAs and their interactions with m6Arelated genes should be confirmed using in vitro and in vivo experiments. Our results may provide some clues for further researchs, concentrating on the mechanism process underlying m6A modification of lncRNAs.