A novel signature based on autophagy-related lncRNA for prognostic prediction and candidate drugs for lung adenocarcinoma.

Background: Autophagy inhibits tumorigenesis by limiting inflammation. Various lncRNAs are associated with tumour biological processes, including lung adenocarcinoma (LUAD), but the role of autophagy-related lncRNAs (ARlncRNAs) in LUAD has not been fully elucidated. Thus, this study aimed to construct a prognostic signature based on ARlncRNAs for LUAD.
Methods: The RNA-seq (FPKM) data and clinical information of LUAD patients were downloaded from The Cancer Genome Atlas (TCGA) database. After differentially expressed lncRNAs in tumour and normal groups were identified, cox regression analyses were performed to construct a prognostic signature which was then assessed through independent prognostic analysis and functional enrichment analysis. Moreover, based on the mRNAs co-expressed with the ARlncRNAs, several potential small molecule drugs were explored in the Connectivity map (Cmap).
Results: A signature consisting of seven ARlncRNAs (FAM83A-AS1, LINC01116, ILF3-DT, EBLN3P, AL161785.1, AC092279.1 and AC026355.2) was constructed to predict overall survival (OS) for LUAD. The signature was identified to be independent by the cox regression analysis and obtained the largest area under the curve (AUC =0.721) in the receiver operating characteristic (ROC). Six small molecule drugs (MS-275, methotrexate, desipramine, benzbromarone, rifampicin and doxazosin) were selected from Cmap. Conclusions: A novel ARlncRNA signature for LUAD prognostic prediction was constructed, which had better efficacy than the TNM stage and used to propose potential therapeutic regimens for LUAD patients. 2022 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is one of the world’s leading causes of cancer mortality (1), with adenocarcinoma being the most common accounting for about 40% of all cases (2). With the widespread use of low-dose computed tomography (CT), lung cancer can be diagnosed at an early stage, greatly reducing the mortality rate (3). Although the therapeutic measures for advanced lung adenocarcinoma (LUAD) have been improved and more individualised, including chemotherapy, radiotherapy, and immunotherapy (4), the prognosis remains poor. Therefore, the construction of a reliable prognostic signature is essential for the prognosis of LUAD patients and to develop effective treatment strategies.

Long non-coding RNAs (lncRNAs) are a series of transcript RNAs longer than 200 nucleotides without the functions of protein translation (5). Several studies have revealed that lncRNAs play important roles in a variety of tumour biological processes, such as neoplastic progression (6), cellular proliferation, migration, and DNA damage (7). Several lncRNAs are associated with LUAD, including FAM83A−AS1, LINC01116, MALAT1 and HOTAIR (8,9). Overexpression of FAM83A−AS1 confers tumour cells with epidermal growth factor receptor-tyrosine kinase inhibitors (EGFR-TKIs) resistance (10). Similarly, upregulation of LINC01116 markedly promoted cell proliferation and migration, inhibited cell apoptosis, and accelerated progression through the cell cycle (11). Taken together, FAM83A−AS1 and LINC01116 are involved in the development and metastasis of LUAD and may be promising biomarkers for LUAD diagnosis, prognosis, and therapy. Compared with the TNM stage, a direct factor associated with the prognosis of LUAD, lncRNAs can be extracted from serum causing no additional trauma to the patients. Thus, a better understanding of lncRNAs in LUAD may provide additional evidence for prognostic biomarkers and therapeutic targets for LUAD.

Autophagy is the degradation of intracellular lysosomes and the circulation of proteins and organelles (12) controlled by autophagy-related genes (ARGs). It has been implicated in several disease processes, including cancer, neurodegenerative disorders, autoimmune diseases, cell death, and cardiovascular diseases (13-16). In cancer development, autophagy is a double-edged sword, functioning as a tumour suppressor in the early stages, while promoting tumorigenesis and causing resistance to therapeutic agents in the advanced stages (17). Autophagy is regulated by several cellular molecules and signalling pathways, including ARGs and lncRNAs. Therefore, this study aimed to construct a risk model based on autophagy-related lncRNAs that could be utilised for clinical prediction and drug selection.

We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-21-1554/rc).

Methods

