The mRNA expression dataset TCGA-LUAD (normal: 59, tumor: 535), as well as the corresponding clinical information data, were accessed from TCGA database (https://portal.gdc.cancer.gov/). GSE31210 (normal: 20, tumor: 226, platform GPL570) and GSE72094 (tumor: 446, platform GPL15048) were downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/). Samples with survival time less than 30 days in each data set were excluded. Finally, 460, 226 and 386 cancer samples in TCGA-LUAD, GSE31210 and GSE72094 were reserved, respectively.

1811 immune-related genes (Supplementary Table S1) and 200 hypoxia-related genes (Supplementary Table S2) were acquired from ImmPort database (https://www.immport.org/home) and Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/), wherein 652 genes were not identified in TCGA-LUAD, GSE31210 and GSE72094 data sets. Therefore, survival analysis was conducted after removing these 652 genes from the data set.

Based on survival information of LUAD patients, PACKAGE survival (Therneau and Grambsch, 2000) was applied to perform univariate Cox regression analysis on immune-hypoxia related genes in TCGA-LUAD and GSE72094 cohorts, and genes significantly correlated with patient overall survival (OS) were screened (p < 0.01). Intersection of survival-related genes in the two cohorts was taken, and overlapping genes were used for further analysis. Next, TCGA-LUAD was taken as the training set. By using R package glmnet, LASSO Cox regression analysis was done on genes that were remarkably associated with OS (Friedman et al., 2010). Appropriate punishment parameter lambda was selected by the ten-fold cross-validation selection method to remove genes with strong correlation to reduce model complexity. Next, by using R package survival, multivariate Cox regression analysis was conducted on genes screened by LASSO Cox analysis to obtain prognostic genes with immune-hypoxia feature in LUAD. The prognostic risk assessment model of LUAD was built according to risk coefficient and expression of feature genes. R i s k s c o r e = ∑ i = 1 n exp i ⁡ ∗ β i n is the number of selected genes, exp _i represents expression of gene i, and β _i indicates regression coefficient of gene i.

The riskscore of each cancer patient in TCGA-LUAD cohort was calculated, and patients were assigned into high- or low-risk groups by the median riskscore value. Survival analysis was conducted, and R package survminer was employed to generate survival curves of the two risk groups (https://cran.r-project.org/web/packages/survminer/index.html). The relevant survival state graph was plotted according to the riskcsore of the samples. Subsequently, ROC curves for prognostic models were plotted by R package survival ROC to predict patients’ OS at 1, 3 and 5 years (https://rdrr.io/cran/survivalROC/man/survivalROC.html). Independent data sets GSE31210 and GSE72094 were used as validation sets to verify accuracy of the model. Finally, effects of prognosis-related gene expression on prognoses of LUAD patients were analyzed on the GEPIA website (http://gepia.cancer-pku.cn/index.html).

Twenty-nine kinds of immune cells and immune-related functional gene sets were got from the published literature (Jiang et al., 2021). We assessed immune cell infiltration level in TME of LUAD patients based on these 29 immune data sets. ssGSEA was on samples in two groups using R package GSVA (Hänzelmann et al., 2013) to study immune cell infiltration level of patients with different riskscores.

GSVA estimates pathway activation variation in sample populations in an unsupervised manner, which is an excellent molecular mapping method (Hänzelmann et al., 2013). Enrichment scores for varying biological pathways in each sample in TCGA-LUAD cohort were calculated using R package GSVA. From MsigDB (http://software.broadinstitute.org/gsea/msigdb/index.jsp), the “c2.cp.kegg.v7.0.symbols.gmt” data were downloaded as reference. Through GSVA analysis, the differential pathways between groups were obtained. The differential enrichment pathways were further identified, with the standard of FDR < 0.01 and p < 0.01.

Univariate and multivariate Cox regression analyses were conducted using riskscore and clinical information (gender, age, T, N, stage) to verify whether riskscore could be an independent prognostic factor for LUAD. To predict 1-, 3-, and 5-year survival rates of LUAD patients, a nomogram was constructed based on riskscore and clinicopathological information using R package rms (Huang et al., 2019). Finally, calibration curves were drawn to evaluate consistency in predicted survival time and actual survival time of patients.

All statistical analyses were done on R software (3.3.1). The results were not statistically significant if data did not follow a normal distribution. Differences between groups were compared by Kruskal–Wallis test or Wilcoxon rank-sum test. Cox Regression model was utilized to perform univariate and multivariate analyses. Survival differences were assessed by logarithmic rank test. When p value < 0.01, the data were statistically significant.

1,328 immune-hypoxia-related genes that were overlapped in TCGA-LUAD and GSE72094 were identified (Figure 1A). Based on these genes as well as the survival information of TCGA-LUAD and GSE72094 cohorts, univariate Cox regression analysis was completed, and 137 and 204 genes associated with survival in the two cohorts were obtained, respectively. By intersecting, 43 genes related to patient survival were identified (Figure 1B). LASSO Cox regression evaluated the regression coefficients of 43 genes associated with survival, and a fitting curve was plotted (Figure 1C). Then, an appropriate punishment parameter lambda was selected according to 10-fold cross-validation. Results showed that when lambda was −4, the model contained 17 survival-related genes, which had the lowest complexity (Figure 1D). The above 17 genes underwent multivariate Cox regression analysis, and complexity and excellence of different models were evaluated according to Akaike Information Criterion (AIC). The model showing the lowest AIC value was chosen as the prognostic model, which contained seven immune-hypoxia-related genes significantly related to prognoses of LUAD patients. Among these genes, S100P, S100A16, PGK1 and TNFSF11 were risk factors and ARRB1, NCR3 and TSLP were protective factors (Figure 1E). As the risk coefficients and expression of seven genes mentioned above, the immune-hypoxia feature prognostic model of LUAD was riskscore = 0.065145908*coef (S100P) − 0.27104694*coef (ARRB1) + 0.104952032*coef (S100A16) + 0.152008632*coef (PGK1) + 0.088969471*coef (TNFSF11) − 0.174725763*coef (NCR3) − 0.056725904 *coef (TSLP). In conclusion, seven immune-hypoxia-related genes in LUAD were identified by LASSO and Cox repression analyses combined with patient survival information.

Following risk coefficient and expression of 7-gene signature, the riskscore of each LUAD patient in TCGA-LUAD cohort was computed, and riskscore distribution map was plotted. Patients were assigned into either high- or low-risk groups (Figure 2A). On the basis of survival time and riskscore of LUAD patients in TCGA-LUAD cohort, we drew a scatter plot of patients’ survival time, and the results displayed a negative correlation of survival time and riskscore (Figure 2B). Meanwhile, we plotted the expression heat maps of the seven feature genes in two risk groups. The heat maps showed that risk factors S100P, S100A16, PGK1 and TNFSF11 expressed highly in high-risk groups, while protective factors ARRB1, NCR3 and TSLP expressed highly in low-risk group (Figure 2C). Survival curve showed that prognoses of patients in low-risk group was notably better than in high-risk group (Figure 2D). Subsequently, ROC curves of 1-, 3- and 5-year survival predicted by the 7-gene signature were plotted. AUC values of 1-, 3- and 5- survival were 0.73, 0.696 and 0.649, respectively (Figure 2E), suggesting that this prognostic model could evaluate the prognoses of LUAD patients well. Overall, the 7-gene signature presented favorable predictive performance in the training set TCGA-LUAD cohort.

To verify predictive potential of the 7-gene signature, riskscore of patients in GSE72094 and GSE31210 data sets was calculated, and samples were assigned into high- or low-risk groups. Survival analysis results showed poor OS in high-risk group (Figure 3A). ROC curves of 1-, 3- and 5-year survival predicted by signature in GSE72094 and GSE31210 data sets were plotted. AUC values of 1-, 3- and 5-year survival in GSE72094 were 0.723, 0.685 and 0.711, respectively, while in GSE31210 were 0.774, 0.686 and 0.757 (Figure 3B). In addition, to explore the influence of seven genes on prognoses of LUAD patients, the survival curve was plotted. The results showed that S100P, TNFSF11, ARRB1 and NCR3 directly affected the prognoses of LUAD patients, while S100A16, PGK1 and TSLP did not (Figures 4A–G). In conclusion, this 7-gene signature based on immune-hypoxia related genes was able to accurately predict prognoses of LUAD patients.

To evaluate the differences in TME immune cell infiltration between high- and low-risk patients, ssGSEA was done. Results exhibited that the expression of immunoreactive gene sets, including APC_co_stimulation, APC_co_inhibition, CCR, Check-point, HLA, Cytolytic_activity, Parainflammation, Inflammation-promoting, T_cell_co-stimulation, T_cell_co-inhibition, Type_II_IFN_Response and Type_I_IFN_Response in low-risk group were prominently higher than those in high-risk group (Figure 5A), and the infiltration levels of immune cells B_cells, aDCs, DCs, iDCs, CD8+_T_cells, Mast_cells, Macrophages, Neutrophils, pDC, Tfh, T_helper_cells, TIL, Th1_cells and Treg were higher (Figure 5B). Based on the above findings, we noted that the immune score of patients in high-risk group was relatively low. Combined with the results of patient survival analysis, we concluded that patients in high-risk group may be immunosuppressed.

Based on enrichment score of each sample, we determined differences in enrichment pathways between two risk groups by GSVA (FDR < 0.01, p < 0.01). Heat map indicated that the pathways related to metabolism and mismatch repair, inluding KEGG_ONE_CARBON_POOL_BY_FOLATE, KEGG_GLYOXYLATE_AND _DICARBOXYLATE_METABOLISM, KEGG_PYRIMIDINE_METABOLISM, KEGG_HOMOLOGOUS_RECOMBINATION, KEGG_MISMATCH_REPAIR and KEGG_DNA_REPLICATION, were substantially activated in high-risk group. In low-risk group, immune-related pathways, including KEGG_INTESTINAL_IMMUNE _NETWORK_FOR_IGA_PRODUCTION, KEGG_PRIMARY_IMMUNODEFICIENCY, KEGG_ALLOGRAFT_REJECTION, KEGG_AUTOIMMUNE_THYROID_DISEASE and KEGG_GRAFT_VERSUS_ HOST_DISEASE, were significantly activated (Figure 6). Based on these results, patients in high-risk group may be immunosuppressed.

Finally, in TCGA-LUAD data set, univariate and multivariate Cox regression analyses were conducted combining clinical information (gender, age, T, N and Stage) with riskscore. The forest plots suggested that riskscore could be an independent factor of LUAD OS (Figures 7A,B). Then we integrated the clinical information to generate a nomogram to predict 1-, 3-, and 5-year OS of LUAD patients (Figure 7C). Finally, nomogram calibration curves presented that the predicted and actual OS values of LUAD patients were highly consistent (Figures 7D,E). In conclusion, the 7-gene signature could effectively predict prognoses of LUAD patients.