In this study, we incorporated NSCLC data from two publicly available databases. For the TCGA, the gene expression data along with the matching clinical data of lung adenocarcinoma (LUAD), as well as lung squamous cell carcinoma (LUSC) were obtained from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/). We combined LUAD samples and LUSC samples as training cohorts, which were called “TCGA” cohort, including 1129 samples.

The gene expression microarray of NSCLC (GSE37745 and GSE50081) with matching clinical data was abstracted from Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo). We integrated GSE37745 and GSE50081 samples as validation cohort, which were called “GEO” cohort, including 377 samples. Gene expression data of both the datasets were normalized using the Robust Multichip Average (RMA) approach from R package “affy”.

A total of 1793 immune-related genes were abstracted from immunology databases, as well as Analysis Portal (ImmPort) data resource (https://www.immport.org/home). A total of 103 ferroptosis-related genes were abstracted from the study by Luo H et al (Luo and Ma, 2021).

The protein–protein interaction (PPI) network of immune-related genes and ferroptosis-related genes was constructed and visualized using Cytoscape (https://cytoscape.org/). PPI data were obtained from the STRING (https://string-db.org/) database.

A univariate Cox proportional regression model was adopted to select OS-linked genes from both immune and ferroptosis gene sets in TCGA training data set. Overall, 42 prognosis-linked signatures were screened with p < 0.05, including 38 immune genes and five ferroptosis-related genes, among which NEDD4 was both the immune- and ferroptosis-related gene.

Next, the least absolute shrinkage and selection operator (LASSO) regression model was constructed to identify significant prognostic genes. A risk score was computed via considering the expression of optimized 32 signatures and correlation: Risk score = (exp gene1 * coef gene1) + (exp gene2 * coef gene2) + … + (exp gene32 * coef gene32). Patients with lung cancer were stratified into the high-risk group or low-risk group by the median of the risk score.

Principal component analysis (PCA) was used according to the expression profile of 32 prognosis-related signatures of the prognostic model in the training (TCGA cohort) and validation sets (GEO cohort). The log-rank test was adopted to assess the difference of the survival time between high-risk patients and low-risk patients. Kaplan–Meier plots were used to present the results.

A one-sided Wilcoxon rank sum test was adopted to explore the difference in the risk score between patients with various clinical characteristics, including sex, patient status, lymph node, tumor recurrence, and clinical pathological stage (TNM categorization of malignant tumors) in TCGA or GEO cohorts.

A chi-square test was implemented to evaluate the relationship of the clinical pathological stage group with the risk score group in TCGA dataset (Table 1).

A multivariable Cox proportional regression model was performed based on the risk score and clinical characteristics. Adjusted p < 0.05 signified statistical significance (Table 2).

We explored tumor immune invasion of TCGA and GEO cohorts using the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) approach by R software with the package “estimate” (Yoshihara et al., 2013). The ESTIMATE method assesses the number of stromal cells along with the invasion level of immune cells in samples. Single-sample GSEA (ssGSEA) in R package gsva was used to calculate the ferroptosis-related enrichment score which we called Ferroptosis Score for each sample based on ferroptosis-related gene sets.

The fraction of 22 tumor-invading immune cells was calculated based on CIBERSORT (https://cibersort.stanford.edu/index.php) (Newman et al., 2015) for TCGA and GEO cohorts. The one-sided Wilcoxon rank sum test was adopted to analyze the differences of infiltrative degree for immune cells, and p < 0.05 denoted statistical significance.

After that, we conducted expression assessment of five immune checkpoint–linked genes consisting of PDCD1 (code PD-1), BTLA and, CD274 (code PD-L1), along with CTLA-4 and CD47. The one-sided Wilcoxon rank sum test was carried out for exploring the differences of the expression of five immune checkpoint–related genes in the high-risk group and low-risk group in TCGA and GEO cohorts, and p < 0.05 denoted statistical significance.

A total of 23 cancer therapeutic-predicted signature sets that we used were obtained from several studies, including “Basal_differentiation”, “EMT_differentiation”, “Immune_differentiation”, “Mismatch_repair”, “Nucleotide_excision_repair”, “p53_signaling_cascade”, “Oocyte_meiosis”, “Proteasome”, “Spliceosome”, “Pyrimidine_metabolism”, “DNA_replication” “Systemic_lupus_erythematosus”, “EGFR_ligands”, “Viral_carcinogenesis”, “FGFR3-coexpressed_genes”, “PPARG_network”, “IDH1”, “KDM6B”, “WNT-β-catenin_network”, “VEGFA”, “Hypoxia”, “Cell_cycle”, and “Progesterone-mediated_oocyte_maturation” (Motz et al., 2014; Peng et al., 2015; Sweis et al., 2016; Mariathasan et al., 2018; Spranger and Gajewski, 2018; Seiler et al., 2019; Kamoun et al., 2020; Necchi et al., 2020). ssGSEA was adopted to calculate the enrichment score of the abovementioned therapeutic signature gene sets. The one-sided Wilcoxon rank sum test was adopted to explore the differences of the enrichment score between high- and low-risk groups.

IMvigor210 was used to predict immunotherapy response (http://research-pub.gene.com/IMvigor210CoreBiologies). It is a study to investigate the anti–PD-L1 antibody atezolizumab in patients with metastatic urothelial cancer (mUCC) (Mariathasan et al., 2018). We evaluated the difference of the risk score between responsive groups [PD (progressive disease), SD (stable disease), PR (partial response), and CR (complete response)].

The transcriptome data in TCGA cohort were used to create a comprehensive indicator from immune- and ferroptosis-related profiling (Figures 1A,B). We constructed PPI networks of immune genes and ferroptosis-related genes. Most immune genes directly interacted with the ferroptosis-related genes according to the STRING database (Figure 1C). After data processing, 1294 immune genes and 94 ferroptosis-related genes were used for subsequent model construction (Figure 1D).

A univariate Cox proportional regression model was used to explore the prognostic value of both immune- and ferroptosis-linked genes. Screened with p < 0.05, 42 prognosis-related genes were obtained (Figure 2A), including 38 immune-related genes and five ferroptosis-related genes, among which, NEDD4 was both an immune- and ferroptosis-related gene (Figure 2A). NEDD4 is an oncogene, which encodes E3 ubiquitin ligase. Next, a LASSO analysis was used to construct a prognostic model with 32 signatures: ACSL3, ACTG1, ANGPTL4, APOD, CD1E, CRHR2, CTF1, DEFB103B, DKK1, EREG, FGA, FGF4, HLA-DOB, IL2, INSL4, ITGA6, LCN1, NEDD4, PDGFB, PF4V1, PTX3, RFXAP, SEMA3C, SEMA7A, SHC1, SLC11A2, STC2, TNFRSF6B, UMODL1, VDAC1, VEGFC, and XCR1 (Figures 2B,C). Then, based on the expression of the optimized 32 signatures and correlation in TCGA cohort, we established a predictive model.

Every patient’ risk score in TCGA and GEO data sets was computed. For analyzing the accuracy of the signatures used for constructing the module, we visualized the expression of 32 genes and found that a majority of 32 genes were differentially expressed between the high-risk and low-risk groups (Figure 3A). Then, PCA was performed to investigate whether lung cancer patients could be distinguished according to the expression of the 32 signatures in TCGA and GEO data sets (Figures 3B,C).

Next, patients were categorized into the high- and low-risk group using the median risk score as the cutoff value. In TCGA data set, patients in the high-risk group had a remarkably worse OS (Figure 3D; p < 0.0001; log-rank test), and the number of alive patients in the low-risk group were more relative to those in the high-risk group (Figures 3E–G). In the independent validation set (GEO dataset), patients in the high-risk group also exhibited a remarkably worse OS (Figure 3H; p = 0.00042 log-rank test), and the alive patients in low-risk group were more than those in the high-risk group (Figures 3I–K).

We next explored the capacity of the prognostic model in clinical pathological characteristics. We first investigated the difference of clinical characteristics between the high-risk and low-risk groups (Table 1). Then, we evaluated the differences of the risk score between patients with different clinical characteristics (dead vs. alive), sex groups (female vs. male), tumor size groups (T1 vs. T2 vs. T3 vs. T4), lymph node (N0 vs. N1-N3), and pathological stage (Stage I vs. Stage II-Stage IV) in TCGA cohort. The risk score in dead patients was remarkably higher in contrast with alive patients (Figure 4A; p < 2.22e-16), and male patients were remarkably higher than female patients (Figure 4B; p = 0.00015). The risk score in stage T2, T3, and T4 was remarkably higher relative to stage T1, and the risk score in stage T3 was remarkably higher than stage T2 (Figure 4C; p < 0.05). Besides, the risk score in stage N1–N3 was remarkably higher in contrast with stage N0 (Figure 4D; p = 0.003), and stage II–stage IV was remarkably higher than stage I (Figure 4E; p = 2.7e-06). Moreover, the risk score of dead patients was remarkably higher relative to alive patients (Figure 4F; p = 0.00094), and male patients were remarkably higher than female patients (Figure 4G; p = 0.012). Stage II–Stage IV was remarkably higher than stage I (Figure 4H; p = 0.022), and recurrent patients were remarkably higher than non-recurrence patients (Figure 4I; p = 0.022) in the GEO cohort.

Further assessment was conducted to explore whether the risk score reveals prognosis in diverse subgroups of clinical features. In the T subgroups (T1, T2, and T3), N subgroups (N0, N1, and N2), M0 subgroup, and pathological stage (stage I, stage II, and stage III) of TCGA cohort, patients in the high-risk group exhibited a poor OS (Figures 5A–J; p < 0.05; log-rank test). In the pathological stage (stage I and stage II) of the GEO cohort, remarkably poorer OS was found in patients in the high-risk group (Figures 5K,L; p < 0.05; log-rank test).

A multivariable Cox proportional regression model was built in TCGA cohort using the risk score and clinical pathological stage groups to verify the prognostic potential and independence of the prognostic model from other clinico-pathologic characteristics. The result suggested that our prognostic model has a potential in clinical application (Table 2).

The discrepancy of immune infiltration and ferroptosis between different risk groups.

To investigate the relationship between the risk score and immune infiltration as well as ferroptosis, we analyzed the distribution of the ESTIMATE score (consists of Immune Score and Stromal Score) and Ferroptosis Score (enrichment score of ferroptosis-related genes) in each sample. Immune Score, Stromal Score, and Ferroptosis Score tended to increase with the escalation of the risk score in both TCGA and GEO cohorts (Figures 6A,B). To further confirm this trend, Pearson’s correlation analysis was calculated between the risk score and Immune Score, Stromal Score, and Ferroptosis Score, respectively. The results showed that they were positively related with the risk score in both TCGA and GEO cohorts (Figures 6C–H; p < 0.05; Pearson’s correlation analysis).

Next, we investigated the difference in the expression of five immune checkpoint–linked genes between high- and low-risk groups. The expression of PD-L1 in the high-risk group was remarkably greater than that in the low-risk group in TCGA and GEO cohorts (Figures 7A,B; p = 0.021, p = 0.0064). Besides, CTLA-4 and PD-1 expressions in the high-risk group were greater than those in the low-risk group in the GEO cohort (Figure 7B; p = 0.0198, p = 0.0297).

To explore immune cell infiltration of tumor samples, we calculated 22 immune cell abundances among TCGA and GEO cohorts by CIBERSORT. Next, we explored the difference of 22 immune cell invasion between high- and low-risk groups. The invasion of M0 Macrophages in the high-risk group was remarkably greater in contrast with that in the low-risk group in TCGA cohort (Figure 7C; p < 0.0001). The infiltration of resting NK cells in the high-risk group was remarkably higher than that in the low risk group in TCGA and GEO cohorts (Figures 7C,D; p < 0.01). Also, significantly greater infiltration of CD8 T cells was observed in the high-risk group relative to that in the low-risk group in TCGA and GEO data sets (Figures 7C,D; p < 0.01).

We next investigated the potential role of the prognostic model in the prediction of response to immunotherapy using the IMvigor210 cohort. We found that the risk score in non-responsive patients [stable disease (SD) and progressive disease (PD)] was significantly higher than responsive patients [complete response (CR) and partial response (PR)] (Supplementary Figure S1A; p < 0.05; one-sided Wilcoxon rank sum test). The number of non-responsive patients in the high-risk group was more than that in the low-risk group (Supplementary Figure S1B). The distribution of CR, PR, SD, and PD patients between the high-risk and low-risk groups was significant (Supplementary Figure S1C; p < 0.05; chi-square test). Besides, patients in the high-risk group had significantly poorer OS (Supplementary Figure S1D; p < 0.01).

For investigating the guiding role of the risk score in cancer treatment, we calculated the enrichment score for each sample according to 23 cancer therapeutic–predicted signature sets by ssGSEA. Next, we analyzed the differences of these scores in the high-risk and low-risk groups in TCGA and GEO cohorts. In TCGA cohort, the enrichment score of the high-risk group was higher relative to the low-risk group in 73.91% (17/23) of cancer therapeutic prediction signature sets, including “Basal_differentiation”, “Cell_cycle”, “DNA_replication”, “EGFR_ligands”, “EMT_differentiation”, “Hypoxia”, “Immune_differentiation”, “Mismatch_repair”, “Nucleotide_excision_repair”, “Oocyte_meiosis”, “p53_signaling_casacde”, “Progesterone-mediated_oocyte_maturation”, “Proteasome”, “Spliceosome”, “Pyrimidine_metabolism”, “Systemic_lupus_erythematosus”, and “Viral_carcinogenesis” (Figure 7E; p < 0.05). In the GEO cohort, the elevated score of “Basal_differentiation”, “EGFR_ligands”, “EMT_differentiation”, “Hypoxia”, “Immune_differentiation”, and “Proteasome” was observed in the high-risk group in contrast with those in the low-risk group (Figure 7F; p < 0.05).