Data acquisition and processing

The LUAD RNA-seq (FPKM) data, tumour mutation burden (TMB), and corresponding clinical information were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The cohort contained 497 tumour tissues and 54 normal tissues, with the complete clinical information of 486 patients (tumour =439, normal =47) extracted as an entire set for further analysis. To increase the study reliability, the entire set (only tumour tissues) was randomly divided into a training set (accounting for 70%, n=307) and a validation set (accounting for 30%, n=132). The “rtracklayer” and “dplyr” R packages were applied to convert gene names from Ensemble IDs to a profile of gene symbols with the Ensemble database.

Identification of autophagy-related lncRNAs

Two hundred and thirty-two ARGs were extracted from the human autophagy database (HADb) analysis. Pearson’s correlation analysis of the expression of lncRNAs and ARGs was performed through the “limma” R package and 312 lncRNAs with high correlation (|R| >0.3, P<0.001) with ARGs were identified as autophagy-related lncRNAs (ARlncRNAs).

Identification of an ARlncRNAs prognostic signature

The tumour and normal tissues were screened for the differentially expressed ARlncRNAs via the “limma” R package with thresholds set at |logFC| >1 and P<0.05. Then, to identify ARlncRNAs associated with survival, univariate cox regression analysis was performed according to the criteria of P<0.01 using the “survival” R package. Subsequently, multivariate cox regression analysis was conducted to construct the prognostic signature based on the lowest Akaike information criterion (AIC). The risk score formula is as follows: Risk score = ∑Coef (ARlncRNAs) * Exp (ARlncRNAs). The Coef (ARlncRNAs) represent the coefficient of each ARlncRNA and Exp (ARlncRNA) is the expression of each ARlncRNA. Based on the median risk score, the patients were divided into high-risk and low-risk groups, with Kaplan-Meier survival analysis performed to estimate the survival rate of both groups using the “survival” and “survminer” R packages. The receiver operating characteristic (ROC) curve was used to assess the prognostic performance through the “timeROC” R package.

Independent prognostic analysis and ROC curve

To assess the relationship of survival prognosis with clinicopathological factors (age, gender, AJCC stage, TNM stage, tobacco history, and anatomical location) and risk score, univariate and multivariate cox regression analyses were performed, respectively, using the “survival” R package. Furthermore, a time-dependent ROC curve was used to evaluate the predictive accuracy for survival time by different clinicopathological factors and risk scores using the “survivalROC” R package.

Construction and assessment of a nomogram

A nomogram incorporating the lncRNA signature and clinical factors was created using the “rms” R package to simplify the predictive model. Calibration curves were also plotted to assess whether the predicted survival in the nomogram agreed with the actual survival.

Principal components analysis (PCA) and gene set enrichment analysis (GSEA)

PCA was used to test the differentiation of patients into low- and high-risk groups, with GSEA performed to evaluate different functional phenotypes between low- and high-risk groups.

Construction of the lncRNA-mRNA co-expression network and functional enrichment analysis

A co-expression network of the autophagy-related lncRNA-mRNAs with prognostic value was established according to the criteria of |R| >0.3 and P<0.001 by Pearson correlation analysis using the “limma” R package. To further explore the functional annotation and pathway analysis of the target mRNAs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler” R package with the criteria of P<0.05 and FDR <0.05.

Identification of potential immunotherapy

To explore the value of our signature for predicting the efficacy of immunotherapy in LUAD, immune cell infiltration was first estimated using CIBERSORT. Then, we assessed the TMB of LUAD patients in the high-risk and low-risk groups. The Cmap database (18) uses gene-expression signatures to predict small molecular compounds for a specific disease. In the present study, the target mRNAs were divided into two groups, up- and downregulated target genes, which were uploaded from the lncRNA-mRNAs network to Cmap. A connectivity score ranging from −1 to 1 was used to reflect the degree of closeness between the expression spectrums. Drugs with negative scores were potential therapeutic molecules, so were further investigated in the PubChem database (19).

Statistical analysis

