An autophagy‐related long non‐coding RNA signature for glioma

Glioma is one of the most common types of malignant primary central nervous system tumor, and prognosis for this disease is poor. As autophagic drugs have been reported to induce glioma cell death, we investigated the potential prognostic role of autophagy‐associated long non‐coding RNA (lncRNA) in glioma patients. In this study, we obtained 879 lncRNAs and 216 autophagy genes from the Chinese Glioma Genome Atlas microarray, and found that 402 lncRNAs are correlated with the autophagy genes. Subsequently, 10 autophagy‐associated lncRNAs with prognostic value (PCBP1‐AS1, TP53TG1, DHRS4‐AS1, ZNF674‐AS1, GABPB1‐AS1, DDX11‐AS1, SBF2‐AS1, MIR4453HG, MAPKAPK5‐AS1 and COX10‐AS1) were identified in glioma patients using multivariate Cox regression analyses. A prognostic signature was then established based on these prognostic lncRNAs, dividing patients into low‐risk and high‐risk groups. The overall survival time was shorter in the high‐risk group than that in the low‐risk group [hazard ratio (HR) = 5.307, 95% CI: 4.195–8.305; P < 0.0001]. Gene set enrichment analysis revealed that the gene sets were significantly enriched in cancer‐related pathways, including interleukin (IL) 6/Janus kinase/signal transducer and activator of transcription (STAT) 3 signaling, tumor necrosis factor α signaling via nuclear factor κB, IL2/STAT5 signaling, the p53 pathway and the KRAS signaling pathway. The Cancer Genome Atlas dataset was used to validate that high‐risk patients have worse survival outcomes than low‐risk patients (HR = 1.544, 95% CI: 1.110–2.231; P = 0.031). In summary, our signature of 10 autophagy‐related lncRNAs has prognostic potential for glioma, and these autophagy‐related lncRNAs may play a key role in glioma biology.

Glioma is one of the most common types of malignant primary central nervous system tumor, and prognosis for this disease is poor. As autophagic drugs have been reported to induce glioma cell death, we investigated the potential prognostic role of autophagy-associated long non-coding RNA (lncRNA) in glioma patients. In this study, we obtained 879 lncRNAs and 216 autophagy genes from the Chinese Glioma Genome Atlas microarray, and found that 402 lncRNAs are correlated with the autophagy genes. Subsequently, 10 autophagy-associated lncRNAs with prognostic value (PCBP1-AS1, TP53TG1, DHRS4-AS1, ZNF674-AS1, GABPB1-AS1, DDX11-AS1, SBF2-AS1, MIR4453HG, MAPKAPK5-AS1 and COX10-AS1) were identified in glioma patients using multivariate Cox regression analyses. A prognostic signature was then established based on these prognostic lncRNAs, dividing patients into low-risk and high-risk groups. The overall survival time was shorter in the high-risk group than that in the low-risk group [hazard ratio (HR) = 5.307, 95% CI: 4.195-8.305; P < 0.0001]. Gene set enrichment analysis revealed that the gene sets were significantly enriched in cancer-related pathways, including interleukin (IL) 6/Janus kinase/signal transducer and activator of transcription (STAT) 3 signaling, tumor necrosis factor a signaling via nuclear factor jB, IL2/ STAT5 signaling, the p53 pathway and the KRAS signaling pathway. The Cancer Genome Atlas dataset was used to validate that high-risk patients have worse survival outcomes than low-risk patients (HR = 1.544, 95% CI: 1.110-2.231; P = 0.031). In summary, our signature of 10 autophagyrelated lncRNAs has prognostic potential for glioma, and these autophagyrelated lncRNAs may play a key role in glioma biology.
Glioma is one of the most common types of malignant primary central nervous system tumors with poor prognosis, comprising approximately 44% of central nervous system tumors [1]. The prognosis of glioblastoma (GBM) is the worst among gliomas, in which the median overall survival (OS) for patients with GBM is 15-23 months and the 5-year survival rate is less than 6% [2,3]. Upon diagnosis, the standard treatment of glioma includes maximal surgical resection, chemotherapy, such as temozolomide, and  outcome is unfavorable due to the complex genetic mechanism.
Autophagy is the physiological process that directs degradation of proteins and whole organelles in cells. The activation of autophagy is divided into normal and pathological conditions. Under normal circumstances, autophagy represents a response to several stresses by providing the necessary circulating metabolic substrates for survival. In addition, autophagy is active in some pathological processes in order to maintain cellular homeostasis, such as neurodegenerative diseases, pathogenic inflammation, aging and cancer [7]. In recent years, many studies have sought to find new potential targeted therapies by investigating autophagy pathways [8][9][10]. In addition, autophagic drugs induce cell autophagic death (type II cell death) and cause glioma cell death. Whether this is an alternative and emerging concept for the study of novel glioma therapies remains largely unknown [11].
Long non-coding RNAs (lncRNAs) have a wide range of functional activities [12]. They play a significant role in physiological processes, including RNA decay, genetic regulation of gene expression, RNA splicing, microRNA (miRNA) regulation and protein folding [13]. lncRNA regulates many proteins that are important for autophagy. Impaired functioning of lncRNAs participates in glioma pathogenesis, such as cellular apoptosis and proliferation [14]. HOX (homeobox) transcript antisense RNA (HOTAIR) is a lncRNA that plays an important role in the regulation of cancer transformation, mainly due to extensive miRNA-HOTAIR interactions and its effect on matrix metalloproteinases (MMPs) [15]. There may be an involvement of HOTAIR-interacting miR-NAs and MMPs in autophagy regulation [16][17][18]. Sufficient evidence shows that lncRNAs mediate transcriptional and post-transcriptional levels of autophagy-related genes to regulate the autophagy regulatory network [19,20]. This paper proposes to construct a coexpression network of autophagyrelated lncRNAs using bioinformatics methods, providing a theoretical basis for the treatment of gliomas [21].
Therefore, autophagy-related lncRNAs may have potential value in the prognosis of glioma patients and may serve as potential therapeutic targets. Here, we aimed to establish an autophagy-related lncRNA signature in glioma and to advance the targeted treatment of glioma.

