Yu-Feng Huang,Jing Xiang,Chun Yuan Chen,Biao Deng,Xia Yang,Lian-Xiang Luo,Zhu Liang*

1The Graduate School,Guangdong Medical University,Zhanjiang 524023,Guangdong,China.2Department of Cardiothoracic Surgery,Affiliated Hospital of Guangdong Medical University,Zhanjiang 524000,Guangdong,China.3Southern Marine Science and Engineering Guangdong Laboratory(Zhanjiang),Zhanjiang 524023,Guangdong,China.4The Marine Biomedical Research Institute,Guangdong Medical University,Zhanjiang 5240235,Guangdong,China.5The Marine Biomedical Research Institute of Guangdong Zhanjiang,Zhanjiang 524023,Guangdong,China.

Abstract Background:Autophagy-associated long non-coding RNAs(aalncRNAs)take an important position in the tumorigeness of lung cancer,but current researches have not systematically investigated autophagy-associated lncRNAs in lung adenocarcinoma(LUAD).Methods:In this research,RNA-sequences of LUAD patients were downloaded from the TCGA and autophagy-associated genes were obtained from the GSEA website.The Pearson's test was conducted to find the correlation between autophagy-associated lncRNAs and autophagy-associated genes.AalncRNAs with prognostic significance were identified by using Cox and LASSO regression analysis in R,gradually.Risk score model was built to estimate prognosis-associated lncRNAs.Results:A risk score model was established according to the expressions of 7 aalncRNAs(RP11-102K13.5,RP11-1029J19.4,LINC00942,KLHL7-AS1,AC092198.1,C20orf197,LINC01116),and low-risk group was found to have a better prognosis(P＜0.001).Next,single gene expression survival analysis show that 4 out of these lncRNAs were significantly associated with the survival of patients.In addition,the AUC value of model reached 0.724,demonstrating the good predictive ability of the model.Conclusion:These aalncRNAs in LUAD might possibly offered biological markers for the diagnosis and therapy of lung adenocarcinoma.

Keywords:lung adenocarcinoma,lncRNA,autophagy,nomogram,TCGA

Introduction

Cancer is the biggest problem for human health and life and a major contributor to the cause of death worldwide[1].Among the different types of cancer,lung cancer ranks first in terms of incidence and mortality[2].About 2.1 million people are diagnosed annually by this cancer,accounting for 1.76 million deaths and its 5-year survival rate is no more than 15%[2].Non-small cell lung cancer(NSCLC)is the most predominant biological type of lung cancer,accounting for more than 80% of all patients,and lung adenocarcinoma(LUAD)and lung squamous cell carcinoma(LUSC)are the most prevalent categories of NSCLC[3].Therefore,it is of great significance to understand the developmental and molecular mechanisms of patients with LUAD and identify potential biological markers to improve the prognosis of patients.

Autophagy is an important course of action involving the degradation of old organelles and proteins to obtain energy in eukaryotic cells and plays important role in many human diseases such as autoimmune diseases,nervous system defects,inflammatory diseases,and many kinds of cancers[4–6].Especially,autophagy is proving to be increasingly relevant in the development of LUAD.Therefore,it is very important to establish an autophagy-associated signature to forecast the prognosis of LUAD.

Long non-coding RNAs(lncRNAs)are RNAs more than 200 bp of nucleotides without protein-coding abilities[7].However,lncRNAs can moderate gene expression at the epigenetic,transcriptional and post-transcriptional levels by interacting with protein,chromatin and RNA targets[8,9].Growing evidence indicates that aberrant expression of lncRNAs may cause adverse effects on biological processes and lead to many diseases particularly different kinds of human diseases,particularly cancer[10–12].Resent study demonstrated that different expression of lncRNAs in lung cancer can promote tumor cell growth,apoptosis,invasion,and metastasis[13].Therefore,it is of great importance to discriminate the autophagy lncRNA correlated with patient prognosis.

