This article is mentioned in:

Introduction

Lung cancer is the most common type of malignant tumor and was the leading cause of cancer-associated death in the USA in 2019 (1). Lung cancer can be categorized into two main histological subtypes, non-small-cell lung cancer (NSCLC) and small-cell lung cancer, accounting for ~85% and ~15% of cases, respectively (2,3). Lung adenocarcinoma (LUAD) is one of the major subtypes of NSCLC with a low 5-year survival rate of ~20% in the USA (4,5). It is currently understood that two of the driver oncogenes involved in tumorigenesis and progression of LUAD are epidermal growth factor receptor and anaplastic lymphoma kinase (6–8). These biomarkers can be used for diagnosis and targeted in the treatment of LUAD; however, prognosis remains unsatisfactory due to the underlying molecular heterogeneity and diverse etiology of LUAD. Therefore, it is important that novel valuable biomarkers, which have associations with prognosis in LUAD are identified.

Long non-coding RNA (lncRNA) is a type of non-coding RNA which is >200 nucleotides in length (9). Previous research has demonstrated that lncRNAs serve diverse functions in cellular processes and are involved in tumorigenesis, progression and metastasis (10,11). However, only a few lncRNAs involved in the development of LUAD have been identified, including histocompatibility leukocyte antigen complex P5, chromatin-associated RNA 10 and metastasis-associated lung adenocarcinoma transcript 1 (12–14). The underlying molecular mechanisms for the functions of numerous lncRNAs remains unclear, to the best of our knowledge, but it is hypothesized that lncRNAs may have valuable clinical applications in the future.

Previous studies have described several novel lncRNAs with potential use as biomarkers differentially expressed in human LUAD tissue (15,16) and further investigation into these lncRNAs may elude their prognostic value. Integrated prognostic analysis of lncRNAs combined with overall outcomes, pathological stages and other clinical parameters may improve the prediction of overall survival time of patients with LUAD.

In the present study, the construction method for a previous lncRNA gene signature was amended (17,18), aiming to improve its accuracy. Based on lncRNA expression profiling, a prognostic model was constructed using a comprehensive approach involving the least absolute shrinkage and selection operator (LASSO) algorithm and multivariate Cox regression. Model robustness and accuracy were assessed using area under the calculated curve (AUC) and Kaplan-Meier (K-M) survival analyses in high- and low-risk groups. Finally, potential biomarkers associated with survival status were identified using K-M survival analysis and possible functions of these biomarkers were predicted using enrichment analysis of co-expressed mRNAs.

Materials and methods

Datasets

RNA sequencing data and corresponding clinical information were obtained from The Cancer Genome Atlas (TCGA) database (cancer.gov/tcga), including data for 535 LUAD tissue samples and 59 control adjacent normal tissues. The LUAD group consisted of 249 males and 286 females with an age range of 33–88 years and a median age of 66 years, whereas the control group consisted of 25 males and 34 females with an age range of 42–86 years and a median age of 66 years. LUAD patients with incomplete survival data were excluded, therefore 500 samples were used for analysis.

Identification of differentially expressed lncRNAs

All analyses were performed using R software (version 3.5.3; r-project.org/) and R packages by R studio version 1.1.463 (19,20). Raw lncRNA expression profiles were constructed using an expression matrix and normalized using the edgeR package (version 3.22.5) (21) and the lncRNAs differentially expressed between LUAD and control samples were identified. lncRNAs which met the threshold of |log2[fold change (FC)]|≥2 and false discovery rate (FDR) <0.05 were screened out for subsequent analysis.

Definition of the gene-related prognostic model

