An autophagy-related long non-coding RNA prognostic signature accurately predicts survival outcomes in bladder urothelial carcinoma patients

In this study, we analyzed the prediction accuracy of an autophagy-related long non-coding RNA (lncRNA) prognostic signature using bladder urothelial carcinoma (BLCA) patient data from The Cancer Genome Atlas (TCGA) database. Univariate and multivariate Cox regression analyses showed significant correlations between five autophagy-related lncRNAs, LINC02178, AC108449.2, Z83843.1, FAM13A-AS1 and USP30−AS1, and overall survival (OS) among BCLA patients. The risk scores based on the autophagy-related lncRNA prognostic signature accurately distinguished high- and low-risk BCLA patients that were stratified according to age; gender; grade; and AJCC, T, and N stages. The autophagy-related lncRNA signature was an independent prognostic predictor with an AUC value of 0.710. The clinical nomogram with the autophagy-related lncRNA prognostic signature showed a high concordance index of 0.73 and accurately predicted 1-, 3-, and 5-year survival times among BCLA patients in the high- and low-risk groups. The lncRNA-mRNA co-expression network contained 77 lncRNA-mRNA links among 5 lncRNAs and 49 related mRNAs. Gene set enrichment analysis showed that cancer- and autophagy-related pathways were significantly enriched in the high-risk group, and immunoregulatory pathways were enriched in the low-risk group. These findings demonstrate that an autophagy-related lncRNA signature accurately predicts the prognosis of BCLA patients.


INTRODUCTION
AGING [7], cardiovascular diseases [8], and inflammatory disorders related to infectious diseases [9]. Autophagy is associated with tumor suppression or oncogenesis depending upon the stage of tumor development [10,11]. Recent studies show that modulation of autophagy improves the sensitivity of BCLA tumors to chemotherapeutic agents [12,13]. Hence, it is critical to discover autophagy-related biomarkers that can serve as valuable the early diagnostic and prognostic biomarkers for BCLA patients.
Genome sequencing studies show that nearly 90% of the human transcriptome represents non-coding RNA or ncRNA [14]. The long non-coding RNAs (lncRNAs) are a type of ncRNAs with transcripts of >200 nucleotides in length without any protein-coding capacity [15]. LncRNAs regulate important biological functions related to cell growth and survival, genomic imprinting, chromatin modifications, and allosteric regulation of enzyme activities [16]. Furthermore, pathogenesis of several human diseases including several different types of cancers involves dysregulation of specific lncRNAs [17]. Some studies have shown that lncRNAs regulate autophagic functions. For example, Ying et al. demonstrated that downregulation of lncRNA MEG3 promotes proliferation of bladder cancer cells by activating autophagy [18]. Another study shows that lncRNA MALAT1 regulates multi-drug resistance of hepatocellular carcinoma cells by altering autophagy [19]. LncRNA PVT1 promotes in vitro and in vivo pancreatic ductal adenocarcinoma progression by activating autophagy through its regulation of the miR-20a-5p/ULK1 axis [20].
New advances in genome sequencing technology and bioinformatics have helped to identify potential prognostic biomarkers that can predict survival outcomes in cancer patients [21,22]. Therefore, we postulated that autophagy-related lncRNAs may be valuable prognostic biomarkers for BLCA patients. In this study, we systematically analyzed the relationship between the expression of autophagy-related lncRNAs and the clinicopathological characteristics of 409 BLCA patients from The Cancer Genome Atlas (TCGA) database. We also constructed a prognostic signature based on 5 autophagy-related lncRNAs and evaluated its ability to independently and accurately predict the prognosis of BLCA patients.