All statistical analyses and plotting were performed in the R language (Version 3.6.3). Univariate and multivariate cox regression analyses were used to identify independent prognostic factors for LUAD, through the criteria of P<0.05 and the lowest AIC respectively. Kaplan-Meier survival analysis was performed in the high-risk and low-risk groups and evaluated using a log-rank test. A two-sided P value <0.05 was considered statistically significant.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Identification of ARlncRNAs with significant prognostic value in LUAD

The flow chart is shown in . First, a total of 14,064 lncRNAs were extracted from the TCGA dataset, 1,639 of which were identified as ARlncRNAs by Pearson correlation analysis (|R| >0.3, P<0.001). After differential expression analysis of the entire set, 312 ARlncRNAs were selected as shown in . Among them, 10 ARlncRNAs were significantly associated with the survival of LUAD patients from the TCGA (P<0.01) by univariate cox regression analysis, namely, ILF3-DT, EBLN3P, AL161785.1, FAM83A-AS1, AL049555.1, CRNDE, CH17-340M24.3, LINC01116, AC092279.1, and AC026355.2. Subsequently, multivariate cox regression analysis further identified seven lncRNAs with prognostic significance based on the lowest AIC (1,037.67) (). Among them, FAM83A−AS1 and LINC01116 were considered as risk factors with a hazard ratio (HR) greater than 1, whereas ILF3-DT, EBLN3P, AL161785.1, AC092279.1, and AC026355.2 were considered as protective factors with HR values less than 1 (). The risk score for each patient was calculated using the following formula for the ARlncRNAs signature: risk score = (–0.039 * expression level of ILF3-DT) + (–0.038 * expression level of EBLN3P) + (–0.079 * expression level of AL161785.1) + (0.033 * expression level of FAM83A-AS1) + (0.069 * expression level of LINC01116) + (–0.165 * expression level of AC092279.1) + (–0.225 * expression level of AC026355.2). Next, based on the risk score formula and the calculated median risk score, the patients in the training cohort (n=307) and validation cohort (n=132) were divided into high-risk and low-risk groups. The risk distribution, survival status, and gene expression pattern of the two cohorts are shown in . The heatmap () showed that two ARlncRNAs are highly expressed in the high-risk group, with the other five ARlncRNAs highly expressed in the low-risk group. The risk curve and scatterplot were applied to illustrate the risk score and the relevant survival status, indicating that the mortality occurrence depended on the risk score (). The Kaplan-Meier curve demonstrated that the OS of the high-risk group was significantly poorer than that of the low-risk group (P=3.283e−05 in the training cohort and P=1.044e-03 in the validation cohort; ), indicating that the risk score has prognostic value. PCA was employed to demonstrate the distribution difference between the high- and low-risk groups based on the risk model (), with satisfactory results. Time-dependent ROC curves revealed that the signature performed well in predicting OS. The AUC at 1-, 3-, and 5-year in the training cohort was 0.744, 0.679, and 0.698, and 0.671, 0.82, and 0.857 in the validation cohort ().

Figure 1

The main flowchart of this study. LUAD, lung adenocarcinoma; ARlncRNAs, autophagy-related long non-coding RNAs; DEARlncRNAs, differentially expressed ARlncRNAs; K-M analysis, Kaplan-Meier survival analysis; TMB, tumour mutation burden; PCA, principal components analysis.

Figure 2

Screening of the 7 ARlncRNAs. (A) Volcano plot of the 312 differentially expressed ARlncRNAs (|logFC|>1, FDR <0.05). (B) Multivariate cox regression analysis results with the lowest AIC (AIC =1,037.67, C-index =0.7). (C) The K-M survival curves of 7 prognostic ARlncRNAs. Five ARlncRNAs (ILF3-DT, EBLN3P, AL161785.1, AC092279.1 and AC026355.2) were protective factors and the other two (FAM83A−AS1 and LINC01116) were risk factors.

Figure 3

