Three Novel Autophagy-Related lncRNAs as Prognostic Biomarkers for Lung Adenocarcinoma

Background: In the present study, autophagy-related long non-coding RNAs (lncRNAs) in lung adenocarcinoma (LUAD) were screened for diagnosis and prognosis, and the molecular mechanisms of LUAD at the genetic level were investigated. Methods: From The Cancer Genome Atlas (TCGA) database, 497 gene expression data and 436 clinical data of LUAD cases were collected for analysis. In addition, 232 autophagy-related genes (ARGs) were extracted from the Human Autophagy Database (HADb). Spearman rank correlation test and the Akaike information criterion (AIC) were performed to screen the data. After filtering, a survival model including three autophagy-related lncRNAs was generated. Based on the following formula: risk score = ΣCoef gene i×Gene i expression, the risk score of all LUAD patients could be calculated. LUAD patients were divided into two groups based on risk score for survival curve using Kaplan-Meier survival analysis. Both univariate and multivariate survival analyses were used to determine whether the three lncRNAs were independent prognostic factors using the survival package in R. Furthermore, the receiver operating characteristic (ROC) curves of clinical data were created to assess the stability of the survival model. Finally, the Gene Set Enrichment Analysis statistically lower than in the low-risk group 41.6% (95% CI: 0.307 - 0.563, P < 0.05). The area under the ROC curve (AUC) of risk score was 0.700, indicating a higher diagnostic accuracy of risk score. The results of GSEA showed enrichment in 36 pathways, including pyrimidine metabolism, pentose phosphate pathway, citric acid cycle, and cell cycle in the high-risk group, and FC-EPSILON-RI signal pathway, intestinal immune network produced by IgA, and ABC transporters in the low-risk group. Conclusion: The prognosis model composed of autophagy-related lncRNAs, AC011477.2, AC099850.3, and TRG-AS1, in LUAD can be used to predict the prognosis of LUAD patients and is expected to improve clinical treatment.


Introduction
Lung cancer is currently the most common and lethal malignancy worldwide.
The 2-year relative survival rate for non-small cell lung cancer (NSCLC) increased from 34% diagnostic rate between 2009 and 2010 to 42% between 2015 and 2016, including an absolute increase of 5% to 6% for each diagnostic stage [1]. Globocan 2020 data indicate the incidence and mortality of lung cancer are 37.0% worldwide and 39.8% in China [2]. Different histological subtypes have various origins, morphological appearance, epigenetic changes, clinical responses, and outcomes [3]. Therefore, investigating diagnostic and prognostic genes or markers for different subtypes could help implement precision therapies and improve patient survival.
Autophagy and apoptosis play important roles in lung cancer progression [4].
Eukaryotic cells respond to different stimuli and stress, such as starvation, hypoxia, drugs, and other harmful stimuli. To maintain cellular homeostasis, the degradation process of eukaryotic cells is regulated by autophagy-related genes (ARGs) [5]. Autophagy is a lysosomal degradation pathway physiologically involved in differentiation, survival, development and homeostasis, and pathologically in infection, cancer, neurodegeneration and aging as well as heart, liver, and kidney disease [6]. Notably, lncRNAs may act as competing endogenous RNAs (ceRNAs) which competitively bind with target miRNAs to regulate mRNAs [7]. In several studies, autophagy was shown to promote the tolerance of tumor cells to radiotherapy and chemotherapy, and lncRNAs play a crucial role in these processes [8] [9].
In the present study, the autophagy-related lncRNAs were mined from The Cancer Genome Atlas (TCGA) by searching for the lncRNA co-expressed genes of ARGs in LUAD and identified the autophagy lncRNAs associated with patient clinical data to investigate the molecular mechanisms of LUAD.