Identification of prognostically significant autophagyrelated lncRNAs in BLCA patient tissue samples
We identified 14153 lncRNAs by analyzing the RNA-seq data of the BLCA patient tissue samples from the TCGA database. We also extracted 232 autophagy-related genes from the Human autophagy database (HADb) analysis. We then identified 49 autophagy-related lncRNAs by performing Pearson correlation analysis between the lncRNAs and the autophagy-related genes using |R| > 0.7 and P < 0.05 as the selection criteria. Univariate Cox regression analysis of the 49 autophagy-related lncRNAs showed that expression of 7 lncRNAs, namely, AC002553.2, Z83843.1, LINC02178, FAM13A−AS1, USP30−AS1, AC108449.2 and AC243960.1 significantly correlated with the overall survival (OS) of BLCA patients (P < 0.05; Figure 1A). Multivariate Cox regression analysis showed that 5 of the 7 autophagyrelated lncRNAs were good candidates for constructing the prognostic signature based on the lowest Akaike information criterion (AIC) ( Table 1). Among the 5 autophagy-related lncRNAs that were included in the prognostic signature, LINC02178 and AC108449.2 were considered as risk factors with HR values greater than 1, whereas the remaining 3 lncRNAs, Z83843.1, FAM13A−AS1 and USP30−AS1, were considered as protective factors with HR values less than 1.

Evaluation of the prognostic signature containing 5 autophagy-related lncRNAs
The risk score for each BCLA patient in the TCGA dataset was calculated using the following formula for the autophagy-related lncRNA signature: risk score = (-0.677 × expression level of Z83843.1) + (0.162 × expression level of LINC02178) + (-0.403 × expression level of FAM13A−AS1) + (-0.307 × expression level of USP30−AS1) + (0.489 × expression level of AC108449.2). Then, BLCA patients were divided into high-risk (n = 196) and low-risk (n = 197) groups using the median risk score (= 1.093) as the cut-off point. Kaplan-Meier survival curve analysis showed that the OS of BCLA patients with high-risk scores was significantly shorter than those with low-risk scores ( Figure 1B). The 3-year survival rates were 39% and 64%, and the 5-year survival rates were 32% and 56% for the high-risk and low-risk patients, respectively. A principal components analysis (PCA) based on the five autophagy-related lncRNAs showed two significantly different distribution patterns between high-risk and low-risk groups ( Figure 1C). Time-dependent receiver operating characteristic (ROC) curve analysis showed that the area under the ROC (AUC) value for the autophagy-related lncRNA prognosis signature was 0.710 ( Figure 1D). BCLA patients were then ranked according to the risk scores calculated based on the autophagy-related lncRNA prognosis signature ( Figure  1E). The scatter dot plot showed that the survival rates of the BCLA patients correlated with the risk score according to the autophagy-related lncRNA prognostic signature; patients with a higher risk score demonstrated AGING lower survival time ( Figure 1F). The heatmap showed distinct differences in the levels of the 5 prognostic signature-related lncRNAs in the high-and low-risk BCLA patients. High-risk patients expressed higher levels of risk factors (AC108449.2 and LINC02178), while low-risk patients expressed higher levels of protective factors (Z83843.1, FAM13A−AS1 and USP30−AS1) ( Figure 1G).

Correlation analysis of the autophagy-related lncRNA prognosis signature with other clinicopathological parameters
We then analyzed the correlation between the risk scores from the autophagy-related lncRNA prognosis signature and the clinicopathological characteristics of the BCLA patients from TCGA database. Patients aged > 65 years showed significantly higher risk scores compared to patients aged ≤ 65 years ( Figure 2A). The risk scores were statistically similar between the male and female BCLA patients ( Figure 2B). Furthermore, the risk scores were statistically similar for BCLA patients belonging to high-and low-grades, probably because majority of the patients analyzed belonged to high-grade group (high-grade, n = 372; low-grade, n = 18; Figure 2C). Moreover, BCLA patients belonging to the higher AJCC stages showed higher risk scores than those with lower AJCC stages ( Figure 2D). These results demonstrate that the autophagy-related lncRNA risk signature is associated with the clinicopathological characteristics of BCLA patients.

AGING
probably because of the smaller sample size. These results suggest that the prognosis signature can accurately determine the prognosis of patients relative to other clinicopathological characteristics.