The validation of the 7 ARlncRNAs prognostic signature in the training and validation cohort. (A,B) Heatmap of the expression of the 7 ARlncRNAs in the high- and low-risk groups; distribution of risk scores between the two groups; distribution of survival status and OS time of each sample in the training cohort (A) and validation cohort (B). (C,D) K-M survival curves of overall survival in the training cohort (C; P=3.283e-05) and validation cohort (D; P=1.044e-03). (E,G) ROC curves of the prognostic signature for predicting 1-, 3- and 5-year OS in the training cohort; (E) AUC of 1-, 3-, and 5-year survival were 0.744, 0.679, and 0.698) and validation cohort; (G) AUC of 1-, 3-, and 5-year survival were 0.671, 0.82, and 0.857). (F,H) PCA based on the risk model of the 7 ARlncRNAs in the training cohort (F) and valiation cohort (H).

Evaluation of the risk model of seven ARlncRNAs as an independent prognostic factor

Univariate and multivariate cox regression analyses were performed in the training set to determine whether the signature was an independent prognostic factor, showing that AJCC stage (P<0.001), T stage (P=0.002), N stage (P<0.001), and risk score (P<0.001) were significantly related to OS in univariate cox regression analysis (), with only the AJCC stage (P<0.001) and risk score (P<0.001) significantly related to OS in multivariate cox regression analysis (). As shown in , the AUC of risk score was 0.721, which was higher than the AUC of the AJCC stage, T stage, and N stage, indicating that the prognostic risk model was relatively reliable. Taken together, the risk model of seven ARlncRNAs is an independent prognostic factor for LUAD patients.

Figure 4

Evaluation of prognostic accuracy of the risk model and other clinicopathological characteristics. (A) Univariate and (B) multivariate cox regression analyses of risk score and other clinical parameters (age, gender, AJCC stage, anatomic location, tobacco history, T, N and M). (C) ROC curves of risk score (AUC =0.721) and other clinical features.

Correlation of the ARlncRNAs prognosis signature to other clinicopathological characteristics

To further evaluate whether the seven ARlncRNAs were involved in the occurrence and development of LUAD, we explored the relationship between the risk score and clinicopathological factors (), showing that patients without lymphatic metastasis had a lower risk than those with lymphatic metastasis (P<0.05).

Table 1

Correlation of the ARlncRNAs prognosis signature with other clinicopathological characteristics

Clinical	Group	n	Mean	SD	t	P
Age	≤65	138	1.452	2.518	1.525013	0.129
	>65	149	1.168	0.793
Gender	Female	164	1.168	0.977	−1.5607	0.12
	Male	123	1.479	2.548
Stage	Stage I-II	223	1.167	0.835	−1.67197	0.098
	Stage III-IV	64	1.85	3.721
Anatomic	Left	122	1.467	2.646	1.219104	0.224
	Right	165	1.202	1.041
Tobacco	Never	46	1.153	0.53	−1.4186	0.157
	Smoke	241	1.333	1.992
T	T1-2	256	1.189	0.833	−1.4276	0.16
	T3-4	31	2.175	4.82
N	N0	195	1.106	0.788	−2.33317	0.021
	N1-3	92	1.749	3.069
M	M0	273	1.284	1.812	−0.84766	0.409
	M1	92	1.848	2.716

Evaluation of the prognostic prediction nomogram

The nomogram visualises the regression analyses results and is commonly used by clinicians to predict the risk of disease or survival (20). A nomogram was constructed to accurately evaluate 1-, 3-, and 5-year survival probabilities using the risk score calculated from the ARlncRNAs prognostic signature and other clinicopathological characteristics, including age, gender, anatomical location, AJCC stage, tobacco history, T stage, N stage, and M stage (). The calibration plot analysis showed a high degree of consistency between the actual and the predicted 1-, 3-, and 5-year OS when compared to the reference line (), suggesting that the nomogram was accurate and reliable.

Figure 5

Construction and assessment of a nomogram based on the risk model. (A) The nomogram of 1-, 3- and 5-year OS based on the risk model and other clinical features. (B-D) Calibration plots used for evaluating the consistency between the actual and the predicted 1-, 3-, and 5-year OS. OS, overall survival.

Gene set enrichment analysis

GSEA was used to conduct the functional annotation revealing that several autophagy-related signalling pathways were significantly enriched, such as “endocytosis”, the “ErbB signalling pathway”, “Notch signalling pathway”, “p53 signalling pathway”, “pathway in cancer”, “TGF-β signalling pathway” etc. ().

Figure 6

