Comprehensive Analysis of Autophagy-Associated lncRNAs Reveal Potential Prognostic Prediction in Pancreatic Cancer

Pancreatic cancer (PC) is a highly malignant tumor in the digestive system. Both long noncoding RNAs (lncRNAs) and autophagy play vital roles in the development and progress of PC. Here, we constructed a prognostic risk score system based on the expression profile of autophagy-associated lncRNAs for prognostic prediction in PC patients. Firstly, we extracted the expression profile of lncRNA and clinical information from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases. The autophagy-associated genes were from The Human Autophagy Database. Through Cox regression and survival analysis, we screened out seven autophagy-associated lncRNAs and built the risk score system in which the patients with PC were distinguished into high- and low-risk groups in both training and validation datasets. PCA plot displayed distinct discrimination, and risk score system displayed independently predictive value for PC patient survival time by multivariate Cox regression. Then, we built a lncRNA and mRNA co-expression network via Cytoscape and Sankey diagram. Finally, we analyzed the function of lncRNAs in high- and low-risk groups by gene set enrichment analysis (GSEA). The results showed that autophagy and metabolism might make significant effects on PC patients of low-risk groups. Taken together, our study provides a new insight to understand the role of autophagy-associated lncRNAs and finds novel therapeutic and prognostic targets in PC.


INTRODUCTION
Pancreatic cancer (PC) is a poorly prognostic malignant tumor. Its incidence and mortality rank second and fifth in digestive system tumors in the United States and China separately, of which the five-year survival rate is about 9% (1,2). Current treatments cannot significantly improve the prognosis of PC patients; meanwhile, the development of pancreatic cancer treatment is relatively slow compared to other tumors, so surgery still represents the most effective treatment to cure resectable pancreatic cancer. Due to the lack of diagnosis at the early stage and the highly malignant characteristics of PC, PC patients frequently exhibit lymph node metastasis and local invasion when the diagnosis is made, leading to approximately 80% of patients losing surgical chances (3). Therefore, the exploration of more effective innovative targets for pancreatic cancer is urgent and necessary.
Autophagy is the homeostatic mechanism through a membrane-mediated process that delivers cytoplasmic organelles and proteins to lysosomes for degradation. There is growing evidence that the level of autophagy can be responded by intracellular and extracellular stresses, such as ER stress, oxidative stress, hypoxia, nutrient shortage, etc., thereby involving tumor progression (4). In pancreatic cancer, autophagy plays a significant tumorigenic role in keeping cancer cell survival and promoting metabolism (5). Hydroxychloroquine (HCQ), which can inhibit autophagy combined with Gemcitabine, is currently being tested in many clinical trials. Consequently, researching new biomarkers related to autophagy to improve early diagnosis and assess prognosis is a promising avenue for PC patients.
Long noncoding RNAs (lncRNAs) mostly have no proteincoding potential of which transcripts are longer than 200 nucleotides. lncRNA can affect different functions at an epigenetic, transcriptional, and post-transcriptional level, and play a vital role in regulating cancer cell behaviors and autophagy (6). There is a recent study indicated that downregulated lncRNA LINC00160 suppressed autophagy and drug resistance in hepatocellular carcinoma by regulating miR-132-targeted PIK3R3 (7). Zhang et al. elaborated lncRNA PVT1 induced cytoprotective autophagy and promoted growth via sponging to miR-20a-5p and regulating ULK1 both in vitro and in vivo in pancreatic ductal adenocarcinoma (8). Another study also demonstrated that silencing lncRNA LINC00160 facilitated autophagy and apoptosis of pancreatic cancer cells (9). Considering that several lncRNAs may influence cancer behaviors through mediating autophagy, it is crucial to explore autophagy-associated lncRNAs to predict the prognosis of PC patients.
In our current study, we analyze the relationship between autophagy-associated lncRNA profiles and clinical information in 178 PC patients from The Cancer Genome Atlas (TCGA) database. The survival analysis showed that seven lncRNAs (AC245041.2, LINC02257, AC006504.8, AC012306.2, AC125494.2, FLVCR1-DT, and AC005332.6) were prognostic biomarkers for patients with PC. Then, the seven lncRNAs were used to develop a risk score system after Cox regression analysis.
Finally, we constructed a prognostic signature that can be applied to independently predict the prognostic of PC patients. These candidate autophagy-associated lncRNAs may become the potential prognostic prediction for PC.