More interestingly,an increasing number of studies have reported that autophagy-associated lncRNAs play an important role in the development of many types of cancer[14–25],especially in tumor resistance.For example,the lncRNA PVT1 promotes gemcitabine resistance in pancreatic cancer by activating Wnt/β-catenin and regulating the autophagic pathway.In gastric cancer,lncRNA CRNDE is a protective aalncRNA that can attenuate chemoresistance in gastric cancer through SRSF6-regulated PICALM alternative splicing;in contrast,aalncRNAs such as MALAT1,ARHGAP5-AS1,and EIF3J-DT can induce chemoresistance in gastric cancer through autophagy.In gallbladder cancer,GBCDRlnc1 induces chemoresistance in cancer cells by activating autophagy.In breast cancer,the lncRNA OTUD6B-AS1 promotes paclitaxel resistance in triple-negative breast cancer by regulating autophagy and genomic instability mediated by the miR-26a-5p/MTDH pathway.In diffuse large B-cell lymphoma,MALAT-1 regulates autophagy-related signaling pathways to induce its chemoresistance.In addition,in lung cancer,aalncRNAs ACTA2-AS1[26],LINC01559[20],SNHG7[27],HIF1A-AS2[28],MIR99AHG[29],APCDD1L-AS1[30],LOC389641[31],LINC00857[32],KTN1-AS1[33],PANDAR[34],BLACAT1[35],NBAT1[36]and LCPAT1[37]can increase or inhibit autophagy to promote or inhibit tumor progression,of which ACTA2-AS1,SNHG7,HIF1A-AS2,APCDD1L-AS1 and BLACAT1 play an important role in lung cancer treatment resistance.This suggests that aalncRNA is very important in the diagnosis and treatment of lung cancer,especially in terms of tumor drug-resistance.This has led to a great interest in the development of new aalncRNAs for the diagnosis and treatment of LUAD.

Consequently,we attempt to develop an aalncRNAs system to predict patient prognosis and provide a theoretical groundwork for the diagnosis and medical therapy of LUAD.

Methods and materials

Collection and downloads of data

The flow chart for our research was displayed in Figure 1.We downloaded 594 patients'RNA sequence data(tumor:535 cases;normal:59 cases)and clinical information from the TCGA(https://cancergenome.nih.gov/).Patients without comprehensive clinical details including gender,age,TNM stage,status as well as survival time less than 30 days were excluded.After removing patients with incomplete clinical information,327 patients were reserved in this research for further analysis.

Figure 1 Flow chart.Flowchart for establishing and evaluating the prognostic model

Identification of autophagy-associated lncRNA in LUAD