Gene Set Enrichment Analysis (GSEA) of the two groups based on the risk model. (A) The red box shows several significant pathways in the high-risk group. (B) The blue box shows enrichment results in the low-risk group.

To explore the regulation of the ARlncRNAs, a co-expression network consisting of seven ARlncRNAs and 51 ARGs was constructed using the Cytoscape software (). The Sankey diagram showed the correlation between the 51 ARGs and 7 ARlncRNAs (). To further investigate the function of the network, GO and KEGG analyses were performed. The top GO terms were “autophagy” in biological process (BP), “plasma membrane-bounded cell projection cytoplasm” in cellular component (CC), and “cysteine-type peptidase activity” in molecular function (MF), respectively (). According to KEGG analysis, “autophagy”, the “longevity regulating pathway”, “FoxO signalling pathway”, “HIF-1 signalling pathway”, “NOD-like receptor signalling pathway”, “PD-L1 expression and PD-1 checkpoint pathway” were the main pathways (). Taken together, the regulatory network may lead to the pathogenesis of LUAD through multiple signalling pathways.

Figure 7

Construction of the lncRNA–mRNA co-expression network and functional enrichment analysis. (A) LncRNA–mRNA co-expression network: grey represents lncRNAs and yellow represents mRNAs. (B) Sankey diagram shows the connection between the 51 mRNAs and 7 ARlncRNAs. (C,D) Functional enrichment analysis of target mRNAs. (C) Gene Ontology (GO) enrichment analysis. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

Immune cell infiltration was estimated in LUAD patients using CIBERSORT to explore the possible relationship between the signature and potential immunotherapy, showing that among the 22 types of immune cells, plasma cells, resting dendritic cells, and resting mast cells had significant negative correlations with the risk score, while M0 macrophages and activated mast cells were positively correlated with the risk score (), implying the potential of the signature to reflect the tumour immune microenvironment (TIME). Moreover, there was a significant difference in the expression of TMB between the high- and low-risk groups (), indicating that immunotherapy may be a potentially effective treatment based on the signature. The 31 upregulated and 20 downregulated ARGs were analysed in Cmap and the results with P<0.05 were regarded as potential small molecule drugs (). The chemical structures of the top six agents, namely, MS-275 (C21H20N4O3), methotrexate (C20H22N8O5), desipramine (C18H22N2), benzbromarone (C17H12Br2O3), rifampicin (C43H58N4O12), and doxazosin (C23H25N5O5) are shown in .

Figure 8

Screening of 6 LUAD candidate small molecule drugs. (A) Correlations of the signature with immune cell infiltration (The P values of plasma cells, macrophages M0, dendritic cells resting, mast cells resting, and mast cells activated are less than 0.05). (B) The correlation between risk score and TMB (P<0.001). (C) The top six most significant small molecule drugs based on Cmap (MS-275, methotrexate, desipramine, benzbromarone, rifampicin, and doxazosin). (D) Chemical structure depiction of the six potential small molecule drugs for LUAD.

Discussion

Lung cancer is the leading cause of cancer morbidity and mortality in the world (1). For early stage LUAD patients, surgery is the recommended treatment (21), whereas chemotherapy, radiotherapy, immunotherapy, and targeted therapy are recommended for advanced LUAD patients (22). Several studies have demonstrated a strong association between TMB and the clinical benefits of receiving immunotherapy (23,24). Moreover, TIME is closely correlated to the efficacy of immunotherapy (25,26), suggesting that immunotherapy may be more effective for high-risk patients based on the risk model.

Autophagy plays a dual role at different phases of cancer development (17). Recently, regulating autophagy to improve cancer therapeutic effects has attracted great interest (27), thus, an autophagy-related signature could be used as a prognostic biomarker and to predict potential small molecular medicine for LUAD patients.

In the present study, we constructed an ARlncRNAs prognostic signature consisting of seven ARlncRNAs, in which, FAM83A−AS1 was upregulated in LUAD tumour tissues and could promote LUAD cell migration and invasion (10) and LINC01116 was associated with gefitinib and cisplatin resistance in LUAD (28,29), indicating that they could be novel biomarkers for LUAD diagnosis, prognosis, and therapy. The risk score was used to classify the training and validation cohorts into high- and low-risk groups, with the low-risk group having a better survival outcome than the high-risk group. The signature was shown to be an independent prognostic factor and more powerful in survival prediction than other clinical indicators, thus having great potential for future clinical application.