Univariate Cox, LASSO and multivariate Cox regression analyses were performed to identify the association between survival time and lncRNA expression using the glmnet (version 2.0.18) (22) and survival (version 2.44.1.1) (23) packages. A univariate Cox model was established to assess the correlation between overall survival (OS) and lncRNA expression levels. P<0.05 was considered to indicate a statistically significant difference. lncRNAs meeting this criterion were included in LASSO regression analysis, which was to further select prognostic genes and to avoid the over-fitting of the signature model. After selection using the minimum λ value, which represented the optimal number of variables in this model, a further refined group of lncRNAs was included in a multivariate Cox analysis to determine the independent contribution of each lncRNA to prognosis and hazard ratios (HR) and 95% confidence intervals (CI) were calculated. The Cox regression coefficient (β) and expression levels of lncRNAs were used to calculate prognosis-related risk scores according to the formula: Σ (explncRNAn × βlncRNAn). Samples were further divided into high- and low-risk groups using the median risk score.

Evaluation of the risk model and identification of potential biomarkers

Time-dependent receiver operator characteristic (ROC) curves were plotted to evaluate performance of the prognostic model. K-M survival analyses were performed to predict survival times for the high and low risk groups and to calculate AUCs value for the prognostic model. Multivariate Cox analysis results were used to identify potential biomarkers. A threshold of P<0.05 was used to identify eligible lncRNAs for survival and K-M curve analyses and to predict the OS associated with each lncRNA. Subsequently, potential biomarkers were identified if this P<0.05 criterion from the survival curves was met.

Functional enrichment analysis

The biological function of the lncRNAs with significant P-values for OS in the comparison between the high and low risk groups was predicted using their co-expressed protein-coding mRNAs. The mRNAs that correlated with prognostic lncRNAs were identified using Pearson correlation analysis. An absolute correlation coefficient value >0.4 was considered significant. The Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Tool (version 6.8) (24,25) was used to perform the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis.

Results

lncRNAs are differentially expressed between LUAD and control samples

A total of 1,684 differentially expressed lncRNAs (Table SI) were identified, meeting the conditions of |log2 (FC)|≥2 and FDR <0.05. Of these, 1,499 lncRNAs were upregulated and 185 lncRNAs were downregulated (Fig. 1). The lncRNAs which had significant differential expression levels between the LUAD and control tissue groups were used for subsequent prognostic analysis.

Figure 1. Identification of 1,684 differentially expressed lncRNAs. Red dots represent 1,499 upregulated lncRNAs, green dots represent 185 downregulated lncRNAs and black dots represent 9,736 non-differentially expressed lncRNAs. lncRNA, long non-coding RNA; FC, fold change; sig, significance; down, downregulated; not, not significant; up, upregulated; FDR, false discovery rate.

Construction of the prognostic model

Univariate Cox regression was used to analyze the correlations between differentially expressed lncRNA profiles and the OS of patients with LUAD. A total of 107 lncRNAs were significantly correlated with patient OS (P<0.05) and these lncRNAs were used in a LASSO regression analysis (Table SII). Overall, 44 key lncRNAs were selected when the λ value was at a minimum (Fig. 2A and B). Finally, a multivariate Cox regression analysis was performed to obtain HR values and 95% CIs for each key lncRNA (Fig. 3). Samples were classified into high- and low-risk groups according to respective median risk scores (Fig. 4A and B).

Figure 2. Construction of the prognosis associated molecular signature by LASSO regression. (A) LASSO coefficient profiles of 107 prognostic lncRNAs. (B) Tuning parameter selection in the LASSO regression, with vertical lines marking the optimal λ value of minimum and 1-SE. LASSO, Least Absolute Shrinkage and Selection Operator; SE, standard error; lncRNA, long non-coding RNA.

Figure 3. Forest plot demonstrating the HR and 95% confidence interval for the key lncRNAs. HRs <1 represents a negative correlation and HR >1 represents a positive correlation. *P<0.05; **P<0.01; ***P<0.001. HR, hazard ratio; lncRNA, long non-coding RNA.

Figure 4. Patient risk scores. (A) Distribution of patient risk scores. Samples were divided into two groups according to the median risk score. (B) Survival status of patients with lung adenocarcinoma. The position of the dots represents the correlation between the survival time and risk scores.