Materials and methods
Information extraction of glioma patients   dataset. TCGA GBM dataset (FPKM level 3) was included in our analysis as a validation dataset with 160 GBM patients.

LncRNA and autophagy gene screening
The profiles of lncRNAs and autophagy genes were obtained from the CGGA ALL mRNAseq dataset.
Specifically, the autophagy gene list was obtained from the Human Autophagy Database (HADb, http://autophagy. lu/clustering/index.html). All of the mRNA expression data were normalized by log2 transformation. Pearson correlation was applied to calculate the correlation between the lncRNAs and autophagy-related genes. A lncRNA with a correlation coefficient |R 2 | > 0.3 and P < 0.05 was considered to be an autophagy-related lncRNA.

Signature development
First, univariate and multivariate Cox regression analyses were performed to evaluate the prognostic value of autophagy-related lncRNAs. The lncRNAs with a P-value < 0.01 by univariate analysis were included in the multivariate stepwise regression Cox analysis to establish the risk score. We used the previous report to determine the risk score for each patient using the following formula: Risk score = b gene1 9 expr gene1 + b gene2 9 expr gene2 + ÁÁÁ + b genen 9 expr genen . Cox analysis was performed to build a signature for predicting survival. For more detail, we assigned risk scores by a linear combination of the expression levels of lncRNAs weighted by regression coefficients (b). The b value was calculated by log transformation of the hazard ratio (HR) from the multivariate Cox regression analysis. High-risk and low-risk groups were established based on the median risk score. The lncRNA expression is defined as expr genen .

Gene set enrichment analysis
Gene set enrichment analysis (GSEA) was used to interpret gene expression data. This method derives its function by analyzing gene sets, so it can be used to determine whether the gene set shows a statistically significant difference between the two biological states. In this study, we verified whether genes that are differentially expressed between two groups are enriched during autophagy.

Statistical analysis
The expression levels of autophagy-related lncRNAs were elevated (P ≤ 0.05). Construction of the autophagy-lncRNA coexpression network was completed using CY-TOSCAPE software [22] (version 3.4.0; The Cytoscape Consortium, San Diego, CA, USA). Pearson correlation analysis and Cox regression analysis were performed using SPSS STATISTICS software (version 24; IBM Corp., Armonk, NY, USA). Survival status was the basis for univariate cox regression analysis. PRISM 7 (GraphPad Software Inc., La Jolla, CA, USA) was used to generate Kaplan-Meier curves. GSEA (http://www.broadinstitute.org/gsea/index. jsp) was used to distinguish between two sets of functional annotations. Statistical significance was set at a threshold of a two-tailed P < 0.05.

Construction of a coexpression network for autophagy-lncRNAs
We identified a total of 878 lncRNAs in the CGGA dataset, which was extracted from the CGGA database. A total of 215 autophagy-related genes were extracted from the Human Autophagy Database (HADb, http:// autophagy.lu/clustering/index.html). We constructed an autophagy-lncRNA coexpression network to identify autophagy-related lncRNAs. Finally, 402 lncRNAs were identified (|R 2 | > 0.3 and P ≤ 0.05).

Identification of a signature of 10 autophagyrelated lncRNAs in patients with glioma
First, we identified autophagy-related lncRNAs by constructing autophagy-lncRNA coexpression networks (P ≤ 0.05). In addition, we used univariate Cox regression analysis based on 402 autophagy-associated lncRNAs to screen prognostic genes. We ranked the prognostic autophagy-related lncRNAs in ascending order by their P values. We used a P value of 0.05 as the cutoff value, and the lncRNAs that satisfied this were used for signature development. Our training set was a collection of 325 glioma patients from the CGGA dataset. A total of 19 lncRNAs have prognostic value for glioma patients (P < 0.01). Subsequently, 10 autophagy-related lncRNAs were found to be independent prognostic factors for glioma patients (Table 1 and Fig. 1), of which five lncRNAs were unfavorable factors (TP53TG1, ZNF674-AS1, COX10-AS1, DDX11-AS1 and SBF2-AS1) and five lncRNAs were confirmed to be favorable prognostic factors for glioma (PCBP1-AS1, DHRS4-AS1, GABPB1-AS1, MAPKAPK5-AS1 and MIR4453HG) ( Table 2 and Fig. 2).
The prognostic impact of an autophagy-related lncRNA signature for glioma Next, we use a risk score method to develop an autophagy-related lncRNA signature. We divided the glioma patients into two groups (low-risk group and high-risk group) by median risk score (Fig. 3). As a result, the risk score could significantly predict the OS of glioma patients, in which the OS period is longer in the lowrisk group than that in the high-risk group (median OS 1211 vs 346 days; log rank P < 0.05). Additionally, the Cox regression analysis also revealed a significant prognostic effect of the risk score on the glioma patients (HR = 5.307, 95% CI: 4.195-8.305; P < 0.0001, Fig. 4). Further, we also explored whether the risk score signature is an independent predictor for the prognosis of glioma patients by multivariate Cox regression analysis. As a consequence, a HR of 2.736 indicated that the risk score could significantly contribute to the prediction of survival of glioma patients, eliminating the influence of other factors such as sex, age, grade, radiotherapy, chemotherapy and the molecular status (IDHDampR, TP53.1, EGFR, ATRX and EZH2) ( Table 3).

Clinical value of the lncRNA signature for glioma patients
Subsequently, we also determined the clinical value of the 10-lncRNA signature regarding the grade, radiotherapy and chemotherapy. As shown in the Table 4, the risk score tends to increase in the higher grades, suggesting that this lncRNA signature might be associated with the progression of glioma. Interestingly, the risk score was lower in patients receiving radiotherapy than that in patients without radiotherapy (t = À2.267, P = 0.025). In contrast to the results of radiotherapy, a higher risk score was found in patients without chemotherapy, while the patients who had received chemotherapy presented a lower risk score. Moreover, we also assessed differences in risk score based on molecular status. As a result, lower risk scores were found in those with the IDH mutation than in those without, indicating a potential association between the lncRNA signature and IDH mutation.

Validation in the TCGA dataset
Next, these results were further validated in the additional dataset (TCGA) using the same b value. In total, 160 GBM patients were enrolled for the validation of the lncRNA signature (Fig. 5). We divided these patients into the high-risk and low-risk groups on the basis of the median value of the risk score. Consistent with the results derived from the CGGA dataset, the high-risk patients had a shorter median OS than that of the low-risk patients (median OS 385 vs 468 days; log rank P = 0.012; Fig. 6). This finding was further validated by Cox regression analysis, in which the high-risk group tended to have a shorter OS time for GBM patients than that of the low-risk group (HR = 1.544, 95% CI: 1.110-2.231; P = 0.031). In light of these results, we could confirm that the lncRNA signature provides a robust prediction for the prognosis of glioma patients.

Gene set enrichment analysis
Further functional annotation was conducted through GSEA. The results revealed that the differentially expressed genes between the two groups were enriched in the autophagy-related and tumor-related pathways. As result, a total of 19 gene sets were significantly enriched at a nominal P-value < 5% (Table 5). Among the gene sets, several pathways are well-established in cancers, including interleukin (IL) 6/Janus kinase/signal transducer and activator of transcription (STAT) 3 signaling, tumor necrosis factor a signaling via nuclear factor-jB, IL2/STAT5 signaling, the p53 pathway and the KRAS signaling pathway (Fig. 7). Moreover, the gene sets were also found to be involved in the vital functions of tumorigenesis and progression of cancer. For instance, epithelial mesenchymal transition, angiogenesis and hypoxia were closely related to the invasion and metastasis of cancer (Fig. 8). Notably, the GSEA revealed that the gene sets were involved in the reactive oxygen species pathway, interferon (IFN)-c response, IFN-a response and inflammatory response, which are strongly associated with autophagy (Fig. 9). Taken together, the defined autophagy-related genes contribute to vital cancer and autophagy pathways, which might provide strong evidence for a cancer-targeted treatment for glioma.

Discussion
Glioma is the most aggressive and common type of primary brain tumor in humans. With the development of clinical management of glioma, some prognostic factors are well characterized, including tumor size, tumor grade and stage. High-throughput biological technologies are being widely used to predict cancer recurrence and tumor metastasis by detecting the alteration of miRNAs or genes [23,24]. The major class of lncRNAs, as a complement to genes or miRNAs, provides a promising opportunity to predict the risk of recurrence of glioma [25]. However, so far, there has been no systematic process to identify lncRNA signature sets for predicting the survival of glioma patients. Therefore, it is necessary to establish a lncRNA signature to predict the prognosis of glioma patients.
In this study, two datasets (CGGA and TCGA) were collected to explore the prognosis of autophagyrelated lncRNAs for glioma patients. In the first step, we identified 402 lncRNAs through the lncRNA-autophagy gene co-expression network. Furthermore, we identified 10 autophagy-associated lncRNA signatures that could divide glioma patients into high-and lowrisk groups based on the median risk score. Additionally, it was found that the OS is longer in the low-risk group than that in the high-risk group. Through univariate and multivariate Cox regression analyses, we can conclude that the signature is an independent factor that is significantly related to OS.
Although little is known about the role of autophagy in cancer therapy to date, recent studies suggest that autophagy therapy will become a new approach to glioma treatment [26,27]. In recent studies, IFN-c was found to influence autophagy and cell growth in human hepatocellular carcinoma (HCC) cells. IFN-c is a cytokine with anti-viral and immune regulation. The cytokine induces autophagosome formation and transformation of microtubule-associated protein 1 light chain 3 proteins and can inhibit cell growth and nonapoptotic cell death in Huh7 cells. In addition, autophagy in Huh7 cells is also activated by the overexpression of interferon-regulatory factor-1. Eventually, induced autophagy will inhibit IFN-c and cell death in HCC [28]. Since autophagy can respond to a variety of stresses to promote the survival of cancer cells, it has protumorigenic functions. Glucose metabolism promotes adhesion-independent conversion driven by oncogene insult-mutationally active Ras. In human cancer cell lines carrying KRAS mutations and cells ectopically expressing oncogenic H-Ras, autophagy is induced after the extracellular matrix is isolated. If autophagy is inhibited by RNA interference-mediated depletion of multiple autophagic regulators or genetic deletion, Ras-mediated conversion and glycolytic capacity proliferation independent of adhesion will be impaired. In addition, when the availability of glucose is decreased, the conversion and proliferation of autophagy-deficient cells expressing oncogenic Ras are unaffected, which is just the opposite of that in autophagy-competent cells. In conclusion, autophagy can promote the unique mechanism of Ras-driven tumor growth in specific metabolic environments [29].
Among the 10 autophagy-related lncRNAs, PCBP1-AS1, DHRS4-AS1, MAPKAPK5-AS1 and GABPB1-AS1 were risk-associated genes, while TP53TG1, ZNF674-AS1, DDX11-AS1, SBF2-AS1, MIR4453HG and COX10-AS1 were protective genes. Specifically, we also found that the high-risk group was enriched in the glycolysis pathway. Consistent with our studies, a recent study revealed that TP53TG1 might affect the expression of glucose metabolism-related genes under glucose deprivation, leading to cell proliferation and migration of glioma cells [30]. Additionally, MAPKAPK5-AS1 regulates gene expression by acting with miRNAs and is significantly associated with the OS of liver cancer [31]. Furthermore, the expression of COX10-AS1 in oral cavity and oropharyngeal squamous cell carcinoma is more than twice that of normal cells [32].
All of the lncRNAs we identified directly or indirectly regulate autophagy, many by regulating miR-NAs; thus, we must perform lncRNA-mRNA coexpression analyses to assess the function of lncRNAs [33][34][35]. Therefore, we can conclude that due to the various functions of lncRNAs, the 10 autophagyrelated lncRNAs we identified will be potential therapeutic targets [12,36].
In conclusion, by constructing an autophagy-lncRNA coexpression network, we identified a signature of 10 autophagy-related lncRNAs, which has prognostic value for glioma patients. In addition, our study classified low-risk and high-risk groups based on the median risk score, and each showed different autophagy states.