The autophagy-related lncRNA signature is an independent prognostic factor
Next, we performed univariate and multivariate Cox regression analyses to determine if the autophagy- AGING related lncRNA prognostic signature was an independent prognostic factor for patients with BLCA. Univariate analyses showed that age (P < 0.001), AJCC stage (P < 0.001), T stage (P < 0.001), N stage (P < 0.001) and autophagy-related lncRNA prognostic risk score (P < 0.001) were significantly associated with OS ( Figure 4A). The HR value tended to infinity within the tumor grade because of uneven distribution of samples (18 cases in low grade and 372 cases in high grade). Multivariate analyses showed that age (P < 0.001) and autophagy-related lncRNA prognostic risk score (P < 0.001) were significantly associated with OS ( Figure 4B). As shown in Figure  4C, the ROC curve analysis demonstrated that the AUC value for the autophagy-related lncRNAs prognostic signature was 0.710, which was higher than the AUC values for age (AUC = 0.627), gender (AUC= 0.526), grade (AUC= 0.537), AJCC stage (AUC=0.688), T stage (AUC=0.605) and N stage (AUC= 0.651). These data demonstrate that the autophagy-related lncRNA prognostic signature is an independent prognostic factor for BLCA patients.

Evaluation of the prognostic prediction nomogram that includes autophagy-related lncRNA prognostic signature risk score
Nomograms are commonly used tools used by clinicians to accurately predict survival time for a patient by calculating the nomogram score based on the points assigned for each prognostic factor included in the nomogram [23]. We constructed a nomogram to accurately estimate the 1-, 3-, and 5-year survival probabilities by using risk score calculated from the autophagy-related lncRNA prognostic signature and other clinicopathological factors, including age, gender, grade, AJCC stage, T stage and N stage ( Figure 5A). The concordance index (C-index) value for the nomogram was 0.715. The calibration curve analysis showed that the actual and the predicted 1-, 3-, and AGING 5-year survival times were in agreement when compared with the reference line ( Figure 5B-5D). These results demonstrated that the nomogram using the autophagy-related lncRNA prognostic signature risk scores was reliable and accurate.

Construction of the lncRNA-mRNA co-expression network and functional enrichment analysis
Next, we investigated the potential functions of the 5 autophagy-related lncRNAs in BLCA by constructing the lncRNA-mRNA co-expression network using Cytoscape. The lncRNA-mRNA co-expression network contained 77 lncRNA-mRNA pairs based on the threshold parameters, Pearson correlation coefficient |R| > 0.3 and P < 0.05 ( Figure 6A). Among these, 49 mRNAs significantly correlated with the 5 lncRNAs in the prognostic signature. The Sankey diagram showed the relationship between the 49 mRNAs and 5 lncRNAs (risk/protective) ( Figure 6B). The top three GO terms for the biological processes were autophagy, process utilizing autophagic mechanism, and macroautophagy ( Figure 6C). The top three GO terms for the cellular components were cytosolic part, PML body, and the nuclear envelope ( Figure 6D). The top three GO terms for molecular functions were protein serine/threonine kinase activity, ubiquitin protein ligase binding, and ubiquitin−like protein ligase binding ( Figure 6E). KEGG pathway analysis confirmed that autophagy was the most significantly enriched pathway ( Figure 6F).

Gene set enrichment analysis
GSEA results showed that the altered genes in the highrisk BCLA patients belonged to pathways related to autophagy and cancer, WNT signaling pathway, renal cell carcinoma, TGF-βsignaling pathway, VEGF signaling pathway, ERBB signaling pathway, PPAR signaling pathway, MAPK signaling pathway, P53 signaling pathway, mTOR signaling pathway, endocytosis, RNA degradation and ubiquitin-mediated proteolysis ( Figure 7A). Immunoregulatory pathways against cancer were significantly enriched in the lowrisk group, including pathways related to antigen processing and presentation, natural killer (NK) cellmediated cytotoxicity, T cell receptor (TCR) signaling, AGING chemokine signaling and B cell receptor (BCR) signaling ( Figure 7B). This suggested that activation of pathways regulation immune function in the low-risk group may contribute to positive prognosis or longer survival outcomes. The top 10 KEGG pathways in the high-risk and low-risk groups based on GSEA are shown in Figure 7C and 7D. These results suggested that a high prognostic signature risk score correlates with autophagy and cancer, whereas low prognostic signature risk score correlates with enhanced immune function. These data provided valuable insights for future investigations into potential individualized treatments for BLCA patients belonging to different risk groups.