Sample Datasets
The RNA-seq data and clinical information of PC patients were downloaded from the TCGA data portal (https://portal.gdc. cancer.gov/) and ICGC (https://icgc.org) databases, respectively. Then, we transformed the RNA sequence data to the lncRNAs and mRNA (protein coding) based on the annotated gene IDs in the Ensembl project. Because the data were extracted from the public database, there was no requirement for ethics committee approval.

Identification of Autophagy-Related lncRNAs
The autophagy gene list was obtained from The Human Autophagy Database (http://www.autophagy.lu/index.html), employing Pearson correlation analysis to screen the relationship between the lncRNAs and autophagy-related genes. An absolute value of correlation coefficient > 0.4 (|R|>0.4) and P < 0.05 were considered statistically significant. Based on the above standard, the autophagyrelated lncRNAs were filtrated for subsequent analysis.

Survival Analysis
Kaplan-Meier (K-M) survival analyses of the autophagyassociated lncRNA were performed using the survival package in R. The patients were classified into high expression and low expression groups using optimal cut-off values determined by the survminer R package (Version:0.4.3). Log-rank P < 0.05 was considered statistical significant.

Construction of Co-Expression Network and Function Analysis
To better understand the relation between lncRNAs and mRNAs, the lncRNA-mRNA co-expression network was visualized by Cytoscape software (http://www.cytoscape.org/). To investigate the functions of these lncRNAs, the co-expression of mRNAs was analyzed by gene ontology (GO) terms enrichment including biological process, molecular function, and cellular component. A P value of < 0.05 was statistically significant.

Construction of the Risk Score System
Firstly, we used the univariate Cox regression analysis to confirm prognostic autophagy-associated lncRNAs. These lncRNAs were significantly associated (P < 0.001) with overall survival (OS). Then, multivariate Cox regression analysis was employed to screen ultimate autophagy-associated lncRNAs and predict the regression coefficients (b) of the model. Finally, a prognosis risk score system based on seven genes was established. Risk score = (b1 × expression level of AC006504.8) + (b2 × expression level of FLVCR1-DT) + (b3 × expression level of AC012306.2) + (b4 × expression level of AC125494.2) + (b5 × expression level of AC005332.6) + (b6 × expression level of AC245041.2) + (b7 ×expression level of LINC02257). Based on an optimal cutoff value, all PC patients were divided into low-and high-risk groups. To estimate the predictive capacities of the risk score system by constructing Kaplan-Meier survival curves and receiver operating characteristic (ROC) curves.

Gene Set Enrichment Analysis (GSEA)
We use GSEA software (http://www.broadinstitute.org/gsea) to identify the underlying different functions between the high-and low-risk groups. The annotated gene sets c5.all.v7.1.symbols.gmt was chosen for the reference gene sets. Enriched gene sets were considered to be statistically significant by a nominal P value < 0.05 that was set as the cut-off criteria.

Cell Culture
The human pancreatic ductal epithelium cell line HPNE and PC cell line PANC-1 were purchased from the ATCC. The cells were cultured at 37°C in a humidified atmosphere containing 5% CO2 with high glucose DMEM supplemented with 10% fetal serum. The culture medium and supplements were purchased from HyClone (Northbrook, IL, USA).

Identification of Seven Prognostic Autophagy-Associated lncRNAs in PC Patients
We extracted a total of 14,142 lncRNAs expression data of tumors from PC tissues in the TCGA database. Two hundred thirty-two autophagy-associated genes were selected from The Human Autophagy Database. We then utilized Pearson correlation analysis to screen the co-expression relationship between the lncRNAs and autophagy-associated genes with the criteria of |R| > 0.4 and P < 0.05. Conclusively, 1,234 autophagy-associated lncRNAs were screened. To identify autophagy-associated lncRNAs related to prognosis, we selected the above filtered autophagy-associated lncRNAs by univariate Cox regression analysis and found that 29 lncRNAs were significantly related to the PC patients' overall survival(OS) ( Table 1). Then, we performed multivariate Cox regression analysis to screen the optimal prognostic lncRNAs. Finally, a total of seven lncRNAs were identified ( Table 2). Among these lncRNAs, AC245041.2 and LINC02257 were risk factors (HR > 1), and AC006504.8, AC012306.2, AC125494.2, FLVCR1-DT, and AC005332.6 were protective factors (HR < 1).

Survival Analysis of Autophagy-Associated lncRNAs
Kaplan-Meier survival analyses and log-rank tests for each autophagy-associated lncRNA were performed to evaluate the prognostic characteristics of patients with PC. K-M survival curves of the seven autophagy-associated lncRNAs are shown that the high expression of AC005332.6, AC006504.8, AC012306.2, AC125494.2, and FLVCR1-DT were positively correlated with the longer overall survival of patients with PC (p < 0.01), indicating protective impacts of these lncRNAs in PC development ( Figure 1A). Furthermore, the high expression of AC245041.2 and LINC02257 was correlated with a short survival time (p < 0.01), which meant that these lncRNAs could play a carcinogenic role in PC ( Figure 1B).

Construction and Validation of the Risk Score Evaluation System of the Autophagy-Associated lncRNAs
Based on seven lncRNAs that were significantly correlated with overall survival, the autophagy-associated lncRNA signature was constructed to predict the outcome of PC patients. The final risk score formula was as follows: Risk score = (-0.3765 × expression level of AC006504.8) + (-0.5525 × expression level of FLVCR1-DT) + (-0.4120 × expression level of AC012306.2) + (-0.5192 × expression level of AC125494.2) + (-0.1040 × expression level of AC005332.6) + (0.2476 × expression level of AC245041.2) + (0.2490 ×expression level of LINC02257). With the above formula, the risk scores of more than 0.1199 were classified into the high-risk group (88 patients); meanwhile, those with less than the cutoff point belonged to the low-risk group (89 patients). The principal components analysis (PCA) showed that the high-risk and low-risk groups were divided into two obvious distribution patterns, which implied that autophagy made distinctly different effects in two groups (Figure 2A). Based on K-M survival analyses, the patients with low-risk scores had longer survival time than those with high-risk scores; OS had statistical significance between the two subgroups ( Figure 2B). We used scatter diagrams to show the risk scores, survival status, and survival time of each PC patient ( Figure 2D). The results demonstrated that the PC patient's survival time and rate gradually deteriorated with increasing risk scores. Moreover, lncRNAs' expression profiles were shown by a heatmap plot ( Figure 2F). We extracted the PACA-CA cohort from the ICGC database to further validated the prognostic stability of the risk evaluation system. Then, we utilized the same formula and cutoff value above to classify the PC patients into low-risk and high-risk groups according to the TCGA dataset. Similarly, the patients from the high-risk group (66 patients) had a lower survival time than the low-risk group (76 patients) ( Figure 2C). Moreover, the scatter diagrams and a heatmap plot of lncRNAs expression profiles were also shown ( Figures 2E, G).
The Autophagy-Related lncRNA Signature Was an Independent Prognostic Factor Subsequently, univariate and multivariate Cox regression analyses were conducted to screen potential biomarkers correlated with OS and total clinical information ( Figures 3A, B). The results showed that only the prognostic value of the risk score was statistically significant. Finally, the receiver operating characteristic (ROC) curve analysis and the AUC value were utilized to assess the prediction accuracy of the above results. The AUC value of risk score based on expression profiles of autophagy-associated lncRNAs was equal to 0.719, which was much higher than age curve (AUC = 0.534), gender curve (AUC= 0.597), grade curve (AUC= 0.607), stage curve (AUC=0.450), T stage curve (AUC=0.504), and N stage curve (AUC= 0.518) ( Figure 3C), suggesting that the risk score was superior to traditional clinical indicators. Thus, the risk score evaluation system derived from the expression levels of the seven lncRNAs was the unique independent prognostic indicator of survival time for PC patients.

Construction of lncRNA-mRNA Network and Enrichment Analysis of GO and KEGG
Based on the abovementioned analysis, to better understand the potential effect of lncRNAs on mRNAs in PC, we built the lncRNA-mRNA network and used Cytoscape and Sankey diagram to visualize the network. We constructed the lncRNA-mRNA coexpression network using the screened seven autophagy-associated lncRNAs with Pearson correlation analysis (|R| > 0.4 and P < 0.05). A total of 61 lncRNA-mRNA pairs were filtrated and the correlation between lncRNAs, mRNAs, and risk score groups by the Sankey diagram ( Figures 4A, B). Furthermore, GO enrichment analysis was performed to clarify the biological processes, cellular components, and molecular function of mRNAs, which were identified from the lncRNA-mRNA co-expression network. As shown in bubble plot revealing top 10 GO terms, We found that the foremost biological processes were "autophagy", "process utilizing autophagic mechanism", and "macroautophagy"; the top three cellular components were "autophagosome", "vacuolar membrane", and "phagophore assembly site"; the top three molecular functions were "ubiquitin protein ligase binding", "ubiquitin-like protein ligase binding", and "protein serine/threonine kinase activity" ( Figure 4C). KRGG enrichment analysis was shown that autophagy, shigellosis, PI3K-Akt signaling pathway, and FoxO signaling pathway were the top four significantly enriched pathways ( Figure 4D).

Gene Set Enrichment Analysis (GSEA)
We carried out the GSEA of the PC samples based on the TCGA to identify the biological pathways associated with the high-risk group and low-risk group. We did not discover a significantly enriched pathway in the high-risk group; moreover, the low-risk group was most significantly enriched for "neuroactive ligand receptor interaction", "tryptophan metabolism", "lysine degradation", "glycosphingolipid biosynthesis ganglio series", "regulation of autophagy", "inositol phosphate metabolism", "glycerophospholipid metabolism", "fatty acid metabolism", etc ( Figure 5). In summary, the GSEA analysis results elaborated that the low-risk score group was closely correlated with autophagy and metabolism. These KEGG data may provide valuable targets to treat for PC.

DISCUSSION
PC is a highly malignant digestive cancer with the lowest fiveyear survival rate in various types of cancer and is predicted to be the second leading cause of cancer-related death in the U.S. by 2030 (10). Most patients with PC cannot be diagnosed early; meanwhile, carbohydrate antigen (CA19-9) as a conventional diagnostic biomarker is not applied to specifically and sensitively diagnose the PC patients (11,12). Only a small part of PC patients can be treated by traditional surgery, and a large number of patients suffer from tumor recurrence and progression. Consequently, the identification of PC regulatory factors has become the focus of recent clinical and basic research. To detect PC early and provide new therapeutic options are of great importance. Currently, LncRNAs have been found to play vital roles in PC and are indispensable for carcinogenetic function, especially autophagy (13,14). An increasing amount of evidence has shown that autophagy-associated lncRNAs make great effects on the occurrence, development, and prognosis of cancer (15).
To better understand the roles of lncRNAs involved with autophagy in the occurrence of PC. Firstly, we analyzed the expression profiles of lncRNAs in PC patients from the TCGA database. We used Pearson analysis to identify the co-expression relationship between lncRNAs and autophagy-related genes in The Human Autophagy Database. Seven autophagy-associated lncRNAs significantly correlated with survival were selected to build the risk score system via the multivariate Cox regression analysis. The high-and low-risk patients can be distinguished according to the median risk score, and PCA analysis displayed a significantly distinct distribution between these two groups. Notably, low-risk patients had a better prognosis than patients in the high-risk group. To assess the potency of the risk score model, we utilized the data from the ICGC database as a validation dataset and got the same result. The AUC value that was calculated from the ROC curves indicated that the risk score had considerable prognostic accuracy for PC patients. Furthermore, the risk score was considered as an independent prognostic factor by univariate and multivariate Cox regression analysis. lncRNAs play an important role in affecting mRNA expression through regulating histone modifications, DNA methylation, and acting as miRNA sponges or precursors of miRNAs, involving the process of transcriptional regulation, post-transcriptional regulation, and epigenetic regulation (16). To elucidate the probable roles of the seven autophagy-associated lncRNAs in PC, the lncRNA-mRNA co-expression network was constructed. The GO and KEGG enrichment analysis was subsequently performed on these mRNAs related to screened lncRNAs, and the results showed that the top enriched GO and KEGG terms were significantly correlated with autophagy. Subsequently, as shown in the GSEA result, we observed that the low-risk group enriched many pathways about lipid, amino acid metabolism, and autophagy, suggesting that metabolism and autophagy were greatly associated with the PC patients of the low-risk score. The above results have shown that the specific autophagy mechanisms were closely related to PC progression.
Besides, there are some limitations that exist in our study. We built the co-expression network between the lncRNAs and mRNA, but how the lncRNAs make specific effects on mRNA was unknown. Furthermore, the autophagy-associated lncRNAs were only detected in PANC-1 and HPNE. The specific molecular mechanisms have not been verified in the experiments. Furthermore, in the existing studies, there are  only a few bioinformatic analyses about some lncRNAs in our risk score system. Chen J et al. found that AC245041.2 was a risk factor in clear cell renal cell carcinoma (ccRCC) and high expression of AC245041.2 was associated with a poor outcome for patients with ccRCC (17). LINC02257 was found to correlate with the prognosis of patients with colorectal cancer and was a risk factor (18). Besides, AC006504.8 was a risk factor in cholangiocarcinoma (19). Concerning the filtered lncRNAs, we need exploratory experiments to prove the functions deeply.
In conclusion, we successfully established the risk score system based on the seven autophagy-associated lncRNAs; meanwhile, it was an independent prognostic factor in PC patients. This approach enhances the prediction accuracy for target lncRNAs, and these autophagy-associated lncRNAs might FIGURE 6 | Expression of 7 lncRNAs in HPNE and PANC-1 cells. The qRT-PCR result showed that AC245041.2 and LINC02257 were high-expressed in PANC-1. AC005332.6, AC012306.2, AC125494.2 were low-expressed in PANC-1. The expression of FLVCR1-DT and AC006504.8 was not statistically different between PANC-1 and HPNE. *P < 0.05, **P < 0.01, NS, no statistically significant. be of great significance for the prediction of prognosis and therapeutic markers for PC patients.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The Human Autophagy Database (http:// www.autophagy.lu/), The Cancer Genome Atlas (https://portal. gdc.cancer.gov/) and International Cancer Genome Consortium (https://icgc.org).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Ethics Committee of Peking Union Medical College Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GC analyzed the data and wrote the manuscript. TZ and YZ conceived the study and obtained financial support. LY and JL applied guiding suggestion, GY, JY, CQ, WL, JQ, and FZ prepared the dataset. All authors contributed to the article and approved the submitted version.