The RNA-sequencing and clinical data of LUAD were downloaded from TCGA database (https://cancergenome.nih.gov/). Our study was based on the open resource data that were free for researching and publishing relevant articles with no ethical issues and other conflicts of interest. The process of this study is presented in Figure 1.

The mRNA expression data were derived from 535 LUAD patients and 59 healthy controls. Meanwhile, corresponding clinico-pathological data, including gender; age; pathologic T, M, and N stage; tumor clinical stage; and overall survival (OS) time were also obtained from TCGA database. Twenty-three out of 522 patients were excluded due to lack of information of OS (or the OS time was zero); therefore, lncRNA expression of 499 patients and their clinico-pathologic data were used for analysis. We presented the basic clinical information of these patients in Table 1.

We performed Gene Set Enrichment Analysis (GSEA) to research two datasets associated with hypoxia (HARRIS_HYPOXIA, WINTER_HYPOXIA_METAGENE) from the MsigDB database. Then, the genes in all samples were analyzed in the abovementioned gene sets. There were 239 hypoxia-related genes found from the statistically significant gene set. Finally, by Pearson’s correlation analysis, we identified hypoxia-related lncRNAs with the criteria of |correlation coefficient| > 0.3 and p-value < 0.001 and then constructed an mRNA–lncRNA coexpression network connected with hypoxia.

The Wilcoxon test was utilized to screen the differential expressed lncRNAs (DELs) between tumor and normal samples. The genes with p-value < 0.05 and |log₂ fold-change (FC)| > 1 were defined as DELs. Then, we utilized the R package “WGCNA” to construct a scale-free coexpression network for the all hypoxia-related lncRNAs by setting the soft threshold power value to 4. Finally, we selected two models highly correlated with cancer samples for followed analysis.

We first took the intersection of 601 DELs and 617 lncRNAs in the two modules (blue: 394, brown: 223). By using univariate Cox analysis, survival-related lncRNAs associated with hypoxia were identified. Then, LASSO regression analysis was used to further screen genes. At this time, we randomly divided the LUAD samples into two cohorts, training and validation cohorts. Finally, the lncRNAs were selected to construct a multivariate Cox regression model, and the risk score was calculated. Based on the median score of the training cohort, the patients were divided into high- and low-risk subgroups. Between the two groups, KM analysis was used to compare the survival time, and the ROC curve was used to evaluate the predictive power of the signature. In this way, the prognostic signature was constructed. In addition, in the other cohort, we performed the same procedure to evaluate the correctness of it.

We utilized GSEA v4.1 (http://www.gsea-msigdb.org/gsea/index.jsp) to perform GSEA between the low- and high-risk groups. After that, we carried out the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the differentially expressed mRNAs (DEGs) between the two groups, using “limma” and “clusterProfiler” packages.

We then downloaded 33 cancer types from the UCSC Xena database (https://xenabrowser.net/datapages/), and they are adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), esophageal carcinoma (ESCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), sarcoma (SARC), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma (UCS), cholangiocarcinoma (CHOL), stomach adenocarcinoma (STAD), and uveal melanoma (UVM). Meanwhile, six types of immune infiltration, namely, C1 (wound healing), C2 (INF-r dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant) were also downloaded from it. Spearman’s analysis was used to calculate the association of the lncRNAs with immune subtypes, tumor mutation burden (TMB), and stemness score (RNAss and DNAss).

All statistical analyses were accomplished with R software (version 4.0.5). The DEGs were identified by the Wilcoxon test. We used Pearson’s correlation analysis to calculate the correlation between lncRNAs and mRNAs associated with hypoxia. Meanwhile, we applied univariate and multivariate Cox regression analyses to evaluate the correlation between the risk score and clinical features. P-value < 0.05 was regarded as a significant outcome.

Through GSEA, we found that the WINTER_HYPOXIA_METAGEN gene set was significantly enriched [(FDR = 0.040), (Supplementary Table S1)], while the HARRIS_HYPOXIA gene set was not [(FDR = 0.164), (Figures 2A,B)]. Next, 1,629 hypoxia-related lncRNAs were screened according to the significant gene set (Supplementary Table S2).

We further studied hypoxia-related lncRNAs by the means of the Wilcoxon test and the weighted gene co-expression network (WGCNA) analysis, and we used the 601 differentially expressed hypoxia-related lncRNAs between normal and tumor samples to intersect with the 617 lncRNAs from the two modules of WGCNA analysis (Figures 2C–F).

After taking the intersection, we got 262 hypoxia-related lncRNAs (Supplementary Table S3). We then used these lncRNAs to construct the prognostic signature. First, 20 lncRNAs were screened by univariate Cox analysis [(p-value < 0.05), (Supplementary Table S4)]. Then, by using LASSO analysis, 20 variables were reduced to 11 potential predictors (Figures 3A,B). Finally, seven lncRNAs were identified by the multivariate Cox regression analysis in the training cohort. LncRNAs, AC022613.1, AC026355.1, GSEC, LINC00941, NKILA, HSPC324, and MYO16-AS1, were used to calculate the risk score (Figures 3C–E). Prognostic risk genes correlated with hypoxia were constructed (Table 2).

The unique risk score of patients was calculated through multivariate Cox analysis and the expression level. Risk score = expression of AC010980.2 × 0.490289217 + expression of AC026355.1 × −0.417505811 + expression of AL606489.1 × 0.293303854 + expression of ITGB1-DT × 0.289653582 + expression of AL034397.3 × −0.277395699 + expression of LINC01116 × 0.192509771 + expression of LINC01150 × −0.506323639. According to the median of the risk score, there were 123 patients in the high- and low-risk groups, respectively. Additionally, in the other cohort, the number was 130 and 114 patients using the same middle score.

By applying ROC curve analysis, in the training cohort, the area under the curve (AUC) was 0.740 at year 3 and 0.736 at year 5. KM analysis also showed that the model could be a valid prognostic indicator for patients (Figure 4A). Likewise, in the other cohort, the AUC values were 0.600 and 0.634, respectively (Figure 4B). From the result of KM analysis, we got the same trend with the training cohort (Figures 4C,D). Meanwhile, we found that in both cohorts, patients in the low-risk group had more survival time than those in the high-risk group (Figures 4E,F). In addition, in the training cohort, the result of the C-index was 0.715 and 0.643 in the other cohort.

Among these lncRNAs, two lncRNAs in the training cohort (AC026355.1 and HSPC324) were upregulated, while MYO16−AS1, GSEC, NKILA, AC022613.1 and LINC00941 were downregulated in the low-risk group. In addition, similar results were obtained in the validation cohort (Figures 5A,B). Moreover, in both cohorts, we used univariate and multivariate Cox regression analyses to assess whether the risk score could serve as an independent prognostic factor. The risk score was an independent factor revealed by the univariate Cox regression, and the HR of it was 1.442. In the multivariate analysis, the risk score also remained an independent prognostic indicator [p < 0.001, HR = 1.391, 95% CI: 1.258–1.540 (Figures 5C,D)] in the training cohort. In the validation cohort, the risk score was an independent prognostic indicator, too (Figures 5E,F).

In addition, to study the applicability of this model, we also conducted validations in different clinical subgroups. The patients were sorted by age (the 65 years or younger group and the more than 65 years group), gender (the male group and the female group), T stage of the tumor (the T1–T2 group and the T3–T4 group), M stage of the tumor (the M0 group and the M1 group), N stage of the tumor (the N0 group and the N1–N3 group), and tumor stage (the stage I–II group and the stage III–IV group). The results showed that the survival rate of high-risk patients with different age, gender, M0, N0, T1–T2, and stage I–II group was significantly different from that of low-risk patients (p-value < 0.05) in the training cohort (Figures 6A–H). In the validation cohort, we obtained the similar result (Supplementary Figure S1).

We then constructed two nomograms that integrated the risk score of the seven-lncRNA models and clinico-pathological features to predict survival probability of patients. Based on these, we predicted the patient’s 1-, 3-, 5-, and 10-year survival probabilities (Figures 6I,J). In both the nomograms, the higher the total points calculated, the worse the prognosis. Meanwhile, the calibration plot for the prediction of 3-year and 5-year survival also indicated the consistency between observation and prediction in both the cohorts (Figures 6K–N).

We used GSEA software to perform KEGG analysis for exploring which pathways were enriched. The results identified that in the high-risk group, processes such as cell cycle, DNA replication, and mismatch repair were enriched by using the training cohort (Figures 7A,B).

We further screened 660 different genes (adjust p-value < 0.05 and |log2FC (fold-change)| > 1) between the high- and low-risk groups and carried out enrichment analysis based on these in the training cohort. GO enrichment analysis indicated that many nuclear biological processes or molecular functions were significantly enriched (Figures 7C,D). Next, KEGG pathway analysis indicated that the biological processes related to cell proliferation were enriched, based on the upregulated genes, while downregulated genes were highly correlated with the immune process (Figures 7E,F and Supplementary Figures S2, S3). It was further verified that hypoxia is related to cell proliferation and immunity.

We calculated the proportion of 22 immune cells of all samples in the training and validation cohorts by using the CIBERSORT algorithm (Figures 8A,C). The violin plot showed that patients in the high-risk group had a higher proportion of activated memory CD4+ T cells, resting NK cells, M0 macrophages, and activated mast cells and a lower proportion of regulatory T cells, resting mast cells, and resting dendritic cells than those in the low-risk group in the training cohort (Figure 8B). In addition, in the validation cohort, CD8+ T cells, activated memory CD4+ T cells, resting NK cells, M0 macrophages, and activated mast cells were higher in the high-risk group than in the other group, whereas monocytes, M2 macrophages, resting mast cells, and resting dendritic cells were lower (Figure 8D). This result indicated that immune-related activities were associated with hypoxia.

In addition, based on the prognostic signature, there was an observably different distribution between high- and low-risk groups through the principal component analysis, which indicated that there was a difference in the hypoxia phenotype of the model (Figures 9A,B).

From the abovementioned results, the seven-model lncRNAs played an important role in LUAD. We then downloaded 33 cancer types to understand the function of the lncRNAs and selected AC022613.1, GSEC, LINC00941, and NKILA for further study. We found that AC022613.1, GSEC, LINC00941, and NKILA were mainly upregulated in tumor samples compared with normal samples (Figures 9D–G). In addition, the expression of these lncRNAs varied in different tumors (Figure 9C). The expressions of LINC00941 and NKILA were highly expressed in CHOL samples than those in the other tumors. In LUAD tissues, the expressions of the four lncRNAs were less than zero (Figure 9H). In addition, we found that LINC00941 and NKILA may have similar functions (Figure 9I).

By plotting KM curves for 32 cancer types, we found that only AC022613.1 significantly affected the survival of patients in many cancer types (Figures 10A–F). Furthermore, we used univariate Cox analysis to investigate the relationship between the expression of AC022613.1 and patient survival. The result showed that the relationship was different in different tumors (Figure 10G). We then explored their role in the six immune subtypes, stromal, ESTIMATE, tumor purity, immune score, tumor stem cells, and TMB. In LUAD patients, we found that AC022613.1 was strongly connected with the tumor stage and GSEC, LINC00941, and NKILA were significantly connected with immune subtypes (Figures 10H,I). Based on the ESTIMATE analysis, we researched the connection between the four lncRNAs and tumor microenvironment in LUAD. The results showed that GSEC had a negative association with stromal, ESTIMATE, and immune score, while LINC00941 had the opposite result (Figure 10J). Meanwhile, the results indicated that in pan-cancer, the four lncRNAs were strongly associated with immune subtypes (Figure 9J). Moreover, we found that they were mainly positively connected with stromal, ESTIMATE, and immune score, while there was a negative correlation between these and tumor purity, RNAss, and DNAss in pan-cancer (Figures 10K–N). Furthermore, we used the radar plots to distribute their association between the four lncRNAs and TMB. Distinctly, we found that NKILA and GSEC had strongly correlation in LUAD patients (Figures 10O–R).