DISCUSSION
The most common malignancy of the urinary system is BLCA, whose incidence rates are constantly increasing worldwide [1,3]. The prognosis of BCLA patients is poor because of late diagnosis and high rate of therapeutic resistance [24]. The role of autophagy in tumorigenesis has been reported for several cancers, including BLCA [25]. Therefore, autophagy-related biomarkers are potential diagnostic biomarkers and therapeutic targets for BCLA patients. Previous studies have focused on the role of specific autophagy-related genes in BCLA progression [26].
LncRNAs are a new class of non-coding RNA molecules that regulate cancer cell growth, progression, and survival [27]. Hence, they are potential biomarkers that can predict cancer risk and survival outcomes. In this study, we systematically analyzed the prognostic prediction accuracy of autophagyrelated lncRNAs in BCLA using bioinformatics and statistical tools.
We first identified 7 autophagy-related lncRNAs that significantly correlated with OS based on the univariate Cox regression analysis of the expression of autophagyrelated lncRNAs in the BCLA patient samples from the TCGA database. Furthermore, 5 autophagy-related lncRNAs, Z83843.1, LINC02178, FAM13A−AS1, USP30−AS1 and AC108449.2 were selected to construct a prognostic signature based on their performance in the multivariate Cox regression analysis. The risk score of each BCLA patient was calculated according to the expression of the five autophagy-related lncRNAs in the prognostic signature and patients were classified into high-and low-risk AGING groups based on the median risk score. BCLA patients with high-risk scores showed shorter survival times compared to those with low-risk scores. Principal components analysis (PCA) based on the confirmed five autophagy-related lncRNAs showed two significantly different distribution patterns between high-risk and lowrisk groups. ROC curve analysis validated the prognostic accuracy of the autophagy-related lncRNA prognostic signature in the BLCA patients. The risk score based on the autophagy-related lncRNA prognostic signature was an independent prognostic factor based on multivariate Cox regression analysis. Stratified correlation analysis showed that the autophagy-related lncRNA prognostic signature accurately predicted survival outcomes for the high-and low-risk BCLA patients.
The autophagy-related lncRNA prognostic signature performed more reliably than the other traditional clinical AGING indicators in prognostic prediction. A nomogram is an effective and reliable clinical tool to predict survival of cancer patients [28]. Therefore, we developed a robust nomogram consisting of several clinical variables (age, gender, grade, AJCC stage, T stage and N stage) and the risk scores based on the autophagy-related lncRNA prognostic signature to improve prognostic prediction of BCLA patients. Older patients (age ≥ 60 years) and those with higher tumor grades and advanced stages are usually associated with worse cancer prognosis, which is consistent with our results. Moreover, calibration plots demonstrated that the actual and predicted 1-, 3-, and 5-year survival rates based on the nomogram were similar. Overall, the autophagyrelated lncRNA prognostic signature accurately predicts survival outcomes of BLCA patients in our study and shows great potential for clinical applications, including individualized prognosis and therapy.
Autophagy is a highly conserved intracellular catabolic process involved in the phagocytosis and degradation of abnormal organelles, proteins and pathogens through the lysosomal pathway [29]. The role of autophagy in cancer is controversial because it can play both tumor suppressor and oncogenic functions [30]. During early stages of tumor development, autophagy-related cell death can suppress tumor progression; autophagic dysregulation can also induce genomic instability and necrosis-induced inflammation, both of which promote tumor growth [31]. Conversely, autophagy sustains tumor metabolism, growth, and survival in nutrientdeprived conditions in the tumor microenvironment and contributes to drug resistance during tumor metastasis [32]. Autophagy is regulated by several signaling pathways, including the PI3K/AKT/mTOR signaling pathway [33] and ubiquitin-proteasome system (UPS). [34] In recent years, several lncRNAs have been implicated in the regulation of cell growth and survival by directly targeting autophagy-related genes. For example, Wang et al. reported that lncRNA ATB induced autophagy by enhancing the expression of autophagy-related protein 5 (ATG5) through activation of the Yes-associated protein (YAP) in hepatocellular carcinoma (HCC) cells [35]. We identified the genes whose expression is regulated by each of the 5 autophagy-associated lncRNAs in BLCA and constructed the lncRNA-mRNA co-expression network. GO, KEGG, and GSEA functional enrichment analyses showed that autophagy-related GO terms or signaling pathways were enriched. GSEA analyses also revealed distinct differences in the autophagy-related signaling pathways between the high-and low-risk groups. Several cancer-and autophagy-related pathways were enriched in the high-risk group, whereas immunomodulatory pathways were enriched in the low-risk group. This suggested that increased immunity correlates with improved prognosis. These results were concordant with the current understanding that autophagy is a critical modulator of BLCA progression [12,13].
There are several limitations in our study. Firstly, our findings need to be further validated in other independent cohorts to determine the robustness of the autophagy-related lncRNA prognostic signature. Secondly, our study was based on a single cohort of 409 patients from the publicly available TCGA database. Moreover, samples belonging to BCLA patients with high-grade tumor (n = 385) were significantly larger than those with low-grade tumors (n = 21), which may have skewed our results and hence need to be further analyzed with larger and more even number of samples in the high-risk and low-risk groups. Finally, further investigations involving biochemical experiments such as immunohistochemistry, quantitative real-time PCR, and flow cytometry, and clinical data analyses are required to further confirm our findings.
In conclusion, our study showed that the autophagyrelated lncRNA prognostic signature accurately predicts the survival outcomes of BCLA patients with BLCA and distinguishes them into high-and low-risk groups. We also established and validated a prognostic nomogram by combining the autophagy-related lncRNA prognostic signature and other clinicopathological features. Our study demonstrated that this nomogram can provide an individualized and accurate survival prediction. Our study also suggests that these 5 autophagy-related lncRNAs are potential prognostic and diagnostic biomarkers as well as promising targets for BCLA therapy.