Cancer initiation and progression is regulated by multiple signalling pathways, including the ErbB signalling pathway. The erbB family of receptor tyrosine kinases include HER1 (EGFR), HER2 (erbB2), HER3 (erbB3), and HER4 (erbB4). The HER2 gene is associated with breast cancer and EGFR is highly expressed in more than 60% of NSCLCs, regulating the proliferation, survival, and motility of the tumour cells. Thus, inhibitors that target the kinase domain of EGFR have been developed and are clinically active (30). The GSEA results showed that the signature was strongly associated with the “erbB signalling pathway”, which could guide the follow-up treatment of LUAD patients.

There are usually co-expression relationships between genes, so a lncRNA-mRNA co-expression network was established for functional enrichment analysis to explore the biological mechanisms of the seven selected ARlncRNAs. Cellular senescence was enriched in the network and has been shown to restrain the cell proliferation and transformation induced by oncogene activation (31). Autophagy was also demonstrated to be necessary for the establishment of cellular senescence, so, it might be a promising therapy for LUAD by targeting the regulation of autophagy.

Importantly, we confirmed that the signature could potentially reflect the TIME and there was a significant difference in the expression of TMB between the two groups. Six small molecular drugs, namely, MS-275, methotrexate, desipramine, benzbromarone, rifampicin and doxazosin were identified from the prediction of Cmap based on the up- and downregulated ARGs. MS-275, also called entinostat, the histone deacetylase inhibitor (HDACi), is a promising drug for the treatment of breast cancer in different clinical settings. Specifically, entinostat may offer better safety through targeting class I HDACs (HDAC 1 and 3) than those nonselective broad-spectrum HDAC inhibitors. So far, there have no HDAC inhibitors been approved for NSCLC, however, entinostat alone or in combination with other molecular targeted therapies for NSCLC may be an option (32-34). Methotrexate (MTX), an anti-folate agent, is one of the most widely used and effective drugs for treating cancer, including lung cancer. However, high doses of MTX usually lead to serious side effects. Nowadays, some carriers have been loaded with MTX to improve drug safety (35-37). Desipramine, an antidepressant, is commonly used in combination regimens, not only helping to relieve depression but also enhancing the cytotoxicity of the platinum agents (38). Benzbromarone, a nonpurine xanthine oxidase inhibitor for the treatment of gout, is promising for the treatment of tumour metastasis, angiogenesis, and angiopathy (39,40). Rifampicin, a common antibiotic used for tuberculosis exerts an anti-tumour effect through direct antitumor activities when highly concentrated and antiangiogenic properties (41). Doxazosin is an α1-adrenoceptor antagonist used to treat hypertension or benign prostatic hyperplasia. Moreover, the potential of doxazosin as an osimertinib sensitiser affecting many cancers via the upregulation of autophagy was also demonstrated (42,43). Taken together, these findings may provide potential therapeutic options for LUAD patients.

Recently, the prognostic signatures of ARlncRNAs were established in several cancers (44-46), however, the role of ARlncRNAs in the prognosis of LUAD remains unclear. Thus, we focused on the essential role of ARlncRNAs in tumour malignancy to construct an ARlncRNA signature for the LUAD cohort, demonstrating that it has independent prognostic value, as well as exploring the clinical utility in potential immunotherapy and targeted therapy for LUAD patients.

Nonetheless, the study has several limitations. First, the ARlncRNAs prognosis signature was constructed based on the publicly available TCGA database. Secondly, our findings need to be proved in other independent cohorts to improve the predictive efficacy. Moreover, biochemical experiments such as immunohistochemistry, quantitative real-time PCR, and flow cytometry need to be conducted to confirm the findings.

Conclusions

A prognostic signature of seven ARlncRNAs was constructed for the LUAD cohort to explore the role of autophagy in tumour malignancy, with various potential therapeutic regimens proposed based on the immune cell infiltration and differential expression of TMB.