To study autophagy-associated genes of LUAD patients,the differentially expressed genes(DEGs)between tumor and normal samples were conducted by utilizing the Bioconductor“edgeR”package(http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)as described in the R environment(version 4.1.1,http://www.r-project.org).DEGs with absolute|log2FC|＞2 and FDR＜0.01 were considered for inclusion in further analysis.Heat maps were produced using“pheatmap”R package in the R environment to present the results of differential expression analysis as previous described.Then,in order to obtain autophagy-associated gene,the autophagy-associated gene set was identified from the MSigDB public online databases(http://www.gsea-msigdb.org/gsea/msigdb/search.jsp)by searching with key words“autophagy”through the"Search Gene Sets"tool,and the autophagy-associated DEGs were isolated by the Venn Diagram R package[38].The correlation coefficient(R2)between aalncRNAs and DEGs was calculate by using Pearson correlation test in the R environment as previous described[39].The|R2|＞0.3 andP＜0.05 was the criteria for screening aalncRNAs.

Establishment of prognostic autophagy-associated lncRNAs signature

The aalncRNAs with prognostic significance were detected by Univariate Cox regression analysis using"survival"R package in the R environment,andP＜0.05 were considered to be prognostically relevant lncRNAs.The"survival"package in R is described in more detail by Christensen and at https://cran.r-project.org/web/packages/survival/survival.pdf.Then,in order to clarify the essential aalncRNAs that were notably related to overall survival,LASSO regression analysis was used to further filtering for the above lncRNAs by using the“glmnet”R package in the R environment[40].Lastly,in order to screen more important signatures,the multivariate Cox regression model was used for the computation of correlation coefficient of these aalncRNAs by using"survival"R package in the R environment.The model including coefficients and lncRNA expressions were applied to obtain risk scores for all patient samples.The calculation formula of risk score is as follows(coef=coefficient,E=expression of aalncRNA):

Next,patients were separated into high-and low-risk groups according to the median risk score.More importantly,survival curves were performed to make a comparison of the prognosis between these groups.Furthermore,the risk scores and clinical variables were analyzed by univariate and multivariate Cox regression analysis by using"survival"R package in the R environment in order to evaluate their prognostic implications.Finally,the receiver-operating characteristic(ROC)curves and area under the curve(AUC)were built up to assess the accuracy and verify the model by using“survivalROC”R package in the R environment[41].

Construction and assessment of prognostic models

The“rms”R package in the R environment was utilized to constructed the nomogram which contained risk scores and clinical factors(age,gender and TMN stage).The nomogram can be used to forecast the possible survival rates of LUAD patients at 1-,3-and 5-year.The receiver-operating characteristic(ROC)curves and area under the curve(AUC)at 1-,3-and 5-year were built up to assess the accuracy and verify by using“survivalROC”R package in the R environment[41].Calibration curve for nomogram was plotted using the R package“survival”in the R environment to visualize the predictive power of nomogram.

Signature-based analysis of biological processes and pathways

Gene Set Enrichment Analysis(GSEA)is a gene set-based enrichment analysis method that can determine the effect of synergistic changes in genes within this gene set on phenotypic changes[42],so we can find the relationship between the model and the signaling pathway.GSEA was conducted with java software GSEA(http://www.gsea-msigdb.org/gsea/downloads.jsp)by using GSEA v.4.0.1 graphical user interface and the hallmarks gene set from MSigDB v7.0,with the following parameters:Signal2Noise metric,weighted scoring,and n=1000 permutations.GO and KEGG enrichment analysis on the differential risk score was analyzed to extrapolate their functions.The data generated by the GSEA software was exported and further visualization were performed in the R software environment.

Results

Differentially expressed mRNAs,lncRNAs co-expression network and differentially expressed lncRNAs

After Data Screening,we enrolled 327 patients for further analysis.The clinical relevance of these LUAD patients was displayed in Table 1.The RNAs expression levels in LUAD samples from the TCGA database were analyzed to identify DEGs may be associated to autophagy in LUAD patients.Firstly,Differential expression analysis between the tumor group and normal group(P＜0.01and|log FC|＞2)identified 44 DEGs(Figure 2a,b)by using“edgeR”R package.Secondly,351 aalncRNAs were abstracted from 1186 lncRNAs.Then,88 different expression aalncRNAs were identified via“edgeR”R package(Figure 2c,d).

Identification of differentially expressed lncRNAs relevant to prognosis in LUAD patients

Firstly,twenty-nine aalncRNAs were identified by the univariate Cox regression analysis.Next,eighteen aalncRNAs were identified by lasso analysis(Figure 3a,b).Lastly,seven aalncRNAs including LINC01116,C20orf197,AC092198.1,KLHL7-AS1,LINC00942,RP11-1029J19.4,RP11-102K13.5,were identified for the establishment of the prognostic model by multivariate Cox regression analysis.The regression coefficients of LINC01116,LINC00942,RP11-1029J19.4,AC092198.1 and RP11-102K13.5 were positive,while the regression coefficients of C20orf197 and KLHL7-AS1 were negative(Table 2).

Table 1 Baseline of the enrolled patients

Table 2 Multivariate Cox regression analysis of characteristics and risk score in LUAD

The effect of autophagy-associated lncRNA on prognosis of LUAD patients

Next,we constructed a signature of these 7 aalncRNAs using a risk score method.The formula for the calculation of the risk score was:riskscore=(0.164031216312927×LINC01116)+(–0.143388804207631×C20orf197)+(0.0789216780216484×AC092198.1)+(–0.160096643362271×KLHL7-AS1)+(0.100419363560448×LINC00942)+(0.148754007261771×RP11-1029J19.4)+(0.154022073635007×RP11-102K13.5).After establishing the risk model,patients with LUAD were divided into low-risk and high-risk group depending on the median risk to pinpoint the difference between these groups.Kaplan-Meier survival analysis were then performed according to a log-rank test,and the prognosis was found to be worse in the high-risk group(Figure 4a).Meanwhile,patients were separated into two groups based on the median expression according to the expression of each aalncRNAs,respectively.We found that high expression of each of these 4 aalncRNAs(RP11-102K13.5(P＜0.0001),LINC00942(P=0.00063),LINC01116(P=0.0005)and AC092198.1(P=0.027))was associated with worse OS in patients with LUAD(Figure 4b–e),all of them indicate a poor prognosis,suggesting that they are oncogenes in LUAD.

Besides,scatterplot was set up to exhibit the survival status and risk score of LUAD patients(Figure 5a,b).Differentially expressed aalncRNAs between the groups were presented in a heat map(Figure 5c).LINC01116,LINC00942,RP11-1029J19.4,AC092198.1 and RP11-102K13.5 were risk elements and their expression was increased in the high-risk group,while C20orf197 and KLHL7-AS1 were protective factors and their expression was reduced in high-risk group.

Clinical evaluation by risk assessment model

The univariate and multivariate Cox regression analysis were conducted to assess the prognostic implications of these lncRNAs signature.As a result,the 7-lncRNA signature was found to be an independent prognostic factor of LUAD patients(Figure 6a,b).Besides,we integrated risk scores and clinicopathological features,included age,gender,tumor stage,and TNM stage into the model.Multiple ROC curves consisting of the risk score and clinicopathologic features showed that the AUC value for this lncRNA signature was 0.724,which was higher than the AUCs of clinicopathologic factors(Figure 6c).

Nomogram building

According to the patients’risk scores and clinical features,a comprehensive prognostic nomogram was built to evaluate LUAD patients’1-,3-and 5-year survival probability.Seven clinical characteristics were integrated into the nomogram(Figure 7a).And the result showed that the AUCs of the 1-,3-and 5-year OS were respectively 0.724,0.731 and 0.719(Figure 7b).The calibration plots showed strong concordance between the predictions and actual observations of the nomograms in terms of survival rates at 1-,3-and 5-year(Figure 7c–e).

Functional enrichment analysis based on the risk score

We applied GSEA to the aalncRNAs prognosis signature depended on the risk score.The GO functional enrichment analysis suggested that these genes were enriched in antigen processing and presentation ofendogenous antigen,bicarbonate transport,chromosome segregation,DNA dependent DNA replication,mitotic nuclear division,mitotic sister chromatid segregation,sister chromatid segregation,GOLGI cis cisterna,MHC protein complex and T cell receptor complex(Figure 8a).KEGG enrichment analysis demonstrated that these genes were implicated in the following pathways:the allograft rejection,asthma,autoimmune thyroid disease,cell cycle,DNA replication,intestinal immune network for IgA production,pentose phosphate pathway,proteasome,ribosome and spliceosome pathway.What we found had the potential to help researchers further explore the pathogenesis of aalncRNAs in LUAD(Figure 8b).

Figure 2 Identification of autophagy differentially expressed genes(DEGs)and DElncRNAs.

Figure 3 Autophagy-associated lncRNA selection utilizing Lasso model.

Figure 4 Survival curve of patients with LUAD in different groups.

Figure 5 Construction of the autophagy-associated lncRNAs prognostic signature model.

Figure 6 Independent prognostic ability evaluation and functional enrichment analysis for the constructed prognostic signature in LUAD.

Figure 7 Construction and evaluation of the nomogram.

Figure 8 Gene set enrichment analysis.

Discussion

The incidence and mortality rate of lung cancer remains high,and the symptoms of early-stage lung cancer patients are not obvious,and many patients are diagnosed to be in advanced stage and have missed the opportunity of surgery[43,44].Chemotherapy,radiotherapy,immunotherapy and other treatment modalities improve the OS and DFS of patients to varying degrees,while targeted therapies such as immunotherapy are more easily accepted due to less toxic side effects[45],but immunotherapy does not work for every patient,and those patients who do work face the risk of tumor resistance and recurrence.Therefore,the development of new diagnostic and therapeutic targets for lung cancer is still a challenge that needs to be resolved.

Autophagy is a degradation process that occurs in eukaryotic cells to maintain cellular homeostasis under the regulation of autophagy-related genes in response to harmful stimuli such as hypoxia,lack of nutrients and drugs.Studies have shown that autophagy is closely associated with tumors,including lung,breast,liver,and colon cancers etc[46–49].Furthermore,recent studies have found that lncRNAs,as a member of non-coding RNAs,can regulate autophagy in tumors,and thereby inhibiting or promoting tumor development.For example,LncRNA SNHG11 facilitate tumor advancement of gastric cancer by promoting autophagy and activating the Wnt/β-Catenin pathway[50].PVT1 triggers cytoprotective autophagy and promotes pancreatic cancer by sponging miR-20a-5p to affect the expression of downstream ULK1[14].PANDAR suppresses autophagy and apoptosis in NSCLC by regulating the expression of BECN1 and apoptosis and inhibits lung cancer development[34].RNA H19 promotes tamoxifen resistance in breast cancer by activating autophagy through the H19/SAHH/DNMT3B axis[51].GAS5 can inhibit colorectal cancer cell migration and invasion and promote autophagy by targeting miR-222-3p through the GAS5/PTEN signaling pathway[52].It is easy to see that the continued development of sequencing technology has revealed more tumor biomarkers and potential therapeutic targets,providing a powerful tool in the human battle of anti-tumor.However,no systematic and definitive procedure has been developed for lncRNA prediction of LUAD prognosis.Therefore,it is essential to establish lncRNA signature genes for LUAD prognosis.

In this study,we obtained RNA-sequence of LUAD from TCGA database,then acquired autophagy-associated genes from GSEA database,and received differential genes and lncRNAs after differential expression analysis,and then constructed a co-expression network of lncRNAs and autophagy genes by Pearson correlation test.After regression analysis,we established a prognostic risk model of LUAD consisting of 7 aalncRNAs.We calculated the risk scores of the patients based on the coefficients of these 7 aalncRNAs,and then divided the patients into high-risk and low-risk groups based on the median risk scores,and the result showed that patients in the low-risk group were relevant to better prognosis.Then,we constructed nomogram with risk scores and clinicopathological features of patients.The results of ROC curves and calibration curves indicated the robustness of the model.

Among these 7 aalncRNAs,LINC01116 has been most extensively studied.LINC01116 is associated with a large number of cancers,like non-small cell lung cancer,liver cancer,breast cancer,prostate cancer,and glioma[53–59],among which studies on non-small cell lung cancer have been reported most frequently.LINC01116 can promote non-small cell lung cancer cell proliferation,invasion and metastasis[58],and silencing LINC01116 can inhibit the development of lung adenocarcinoma through the AKT signaling pathway[59].In addition,LINC01116 can also lead to gefitinib and cisplatin resistance in lung adenocarcinoma[60,61].In addition,lncRNA LINC00942 can promote METTL14-mediated m6A methylation in the proliferation and progression of breast cancer cells[62].Additionally,lncRNA KLHL7-AS1,which has only been mentioned in Parkinson's disease,it can affect GPNMB expression and contributes to the development of Parkinson's disease[63].All these results indicates that the aalncRNAs screened in this study can have a great impact on the treatment and diagnosis of clinical diseases,particularly tumors,and further demonstrates the reliability of our results.

Furthermore,the GO functional enrichment analysis demonstrated that these aalncRNA were enriched in antigen processing and presentation of endogenous antigen,MHC protein complex and T cell receptor complex.The KEGG signaling pathway analysis demonstrate that these genes were in relation to the allograft rejection,asthma,autoimmune thyroid disease and intestinal immune network for IgA production.These results suggest that aalncRNAs in LUAD,which may also play a role in immunity,may have potential therapeutic value in the grand era of tumor immunotherapy.

Our study also has some limitations,firstly our RNA-seq was obtained from TCGA database and not validated from other databases,so we cannot be confident of the universality of the model we established.In addition,our study did not perform further experimental validation of these aalncRNAs and cannot yet reveal the more precise biological regulatory mechanisms of these aalncRNAs associated with autophagy in LUAD.We will continue to remedy the deficiencies of this study in future studies.

Nowadays,more and more lncRNAs are being developed as biological markers for lung cancer diagnosis and therapeutic targets,and it has been shown that the expression of UCA1 was significantly elevated in the plasma of lung cancer patients and was consistent with the expression in tumor tissues[64].Tang et al.used lncRNA microarray analysis to search for potential lung cancer biomarkers in blood and found that the expression of three lncRNAs(RP11-397D12.4,AC007403.1,ERICH1-AS1)was upregulated,and the combined positive and negative predictive values of the three lncRNAs were 0.72 and 0.87,respectively[65].Therefore,the aberrantly expressed lncRNAs could be used as biomarkers for the diagnosis of lung cancer.In addition,MALAT1 can promote tumor cell migration,invasion and tumor growth,and plays an important role in the development of lung cancer.A small molecule inhibitor,JMJD1A,was used to inhibit the expression of MALAT1 by binding to the promoter region of MALAT1 gene,thereby suppressing the migration and invasion ability of tumor cells[66].This demonstrates that the development of targeted drugs against lncRNAs in lung cancer could be a new direction in the treatment of lung cancer,and could even be of benefit to patients who are clinically untargeted or who are resistant or relapsed during treatment.Therefore,the aalncRNAs screened in this study also have the potential to become biological markers for the diagnosis and treatment of LUAD.

Conclusion

In conclusion,this study constructed a signature of 7 aalncRNAs to predict the prognostic value of LUAD patients.And more importantly,the functions and methods related to our signature may help to develop new diagnostic and therapeutic strategies for LUAD patients.So further biological experiments are needed to validate the results of our study.