Identification of novel autophagy-related lncRNAs associated with a poor prognosis of colon adenocarcinoma through bioinformatics analysis

LncRNAs play a pivotal role in tumorigenesis and development. However, the potential involvement of lncRNAs in colon adenocarcinoma (COAD) needs to be further explored. All the data used in this study were obtained from The Cancer Genome Atlas database, and all analyses were conducted using R software. Basing on the seven prognosis-related lncRNAs finally selected, we developed a prognosis-predicting model with powerful effectiveness (training cohort, 1 year: AUC = 0.70, 95% Cl = 0.57–0.78; 3 years: AUC = 0.71, 95% Cl = 0.6–0.8; 5 years: AUC = 0.76, 95% Cl = 0.66–0.87; validation cohort, 1 year: AUC = 0.70, 95% Cl = 0.58–0.8; 3 years: AUC = 0.73, 95% Cl = 0.63–0.82; 5 years: AUC = 0.68, 95% Cl = 0.5–0.85). The VEGF and Notch pathway were analyzed through GSEA analysis, and low immune and stromal scores were found in high-risk patients (immune score, cor =  − 0.15, P < 0.001; stromal score, cor =  − 0.18, P < 0.001) , which may partially explain the poor prognosis of patients in the high-risk group. We screened lncRNAs that are significantly associated with the survival of patients with COAD and possibly participate in autophagy regulation. This study may provide direction for future research.