The Gene Expression and Clinical Data of LUAD Patients
The gene expression and clinical data of LUAD patients were downloaded from TCGA database (https://portal.gdc.cancer.gov/). In this study, among 551 LUADs, 497 tumor samples were compared with 54 adjacent tissues, and their gene expression and clinical data were collected. In a previous study, the M stage in the TNM staging system was suggested to be an unreliable independent prognostic factor [10], therefore, the missing and the M stage clinical data were excluded. In addition, after filtering, 497 cases of gene expression and 436 cases of clinical data were analyzed.

The Autophagy-Related lncRNA Prognostic Model
The autophagy-related lncRNAs were combined with survival data and then a univariate Cox's proportional hazard analysis was performed as well as Kaplan-Meier analysis using the survival package in R. A P-value < 0.05 was considered statistically significant. To determine the optimum prognostic autophagy-related lncRNAs, the minimum Akaike information criterion (AIC) value was screened. Next, the OS curves of the optimum autophagy-related lncRNAs were drawn based on the gene expression levels. The LUAD patients were divided into high-risk and low-risk groups based on the median risk score calculated using the following formula: risk score = ΣCoef gene i × Gene i expression. Survival curves for each group were constructed based on Kaplan-Meier analysis.

Stability Evaluation of the Prognostic Model
To assess the stability of the constructed prognostic model, the clinical data were combined with the risk score of each LUAD patient. The clinical data included age, gender, clinical stage, T stage, and N stage. Then, univariate and multivariate prognostic analyses were performed to determine whether the three lncRNAs were independent prognostic factors using the survival package in R.
The area under the ROC curve (AUC) of risk score and different clinical data were created to assess the stability for 1-year survival time using the survival ROC package in R.

Gene Set Enrichment Analysis (GSEA)
The Gene Set Enrichment Analysis (GSEA) was used to analyze the related pathways by choosing "c2.cp.kegg.v7.0.symbols.gmt" gene sets as reference gene sets. The above analysis was performed using GSEA 4.1.0 software.

The Optimum Three Prognostic Autophagy-Related lncRNAs
The raw data contained 19

Construction of Prognostic Model
The risk score of all LUAD patients was calculated based on the following formula: risk score = ΣCoef gene i × Gene i expression.
Based on the median risk score (cutoff value = 1.02407), the LUAD patients were divided into high-risk and low-risk groups. A Sankey diagram of autophagy-related mRNA-lncRNA co-expression network was established (Table 3 and    To visualize the data, the risk curve (Figure 3(a)), survival status scatter plot ( Figure 3(b)), and the heat map of the three lncRNAs' expression ( Figure 3(c)) were drawn.

Stability Validation of the Prognostic Model Composed of the Three Autophagy-Related lncRNAs
To evaluate the stability of the prognostic model composed of the three lncRNAs,   score was 0.700, indicating a higher diagnostic accuracy of risk score. Furthermore, the AUC of the clinical staging was 0.715, indicating that applying clinical staging for diagnosis was equally highly valuable.
In conclusion, the stability of this prognostic model was acceptable and could be applied for the prognosis of LUAD patients.
In addition, clinical data were divided into different groups based on correlation with risk score to perform differential analysis using t-test (Table 4).
Women had lower risk than men, the risk of clinical stage I -II was less than of stage III -IV, the risk of T1 -T2 was less than of T3 -T4, and the risk of N0 was less than of N1 -N3 (P < 0.05).

GSEA
GSEA results showed different pathways were enriched in the high-risk and low-risk groups ( Figure 6). KEGG reference gene sets were performed and results indicated enrichment in 36 pathways including pyrimidine metabolism, pentose phosphate pathway, citric acid cycle, and cell cycle in the high-risk group (P < 0.05) and enrichment in FC-EPSILON-RI signal pathway, intestinal immune network produced by IgA, and ABC transporters in the low-risk group (P < 0.05).

Discussion
Over the past few decades, researchers could investigate genomes more systematically due to application of high-throughput technologies [11]. Bioinformatics, interdisciplinary medicine, statistics, computer science, and applied mathematics, are disciplines that contributed to the development of the human genome project [12].
In several studies, autophagy has been shown to affect gene expression and even clinical treatment. Autophagy was shown to affect hepatocarcinogenesis and progression in several studies by regulating liver cancer immunity, oxidative stress, and cellular metabolism, including suppressing hepatocarcinogenesis at an early stage, promoting tumor growth at a progressive stage, and increasing chemoradiotherapy resistance of cancer cells [13] [14] [15]. Zhang et al. [16] found that lncRNA CASC2 might be useful for evaluating the expression of autophagy-related LC3-I, LC3-II, and p62. Several lncRNA-related models have been applied for the prognosis of patients with malignant tumors, such as lung, colon, liver, and breast cancers [17] [18] [19]. However, autophagy-related lncRNAs have not been extensively investigated.  In the present study, a stable prognostic model of three autophagy-related lncRNAs (AC011477.2, AC099850.3, and TRG-AS1) in LUAD was established.
AC099850.3 has been included in some hepatocellular carcinoma (HCC) prognosis studies, and the results showed that AC099850.3 may play an important role in the migration and proliferation of HCC cells. Wu et al. [20] found that AC099850.3 could promote the expression of cell cycle molecules BUB1, CDK1, PLK1, and TTK based on RT-PCR, and could inhibit CD155 and PD-L1 expression in the interference group based on Western blot analysis. The results indicated that AC099850.3 played a specific role in the immune response and warranted further research. Zhou et al. [21] showed that AC099850.3 was upregulated (logFC > 1) in tongue squamous cell carcinoma (TSCC).
He et al. [22] found that TRG-AS1 was highly expressed in TSCC patients and strongly associated with advanced TNM staging and poor OS. Sun et al. [23] demonstrated that TRG-AS1 was significantly overexpressed in HCC. Reportedly, TRG-AS1 silencing might inhibit the proliferation, migration, invasion, and epithelial mesenchymal transition (EMT). Xie et al. [24] found that TRG-AS1 regulated SUZ12 expression through the ceRNA of miR-877-5p, and TRG-AS1 might become a novel target for glioblastoma treatment. Notably, in the present study, TRG-AS1, a low-risk lncRNA for LUAD, had the opposite results in TSCC, liver cancer, and glioblastoma, which requires further research.
Currently, related studies in which lncRNA AC011477.2 has been investigated have not been performed. However, the present study can be used as a basis for future research on lncRNA AC011477.2 and clarify whether it has the potential to become a molecular diagnostic marker and tumor therapeutic target for LUAD.
GSEA enrichment analysis showed the high-risk group was enriched in the cell cycle, which could be associated with some of the above-mentioned effects of AC099850.3. However, further experimental proof is required. The results of the present study can be used for future prospective experiments to verify whether this prognostic model can be applicable for evaluating the prognosis of patients and identifying the specific mechanism(s) of these three lncRNAs in LUAD patients.

Conclusion
In the present study, an autophagy-related lncRNA prognostic model for LUAD consisting of AC011477.2, AC099850.3, and TRG-AS1 was established. The model showed acceptable stability for predicting the prognosis of LUAD patients and can be used for subsequent prospective experiments to improve clinical treatment.

Funding Information
Funding was provided by The National Natural Science Foundation of China