The lncRNA and mRNA expression dataset in the FPKM format as well as the clinical characters for 535 LUAD patients and 59 normal patients were directly downloaded from the TCGA (https://portal.gdc.cancer.gov/), updated until 5 December 2021. GEO database (https://www.ncbi.nlm.nih.gov/geo/) was used to perform external validation.

DELs between the LUAD and normal samples were isolated from all lncRNAs using the R software. The p-value of each lncRNA in the LUAD and normal samples was calculated by the rank-sum test, and the p-values were rectified by the False Discovery Rate (FDR) method. Only the lncRNAs with adjusted p-values < 0.05 and log2 | fold change | values > 2 were defined as differentially expressed lncRNAs. Volcano plot and heatmaps were visualized by the “plot” function and the “heatmap” package of the R, respectively.

First, we removed LUAD patients with OS < 0 days. Next, we employed univariate Cox proportional hazards regression (CPHR) analysis and Kaplan-Meier analysis to assess the presence of any significant correlations between the expression of each DELs and the OS of LUAD patients. Only lncRNAs with p < 0.01 from both the analyses were considered with a logical agreement in expression and prognostic effect and selected as the candidate OS-related lncRNAs. Then, half of the patients were randomly assigned as the “primary dataset” after removing patients with incomplete clinical information; the original complete dataset was called the “entire dataset”. In addition to randomization, the criteria for grouping included no statistical differences in the clinical characteristics between the “primary” and “entire datasets. In order to fit the prediction model with the best-prediction effect, multivariate CPHR analysis (stepwise model) of candidate OS-associated lncRNAs was performed with the R software in the “primary dataset”. To ensure the goodness of fitting and to avoid overfitting, the Akaike information criterion (AIC) was computed, and the prediction model with the lowest AIC was considered as the most ideal. LncRNAs included in the best prediction model were selected as OS-related lncRNAs.

We determined the coefficients for each lncRNA by another multivariate CPHR analysis in the “primary dataset”. Until this point, we confirmed a risk score formula with the expressions of the OS-related lncRNAs as the independent variables and weighted by the regression coefficients corresponding to the lncRNAs. The risk scoring formula used is given below: R i s k   S c o r e = β 1 × E x p r e s s i o n   g e n e 1 + β 2 × E x p r e s s i o n   g e n e 2 + ⋯ + β n × E x p r e s s i o n   g e n e n where βi correspond to the correlation coefficient.

To determine whether the OS-related lncRNA signature was an independent predictor of OS, we applied both univariate and multivariate CPHR analyses of OS-related lncRNA signature and the routine clinical risk factors (such as sex, age, TNM stage, tumor stage, lymph node metastasis, and distant metastasis) in the LUAD patients. Next, we assessed whether the predictive effect of the OS-related lncRNA signature on OS was independent of the routine clinical risk factors by stratified analysis. Meanwhile, to evaluate the prognostic effect of the lncRNA-based classifiers across different time ranges, we plotted the time-dependent receiver operating characteristic (ROC) curves and then calculated the area under the time-dependent ROC curve (AUC) values for each dataset. Finally, the predictive effects of the 5-lncRNA classifier and the classifiers based on the other clinical risk factors were compared by AUC.

Whether the OS-related lncRNA was differentially expressed between the male and female patients was determined by the rank-sum test using p < 0.05 as the significance threshold. In both the male and female groups, the patients were assigned into two groups of high or low expression bounded by the median expression of an OS-related lncRNA, and the Kaplan–Meier curve was applied to analyze whether there were differences in survival time between the high and low expression groups. The lncRNA was considered to be with gender dimorphism if an OS-related lncRNA was differentially expressed in males and females while showing different prognostic association in males and females.

The co-expression degree of OS-related lncRNAs and mRNA was determined by Pearson’s correlational analysis. The mRNAs with a positive correlation coefficient >0.5 with OS-related lncRNAs were employed in the next step of enrichment analysis. The “cluster profile” package in R software was used for the functional enrichment analysis using the Gene Ontology (GO) terms and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathways (Kanehisa and Goto, 2000; Kanehisa, 2019; Kanehisa et al., 2021), with p < 0.01 set as a significance threshold.

For the survival analysis, the survival curves were plotted by the Kaplan-Meier method, and the differential p-values were calculated by the log-rank test. The t-test was used to compare the presence of any significant differences between the “primary” and “entire datasets”. Unless otherwise specified, p < 0.05 was considered to indicate a statistical difference.

The flow chart illustrated in Figure 1 shows the overall design of this study and some of the main results. After data collation, we obtained the expression data of 14,142 lncRNAs and 19,658 mRNAs for 535 LUAD samples and 59 normal samples from the TCGA-LUAD database. Through statistical comparison, 1,223 DELs in tumor samples and normal samples were identified with a log2 | fold change |> 2 and adjusted p < 0.05. Of these 1,223 DELs, 1,044 lncRNA were upregulated and 179 were downregulated in the LUAD patients. Next, volcano plots and heatmaps of the differential genes were drawn using the “plot” function and the “pheatmap” package in the R software, the results of which are illustrated in Figures 2A,B.

After the exclusion of 45 LUAD samples with incomplete survival data, 490 LUAD samples were finally enrolled in the study. In these 490 samples, 1,223 DELs were analyzed by the Kaplan-Meier method and univariate CPHR analysis, where OS served as the dependent variable and lncRNA expression as the independent variable. The results of the univariate CPHR analysis are depicted in Supplementary Table S1, and a total of 15 lncRNAs were found to be statistically significantly associated with OS in LUAD patients (all p < 0.01). Of this 15 lncRNA, the high expression of 13 lncRNAs (namely, LINC02081, AC010343.3, LINC02086, AC068228.1, AC022784.1, SATB2-AS1, AL138789.1, LINC01843, LINC00519, AL606489.1, DEPDC1-AS1, AC087588.2, and FAM83A-AS1) was associated with a shorter OS. In contrast, the high expression of AC026355.1 and AL031600.2 was associated with a higher OS. Moreover, as shown in Supplementary Figure S1, the results of the Kaplan-Meier analysis conformed to those of the univariate CPHR analyses. To this point, 15 lncRNAs with some correlation between the gene expression volume and prognosis were included as the candidate OS-related lncRNAs.

After removing 11 samples without complete clinical features (such as TNM stage or age), 479 LUAD samples formed the “entire dataset”, of which 239 groups were randomly selected as the “primary dataset”. The differential analysis revealed no statistical differences in the baseline clinical risk factors and OS between the “entire” and “primary datasets” (all p > 0.05; Table 1).

Candidate prognosis lncRNAs were further screened by multivariate-CPHR analysis (stepwise model) in the “primary dataset” using AIC to avoid overfitting. Five OS-related lncRNAs were picked with the largest fit and the lowest AIC values (Table 2), namely, is AC068228.1, SATB2-AS1, LINC01843, AC026355.1, and AL606489.1. Next, these 5 OS-related lncRNAs and their risk coefficients were integrated into the predictive signature to obtain a risk scoring using the following formula: R i s k   S c o r e = 0.4537 × expression   AC 068228.1 + 2.1929 × e x p r e s s i o n   SATB 2 – AS 1 + 0.0824 × e x p r e s s i o n   LINC 01843 + – 0.3963 × e x p r e s s i o n   AC 026355.1 + 0.1383 × e x p r e s s i o n   AL 606489.1

Next, we computed the risk score for LUAD patients in the “primary dataset” according to the 5 lncRNA signatures. Using the median risk score (0.09382005) as the cut-off value, 239 LUAD patients were classified into high- (n = 119) or low- (n = 120) risk groups. The risk score distributions, OS status, and the 5 lncRNA expression profiles in the “primary datasets” are depicted in Figure 3 (A–C). OS-related lncRNAs expression heatmaps revealed that the 4 upregulated lncRNA (i.e., AC068228.1, SATB2-AS1, LINC01843, and AL606489.1) demonstrated higher expression levels in the high-risk group, and the AC026355.1 expression levels were lower in the high-risk groups. As shown in Figure 3D, the Kaplan–Meier curve obviously showed that the OS time in the high-risk group was less than that in the low-risk group (p = 1.071E-04, log-rank test). Subsequently, in the “primary dataset”, as shown in Figures 3E–G, the AUC of the time-dependent ROC curve was 0.768 at 1 year, 0.668 at 3 years, and 0.702 at 5 years.

To verify the prediction of 5-lncRNA signatures obtained from the “primary dataset”, we applied 5-lncRNA signatures to the “entire dataset” (n = 479). Similarly, 479 patients were classified into the high-risk (n = 244) and low-risk (n = 235) groups according to the median risk score in the “primary dataset”. The risk score distributions, OS status, and the 5 lncRNA expression profiles in the “entire dataset” are depicted in Figures 4A–C. The results from the “entire dataset” are consistent with those from the “primary dataset”. Meanwhile, the Kaplan–Meier curve (Figure 4D) showed that the OS in the high-risk group (n = 244) was significantly shorter than that in the low-risk group (n = 235) (p = 5.587E-07, log-rank test). As shown in Figures 4E–G, the AUC of the time-dependent ROC curve was 0.738 at 1 year, 0.661 at 3 years, and 0.709 at 5 years. The 5-lncRNA signature showed a good prediction performance both in the “primary dataset” and the “entire dataset” of the LUAD patients. The prediction results of 5-lncRNA signature in “primary dataset” and the “entire dataset” were shown in Supplementary Table S2.

Next, to examine whether the prognostic performance of the 5-lncRNA features was independent of other conventional clinical risk factors, we performed multivariate CPHR analyses. The hazard ratio (HR) in the “entire dataset” (Table 3) was 1.085 (p < 0.001, 95% CI = 1.052–1.118), and in the “primary dataset” (Supplementary Table S3) was 1.065 (p < 0.001, 95% CI = 1.028–1.102). The abovementioned data indicates that these 5 lncRNA signatures could independently predict the prognosis of LUAD patients as an independent prognostic factor for LUAD.

To validate the scope of applicability of the risk score prediction, we conducted a stratified analysis of the “entire dataset”. First, considering the number of people, 479 LUAD patients were classified into stage I (n = 259; Figure 5A) and stage Ⅱ–Ⅳ (n = 220; Figure 5B) based on the TNM stage. Each subgroup was classified as the high-risk and low-risk groups and then Kaplan–Meier curves were accordingly plotted. Second, 479 patients were classified into no (n = 311, Figure 5C) or yes (n = 159, Figure 5D) subgroups according to the absence or presence of lymphoid tract metastasis, respectively. Next, 479 patients with LUAD were assigned into male (n = 219, Figure 5E) and female subgroups (n = 260, Figure 5F). Then, 479 patients with LUAD were assigned to the age ≥ 65 years (n = 266, Figure 5G) and <65 years subgroups (n = 213, Figure 5H). Finally, we noted that, in all subgroups, the survival time was significantly lower in the high-risk group than that in the low-risk groups, albeit it was not statistically significant in the female (p = 0.08) subgroup. To further validate the association between OS and 5-lncRNA, GEO database (GSE3141 and GSE19188) was used to perform external validation. The Kaplan-Meier curves for OS associated with the SATB2-AS1 expression were shown in Supplementary Figure S2 (GSE3141: p = 0.6312, GSE19188: p = 0.0914, GSE3141 + GSE19188: p = 0.1322). It is a pity that all three statistics did not show a significant effect. However, SATB2-AS1 still showed a clear trend towards promoting oncogenes, which is consistent with our findings in the TCGA database. All case ID involved in this study were shown in Supplementary Table S4.

We employed the time-dependent ROC curves to compare the predictive effects of different prognostic factors using the AUC as a comparison indicator. As shown in Figure 6, the stable predictive performance of the 5-lncRNA signature is more outstanding than the conventional clinical characters such as the TNM stage, and are efficient to predict the prognosis of LUAD patients.

Among the 5 OS-related genes, AL606489.1, SATB2-AS1 and AC068228.1 was differentially expressed between males and females (Figures 7A–C. This significant difference was not shown in LINC01843 (p = 0.5833) and AC026355.1 (p = 0.5177), as shown in Supplementary Figures S3A,B. The Kaplan–Meier curves for the OS related with AL606489.1 expression in males (low = 109, high = 110) and females (low = 130, high = 130) are depicted in Figures 7D,G respectively. For SATB2-AS1, the Kaplan–Meier curves in males and females are depicted in Figures 7E,H respectively. For AC068228.1, the Kaplan–Meier curves in males and females are depicted in Figures 7F,I respectively. In males, the high expression of AL606489.1, SATB2-AS1 or AC068228.1 was associated with the shorter OS. In females, the high expression of SATB2-AS1 or AC068228.1 was associated with the shorter OS. Dissimilarly, the high expression of AL606489.1 in females was not significantly associated with the OS (p = 0.2704). Finally, to verify whether this discrepancy was attributable to AL606489.1 association with the gender, we noted no significant difference in the overall survival between males and females by the Kaplan-Meier analysis (Supplementary Figure S3).

To determine the possible function of 5 OS-related lncRNAs in the tumorigenic development of LUAD tumors, we conducted an function enrichment analysis on mRNAs co-expressed with OS-associated lncRNAs in 490 LUAD samples. The levels of the 928 mRNA expressions were positively associated with the level of at least one OS-related lncRNA (co-expression coefficient >0.50). The GO analysis indicated that these co-expressed mRNAs were enriched in 52 GO terms (Supplementary Table S5). These GO terms were mainly enriched in regulating the mRNA metabolic processes, RNA splicing, and ubiquitin-specific protease activity (Figure 8A). Similar findings were obtained from the KEGG pathway enrichment analysis (Figure 8B), such as the ubiquitin-mediated proteolysis pathway. Therefore, the characteristics of 5-lncRNA mainly affected the gene expression within the nucleus and may be related to cell cycle regulation.