Evaluation of the prognostic model

K-M survival curves for OS outcomes were performed according to median risk score values, demonstrating that the predicted survival time of the high-risk group was significantly shorter compared with the low risk group (P< 0.001; Fig. 5). The AUCs for 3- and 5-year survivals were 0.836 and 0.818, respectively, of a time-dependent ROC curve demonstrating that the risk score model had stable performance (Fig. 6).

Figure 5. Kaplan-Meier survival curves for overall survival outcomes based on the lncRNA signature. There was a significant difference between the high and low risk groups. lncRNA, long non-coding RNA.

Figure 6. Time-dependent receiver operator characteristic curves for overall survival prediction at 3- and 5-years. The signature had a moderately stable performance for prognosis prediction. AUC, area under curve.

Identification of potential biomarkers

Multivariate Cox regression showed that lncRNAs with P<0.05 were correlated with survival status (Fig. 3). K-M survival curves enabled screening out of two lncRNAs, AC124804.1 and MIR34AHG, with significantly different predicted survival times between the high- and low-risk score groups (Fig. 7A and B). These two lncRNAs had the greatest correlation with OS time and may have potential as prognostic biomarkers for patients with LUAD.

Figure 7. Kaplan-Meier survival curves for overall survival outcomes. Predicted survival time of: (A) AC124804.1 and (B) MIR34AHG long non-coding RNAs.

Functional enrichment

To investigate the underlying molecular mechanisms of the functions of AC124804.1 and MIR34AHG, the co-expressed mRNAs of these lncRNAs were identified using Pearson's correlation analysis. For MIR34AHG, no mRNAs met the screening condition of a correlation coefficient value of >0.4. For AC124804.1, a total of 818 protein coding mRNAs were identified for further functional enrichment. Overall, four GO terms and four KEGG pathways had confirmed association with AC124804.1. The 818 co-expressed mRNAs were primarily clustered in the poly(A) RNA binding and protein binding of molecular functions, and to the mitochondrion and nucleoplasm cellular components (P<0.05; Fig. 8A). KEGG pathway enrichment analysis demonstrated that these mRNAs were primarily associated with Proteasome, Alzheimer's disease, SNARE interactions in vesicular transport and Peroxisome pathways (Fig. 8B).

Figure 8. Functional enrichment analysis for the co-expressed protein-coding mRNAs of AC124804.1. (A) Significantly enriched GO terms. The bar lengths indicate the counts enriched in this category. (B) Significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways. The bar lengths indicate the counts enriched in these pathways. GO, Gene Ontology; MF, molecular function; CC, cellular component.

Discussion

LUAD is a large global health challenge (1,26) and our knowledge of its underlying mechanisms may be improved by preclinical next-generation sequencing bioinformatics technology (27,28). However, despite advances in the diagnosis and treatment of LUAD, the OS time of patients with LUAD remains poor (29). Recent studies have identified molecular biomarkers throughout the various clinical stages and pathological tissue types (30,31). However, patients with the same clinical and pathological features at diagnosis often have a completely different final prognoses (32). Molecular gene signatures associated with prognosis may aid the identification of novel molecular biomarkers for use in diagnosis and treatment of LUAD in combination with patients' genetic profiles and clinical parameters (17,18). Recent studies pertaining to LUAD have investigated gene signatures based on lncRNA data and the selected lncRNAs may be important risk factors and have potential predictive power for patient survival (17,18). However, due to complex underlying mechanisms, the biological functions of the majority of lncRNAs have not been resolved (33). The prognostic lncRNA signature constructed in the present study, demonstrating correlation with clinical features of LUAD, may be aid the screening out of potential lncRNAs for further investigation.