Patient data acquisition
The raw RNA-sequencing (RNA-seq) data and clinical information of 409 BLCA patients was downloaded from The Cancer Genome Atlas (TCGA) data portal (https://www.cancer.gov/about-nci/organization/ccg/ research/structural-genomics/tcga/using-tcga/types).The Ensembl human genome browser, GRCh38.p13 (http://asia.ensembl.org/index.html) was used to annotate and classify the lncRNAs and protein-coding genes [36]. Patient samples were excluded (n = 16) if survival times of patients were less than 30 days to eliminate non-cancer related deaths. In addition, patients with incomplete clinical data (grade stage, n = 3; AJCC stage, n = 2) were excluded from the study. Since the data was obtained from a public database, approval from the Ethics committee or written informed consent from patients was not required.

Identification of autophagy-related lncRNAs
We first identified 232 autophagy-associated genes from the Human Autophagy Database (HADb; http://www. autophagy.lu/index.html), which contains exhaustive, upto-date list of human autophagy-related genes [37]. We calculated Pearson correlation coefficients to determine the correlation between the expression of the lncRNAs and the corresponding autophagy-related genes. The autophagy-related lncRNAs were selected based on the criteria that the absolute value of correlation coefficient was greater than 0.7 (|R|>0.7) and the P value was less than 0.05 (P < 0.05).

Construction of the prognostic signature
The univariate Cox regression model was used to identify autophagy-related lncRNAs whose expression levels were significantly associated (P < 0.05) with the overall survival (OS) of the BLCA patient cohort. The hazard ratios (HRs) were used to identify risk-related lncRNAs (HR > 1) and protective lncRNAs (HR < 1). Subsequently, the candidate autophagy-related lncRNAs were subjected to multivariate Cox regression analysis to evaluate their contribution as independent prognostic factors in patient survival. Thus, we identified five target autophagy-related lncRNAs as candidates for the prognostic signature model, which was constructed based on a linear combination of the lncRNA expression levels and regression coefficients obtained from the multivariate Cox regression model. The optimal lncRNA prognostic signature was selected based on the lowest Akaike information criterion (AIC) value for further analysis. The computational formula used to determine the risk score for each patient based on this prognostic signature model was as follows: where Coef (i) and x(i) represent the estimated regression coefficient and the expression value of each autophagy-related lncRNA, respectively.

Evaluation of the prognostic signature
The BCLA patients were classified into high-risk or low-risk groups based on their prognostic risk score by using the median risk score as a cut-off point. The Kaplan-Meier survival curve and two-sided log-rank test was used to compare the overall survival (OS) of the high-and low-risk group patients. Principal component analysis (PCA) was performed to visualize gene expression patterns in the patient samples from the two groups. The receiver-operating characteristic (ROC) curves were applied to evaluate the diagnostic efficacy of each clinicopathological characteristic and the prognostic signature. Stratified survival analysis was performed to examine the accuracy of the prognostic signature in predicting patient survival outcomes. Furthermore, univariate and multivariate Cox regression analyses were performed to evaluate whether the risk score was independent of other clinical variables such as age, gender, grade, AJCC stage, T stage and N stage in determining the prognosis of the BCLA patients. M stage was not analyzed because the data was missing for several patients. P < 0.05 was considered statistically significant.

Establishment and validation of nomogram
We constructed a nomogram by integrating traditional clinical variables such as age, gender, grade, AJCC stage, T stage and N stage as well as the risk score derived from the prognostic signature to analyze the probable 1-, 3-, and 5-year OS of the BLCA patients. We then used the concordance index (C-index) to evaluate the discrimination and predictive ability of the nomogram. The range of the C-index value was 0.5 to 1.0. A higher C-index indicates greater discrimination ability of the predicting model. Furthermore, calibration curves of the nomogram were generated to examine the concordance between predicted survival and observed survival after bias correction.

Construction of the LncRNA-mRNA co-expression network
The mRNA-lncRNA co-expression network was constructed to analyze the correlation between the autophagy-related lncRNAs and their target mRNAs. Pearson correlation coefficients were calculated to identify the mRNAs that are significantly associated with their target lncRNAs based on the absolute threshold coefficient value > 0.3. The lncRNA-mRNA co-expression network was constructed and visualized using the Cytoscape software (version 3.7.2, http://www.cytoscape.org/).

Functional enrichment analysis
The lncRNA-related mRNAs were subjected to gene ontology (GO) enrichment analysis to identify the biological processes, molecular functions, and cellular components associated with the lncRNAs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to determine the main signaling pathways regulated by these lncRNAs. P < 0.05 was considered statistically significant.

Gene set enrichment analysis (GSEA)
The genome wide expression profiles of the BCLA patients were subjected to gene set enrichment analysis (GSEA; http://www.broadinstitute.org/gsea) to determine the genes that are differentially expressed between the high-and low-risk group patients [38]. The gene sets were filtered using the maximum and minimum gene set size of 500 and 15 genes, respectively. The enriched gene sets were obtained based on a P value < 0.05 and a false discovery rate (FDR) value < 0.25 after performing 1,000 permutations.

Statistical analysis
The data was processed using the PERL programming language (Version 5.30.2, http://www.perl.org). All statistical analyses were performed using the R software (version 3.6.2, https://www.r-project.org/). P < 0.05 was regarded as statistically significant.

AUTHOR CONTRIBUTIONS
SZL developed the study concept and design, performed data acquisition and analysis, and drafted the manuscript; JCY performed the bioinformatics and the statistical data analysis; XCT constructed the figures and tables; LTC was responsible for the integrity of the entire study and manuscript review. All authors read and approved the final manuscript.