identified, the co-expression network linking these autophagy-related lncRNAs and autophagy-related genes was established. Cytoscape (version 3.4; The Cytoscape Consortium (San Diego, CA, USA) was then used in visualizing the PPI networks 14 . The "clusterProfiler" package in R software was used in GO and KEGG enrichment analyses (www. kegg. jp/ kegg/ kegg1. html) 15 . A P value of < 0.05 was regarded as statistically significant, and the top ten results of enrichment analyses were selected for visualization.

Nomogram and calibration.
We included the clinical features of age, gender, stage, and TNM classification in our analysis and considered the practical utility of our model in clinical processes. A nomogram plot was constructed using the "rms" package in R software for the prediction of 1-, 3-, and 5-year survival rates of patients with COAD. Calibration plots were used in assessing the prognostic accuracy of the nomogram. Specifically, nomogram-predicted and actual probabilites of the patients were compared.
Gene set enrichment analysis and gene set variation analysis. After the samples were divided into high-and low-risk score groups, GSEA was conducted for the link the genes to feasible pathways 16 . Gene set permutations were performed 1000 times for each analysis. The enriched pathways were selected under the following screening conditions: FDR < 0.25 and NOM P value < 0.05. Gene set variation analysis (GSVA) was performed using the GSVA package R software with the Hallmark dataset.

Clinical correlation and tumor microenvironment analysis. Wilcox test was used in investigat-
ing the clinical relevance of the seven lncRNAs, risk scores, and clinical features. The "Estimate" package and ssGSEA algorithm were used in the calculation of stromal and immune scores in tumor microenvironment.

Result
Screening for autophagy-genes and -lncRNAs. A total of 231 autophagy-related genes were identified from the website of HADb (Table S1 and Fig. 1A), then 1285 autophagy-related lncRNAs were screened at a set threshold of correlation coefficient |R 2 | of > 0.3 and P value of < 0.001 (Fig. 1B).

Co-expression network and enrichment analysis.
A total of 47 autophagy-related genes were found to be associated with the seven prognosis-related lncRNAs. Furthermore, the co-expression network based on the relationships of these genes with the lncRNAs were established with 54 nodes and 51 edges ( Fig. 4A and Figure S2). The top 20 significant nodes were selected according to the MCC value in the cytohubba analysis, and the lncRNA AL161729.4 was the most important node (Fig. 4B). The Sankey diagram intuitively showed the association of each node with patients' prognosis ( Fig. 4C). GO and KEGG analyses were conducted for the functional enrichment of the nodes in this co-expression network (Fig. 4D-E). The result revealed that for biological processes, the nodes were mainly enriched in "autophagy, " "macro autophagy, " "response to nutrient lev- www.nature.com/scientificreports/ els, " and "regulation of autophagy. " Changes in cellular components were markedly enriched in "vacuolar membrane, " "phagocytic vesicle, " "endocytic vesicle, " and "membrane raft. " Changes in the DEG molecular function were primarily enriched in "protein serine/threonine kinase activity, " "protein serine/threonine/tyrosine kinase activity, " "chaperone binding, " and "heat shock protein binding. " KEGG analysis results showed that DEGs were strikingly enriched in "Autophagy-animal, " "Alzheimer's disease, " "Shigellosis, " and "Kaposi sarcoma-associated herpes virus infection. "  GSEA and GSVA analysis. The biological functions of the high-and low-risk groups were further explored through GSEA analysis. As is shown in Fig. 6, the VEGF and Notch signaling pathways were enriched in the high-risk phenotype. In the low-expression phenotype, the top five enriched gene sets were glycan biosynthesis, alanine aspartate and glutamate metabolism, pentose and glucuronate interconversions, ascorbate and aldarate metabolism, and metabolism of seven amino acids. The underlying signaling pathway of seven lncRNAs was analyzed through GSVA analysis ( Figure S1).  The mRNA expression of seven lncRNAs in COAD cell lines. We evaluated the mRNA level of seven lncRNAs (AC027307.2, AC073611.1, AC156455.1, AL161729.4, LINC01063, MIR210HG and PCAT6) in five COAD cell lines and normal colon mucosal epithelial cell line (Fig. 8). The result revealed that AC027307.2 has a lower expression level in HCT116 and DLD-1 cell lines than CCD-18Co (Fig. 8A). The mRNA level of AC073611.1, MIR210HG, AL161729.4 and LINC01063 has no statistically significant difference between cancer and normal cell lines (Fig. 8B-E). AC156455.1 mRNA level is significantly down-regulated in SW480, LS174T and HT29 cell lines (Fig. 8F). A higher mRNA level of PCAT6 is observed in SW480, LS174T and HCT116 cell lines (Fig. 8G).

Discussion
COAD is one of the most common cancers worldwide and responsible for more than 600,000 deaths each year 17 . Many factors, such as high-fat and low-fiber diets and genetics, are now widely recognized risk factors for COAD [18][19][20] . Despite the progress in the development of effective diagnostic and therapeutic strategies for COAD, the lack of detection in early status and subsequent invasiveness have rendered the disease a persistent problem for susceptible populations 21 . Therefore, useful biomarkers for the early diagnosis and prognosis prediction of COAD are crucial. Although the role of autophagy's in cancer treatment remain unclear, substantial evidence suggests the immense potential of autophagy therapy as a novel approach for COAD treatment 22 . Meanwhile, with a wide range of functional activities, lncRNAs may play a pivotal role in physiological processes, such as RNA decay, genetic regulation of gene expression, RNA splicing, and protein folding 23,24 . They can regulate many proteins that are essential to autophagy. In this study, basing on the open-access data obtained from the TCGA database (TCGA-COAD), we systematically studied the association between autophagy-related lncRNA and COAD through bioinformatics analysis. We aimed to screen signatures that are useful in predicting the development of COAD and guiding therapy strategy becaue these signature might be novel prognostic markers. To the best of our knowledge, this is the first study that focuses on the role of autophagy-related lncRNAs in COAD treatment and tumor microenvironment regulation.
First, we identified seven prognostic-related lncRNAs by respectively conducting univariate cox analysis, LASSO regression, and multivariate cox analysis. We divided the patients with COAD into high-and low-risk  www.nature.com/scientificreports/ groups according to their median risk scores. The low-risk group had a longer OS. Through univariate and multivariate COX regression analysis, we were able to conclude that our prognostic model is efficient and independent of other clinical factors, such as TNM classification, clinical stage, age, and gender. Next, we constructed a co-expression network of autophagy-related lncRNAs and mRNAs. Most of the nodes of the network were reported for the first time. Currently, autophagy is still an emerging field in preliminary basic studies, and the implication of lncRNAs in autophagy is extremely important 25 . For example, Liu et al. revealed that under energy stress, the lncRNA NBR2 can interact with AMPK and promotes AMPK kinase activity, subsequently activating autophagy in cancer cells 26 . Moreover, the lncRNAs HOTAIRM1, PTENP1, and MALAT1 were found to be involved in the activation of autophagy and regulation of several physiological processes and malignant phenotype of cancer cells [27][28][29][30] . By contrast, the down-regulated lncRNA Risa can improve insulin sensitivity by enhancing autophagy 31 . In COAD, on the one hand, autophagy is involved in tumor development and drug resistance 32 . On the other hand, as non-canonical regulators, lncRNAs play a pivotal role in the physiological  Our results showed that patients in the M1 or N1-2 stage a likely to have high risk scores. This condition partially explains the observed poor prognosis impact. Furthermore, we performed GSEA analysis to explore the difference in biological function and pathways between the high-and low-risk groups. Notably, in the high-risk phenotype, the VEGF and Notch signaling pathways were enriched. Many researchers have observed that tumor development pro-angiogenic factors, such as VEGF, basic and acidic fibroblast growth factor, tumor necrosis factor-α, and interleukin-1 are enhanced in multiple pathological processes. The persistent growth of tumordirected capillary networks creates a favorable microenvironment, promoting cancer growth, progression, and metastasis 38 . Notch signaling is a significant regulator of sprouting angiogenesis, which is controlled by a tightly regulated balance between endothelial tip and stalk cells 39 . Moreover, the Notch site activated by adjacent vascular cells in cancer cells increases migration across endothelial cells, thus increasing metastasis in CRC patients 40 . Finally, we evaluated the effect of risk scores on the tumor microenvironment (stromal score and immune score), which have essential roles in COAD development 42 . Differences among the components of immune and stromal cells can substantially affect patients' survival [42][43][44] . In our study, a negative relation was found between risk and tumor environment scores (stromal and immune). Meanwhile, the patients in the TGCA cohort with low immune and stromal scores may have poor prognoses. These results provide a new perspective for the poor prognosis of high-risk COAD patients.
Some limitation still exists in our study. First, the data series obtained for analysis are primarily from Western countries, and thus the results of the study may not fully apply to patients in Asian countries, given the difference in genetics between races. Second, the amount of data published in the public database is limited, and thus the clinical pathology parameters used for analysis in this study are not comprehensive and may lead to potential errors or biases.

Conclusions
Through serial bioinformatics analysis, we identified autophagy-related lncRNAs that markedly affect the prognoses of patients. Basing on these lncRNAs, we established an effective prognosis model for predicting the OS of patients with COAD. Moreover, the co-expression network-linked autophagy-related lncRNAs and genes can provide direction for future research.

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. The expression profile of 1285 autophagy-related lncRNAs has been uploaded on the website https:// figsh are. com/ artic les/ datas et/ lncRN As/ 13490 610 (figshare).