To improve the accuracy and robustness of the present model compared with previous studies (15,17,18), a different approach was taken to data preprocessing and construction of the model. First, the present study utilized larger sample sizes compared with previous studies (15,17,18) resulting from differences in inclusion criteria, including a total of 535 LUAD tissues and 59 adjacent normal tissues (controls). Second, the present study used different methods of data preprocessing and normalization. For example, previous studies using some online tools for analysis, whereas the present study performed analysis using edgeR package in R software, which could analyze the latest raw counts of RNA-seq data of TCGA. Third, the present group selected variables using univariate Cox analysis and LASSO regression, which were not employed in previous studies. In addition, the eligible lncRNAs in the present study were confirmed using P-values of multivariate Cox regression and survival analyses.

In the present study, lncRNA profiles were used to construct the prognostic model, using data from tissue samples and clinical information of 500 patients held in the TCGA database. A total of 1,684 differentially expressed lncRNAs were investigated as potential biomarkers using univariate Cox regression to screen out lncRNAs that were significantly associated with the clinical characteristics. The results of univariate Cox regression were further filtered using the LASSO algorithm, which weakens collinearity of the risk model and enables the most influential variables to be selected (34,35). This approach is suitable for dealing with datasets containing a large volume of variables, as occurs in genomics (36). The final 44 lncRNAs identified using LASSO regression were then analyzed using multivariate Cox regression. The correlation of each lncRNA with survival of patients with LUAD was presented as an HR value. Patients were divided into high and low risk groups according to the median risk scores which demonstrated that 3 and 5-year survival times were significantly different between the high and low risk score groups. These data combined with the AUC score results confirmed the accuracy and robustness of this model for predicting prognosis of patient with LUAD.

A multivariate Cox regression forest plot identified lncRNAs with significant P-values as potential prognostic biomarkers for LUAD. The respective survival times correlated with these lncRNAs were calculated prospectively using K-M survival analysis, identifying candidate prognostic factors. According to the prognostic signature constructed in the present study, the expression levels of AC124804.1 and MIR34AHG were significantly correlated with differential clinical outcomes and predicted OS times of the two risk groups. These two lncRNAs may have value as novel biomarkers for estimating prognosis of patients with LUAD. To the best of our knowledge, there is no published work describing the relationships between these lncRNAs and LUAD to date. To explore the potential functions of AC124804.1 and MIR34AHG, the co-expressed protein-coding mRNAs of these lncRNAs were investigated using a Pearson correlation coefficient test followed by GO and KEGG functional enrichment analyses. These analyses demonstrated that the potential underlying molecular mechanisms of AC124804.1 function were pathways involving proteasome, peroxisome and SNARE interactions in vesicular transport, which have confirmed roles in tumorigenesis and progression in previous studies (37–39). MIR34AHG is the host gene of miR-34a and consequently may have biological functions associated with the miR-34 family, which may have suppressive functions in LUAD as previously demonstrated (40,41). The relationship between these two lncRNAs and LUAD should be studied further.

There are several limitations to the present study. First, the data for the lncRNA signatures were from the TCGA database and these samples primarily originated from Caucasians and African Americans, limiting extrapolation of the results of the present study to other ethnic groups. Second, the present study used a database mining design without validation in fresh samples and prospective experimental studies. Third, the clinical parameters included only OS time, with histological subtypes and other risk factors for LUAD were excluded from the analysis. The primary purpose of the present study was to identify novel biomarkers for the diagnosis and treatment of LUAD. The lncRNA prognosis signature should be further investigated incorporating more specific clinical characteristics to fully understand the associations involved.

In conclusion, a robust lncRNA signature that could stratify and predict survival time of patients with LUAD was constructed. Furthermore, Cox regression was used to identify two novel survival associated lncRNA biomarkers and provided potential targets for diagnostic and therapeutic application. The biological functions of the two lncRNAs identified in the present study need confirmation in further experimental studies.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was funded by The National Natural Science Foundation of China (grant no. 81603577).

Availability of data and materials

The datasets used and/or analyzed during the present study available from the corresponding author on reasonable request.

Authors' contributions

YZ conceived and designed the study. XY conducted the analysis and wrote the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References