Public lung adenocarcinoma transcriptome datasets were obtained from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). Patients without prognostic information were excluded from the study. A total of six datasets (TCGA-LUAD, GSE13213, GSE30219, GSE31210, GSE3141, and GSE41271) were gathered for further analysis (Bild et al., 2006; Tomida et al., 2009; Okayama et al., 2012; Rousseaux et al., 2013; Sato et al., 2013). For RNA-Seq data from the TCGA, FPKM values were downloaded from the TCGA Data Portal¹ and transformed into TPM values before analysis. For datasets in the GEO, normalized probe expression values generated from the microarray were acquired using the R package GEOquery (Davis and Meltzer, 2007). The probes in each dataset were annotated according to the appropriate platform annotation file. When multiple probes mapped to the same gene, the average expression values of these genes were used as the final expression value. To calculate tumor mutation burden (TMB) in TCGA-LUAD, gene-level mutation (VarScan) data were extracted from the TCGA Data Portal. In the present study, TMB was defined as the total number of non-synonymous mutations divided by the size of the exome.

A comprehensive list of ferroptosis drivers and suppressors was acquired from the FerrDb database² (Zhou and Bao, 2020). After removing duplicated entries and performing univariate Cox regression for each gene, 37 unique prognosis-associated genes were obtained for unsupervised clustering. The identification of ferroptosis subtypes was conducted with the TCGA cohort by applying the partitioning around medoids cluster algorithm in the R package ConsensusClusterPlus (Wilkerson and Hayes, 2010). The number of repetitions was set to 100 to ensure the stability of the result. The other parameters of the cluster algorithm were set to the default values.

To determine the biological function associated with different ferroptosis subtypes, GSVA was performed to detect subtle pathway activity in each sample using the R package GSVA (Hänzelmann et al., 2013). GSVA is a non-parametric, unsupervised method for calculating gene set enrichment through expression profiles. The gene set file named “c2.cp.kegg.v7.4.symbols.gmt,” including 186 canonical pathways, was downloaded from MSigDB and was chosen for pathway enrichment (Liberzon et al., 2011). A pathway comparison was conducted using the R package limma (Smyth, 2005), and an adjusted p value of less than 0.05 was considered statistically significant.

Two widely validated algorithms, CIBERSORT and ImmuCellAI (Newman et al., 2015; Miao et al., 2020), were implemented to estimate the abundance of specific immune components. CIBERSORT is a computational approach reinforced by support vector regression to quantify 22 immune cell fractions. ImmuCellAI was used for abundance estimation of 24 immune cell types, including 18 T cell subsets, from gene expression matrices. The absolute fraction and abundance of immune cells were estimated by CIBERSORT and ImmuCellAI, respectively. Only intergroup comparisons were conducted and the differences in immune cell abundance were depicted between risk groups and among subtypes.

To identify common DEGs among the three ferroptosis subtypes, gene expression was quantified by TPM value and empirical Bayes statistics was applied for differential expression analysis using limma package (Smyth, 2005). The obtained p-values were further adjusted by the Benjamini–Hochberg method. The DEGs were defined as genes with an adjusted p-value less than 0.001. By taking the intersection of the outputs, common DEGs of the three subtypes were obtained and illustrated using a Venn diagram. Next, the Boruta feature selection algorithm was applied to narrow the genes for model construction using the R package Boruta (Kursa and Rudnicki, 2010). Principal component analysis (PCA) was p to verify the ability to distinguish subtypes using selected genes. Furthermore, genes related to prognosis were screened out through univariate Cox regression and served as input for model training. To generate a robust ferroptosis-related prognostic (FRP) signature, LASSO Cox regression with 10-fold cross-validation was performed using a random seed, and the whole procedure was repeated 1,000 times (Tibshirani, 1997). The model with the highest frequency was selected as the final signature. Patients were subsequently stratified into high-risk and low-risk groups according to the median of the risk score. Both Kaplan–Meier curves and multivariate Cox regression were performed in the training cohort.

To quantify the activity of ferroptosis, an established indicator named the ferroptosis potential index (FPI) was obtained from a previous study and calculated for each sample based on gene expression data (Liu et al., 2020). During computation, key ferroptosis regulator genes were first divided into positive and negative groups according to their potential function. Then, the enrichment score of these two groups was calculated separately using the R package GSVA (Hänzelmann et al., 2013). The FPI was defined as the enrichment score of the positive group minus the negative group and was considered a reliable indicator of ferroptosis level.

To further validate the prognostic value of the FRP signature, the risk score was first calculated for each sample based on the formula derived from the training cohort. Then, samples were divided into low-risk and high-risk groups according to the median risk score in each cohort. Kaplan–Meier curves were generated, and univariate Cox regression was performed, and the concordance index (C-index) was applied to evaluate the performance of the FRP signature.

To infer the potential benefit of immunotherapy for each patient, the number of neoantigens and the PD-L1 protein expression of patients in the TCGA-LUAD cohort were extracted from The Cancer Immunome Atlas³ and cBioPortal for Cancer Genomics⁴, respectively (Gao et al., 2013; Charoentong et al., 2017). To determine the chemotherapeutic implications, experimental data for cancer cell drug sensitivity were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC), and the whole prediction process was conducted with the R package pRRophetic (Geeleher et al., 2014). “allSoldTumours” was selected for the tissue type parameter, and the remaining parameters were set to the default values.

BEAS-2B and HCC827 cell lines were purchased from ATCC. BEAS-2B cells were cultured in BEBM basal medium supplemented with 10% fetal bovine serum (FBS), 1,000 U ml^–1 penicillin and 100 μg ml^–1 streptomycin. HCC827 cells were maintained in RPMI 1640 supplemented with 10% FBS, 1,000 U ml^–1 penicillin and 100 μg ml^–1 streptomycin. The cells were cultured with 5% CO₂ at 37°C in a humidified incubator.

Total RNA was extracted from cells using TRIzol reagent (Invitrogen, United States) according to the manufacturer’s instructions. Equal amounts of RNA samples were used for cDNA synthesis with a RevertAid First Strand cDNA synthesis kit (Thermo). SYBR Green-based qRT-PCR on a 7900HT fast real-time PCR system (Applied Biosystems/Life Technologies, Waltham, United States) was performed. The relative mRNA expression levels were calculated by the 2^–ΔΔCt method with normalization to ACTB; the PCR primers are listed in Supplementary Table 1.

All bioinformatic analyses were conducted in R software (Version 3.6.3). Wilcoxon test and Kruskal–Wallis test were performed to statistically test the difference in continuous data for defined subtypes. The Spearman coefficient was used to evaluate the association between variables. Survival analyses, including Kaplan–Meier curves and Cox regression, were performed using the R package survival. The C-index was calculated by R package survcomp. A P-value less than 0.05 was considered statistically significant.

The flow chart of this study is shown in Figure 1. A comprehensive list of genes including 108 that encode ferroptosis drivers and 69 that encode suppressors was extracted from the FerrDb. After removing multirole regulators that would be repeatedly counted, a total of 161 genes remained for screening. Univariate Cox regression was performed, and 37 genes with prognostic significance (p < 0.05) were chosen from the list. Based on the expression level of these genes, three subtypes were identified in the TCGA-LUAD cohort. The clinical information and gene expression patterns are displayed in Figure 1. In subtype A, most of the genes, including both drivers and suppressors, were highly expressed, and in contrast, a substantial portion of ferroptosis regulators was expressed at low levels in subtype B (Figure 2A). The expression pattern in subtype C was somewhere in between the other two subtypes. An association analysis suggested that there was a close connection between the expression levels of the 37 ferroptosis regulators (Figure 2B). Of 24 ferroptosis drivers, nearly one-third (7/24) of the genes was revealed to be favorable prognostic factors, and the remaining genes (17/24) were risk factors for overall survival (OS; Figure 2B). All 13 ferroptosis suppressors were unfavorable prognostic factors (Figure 2B). These results suggested that inhibition of ferroptosis suppressors may improve the survival of LUAD patients. In addition, a significant survival difference (p < 0.0001) was observed among the three ferroptosis subtypes (Figure 2C). Patients with subtype B had the best prognosis compared with patients with the other subtypes, while those with subtype C had the worst prognosis (Figure 2C). The clinicopathological characteristics of the three subtypes were summarized in Table 1.

To further explore the biological pathway underlying each distinct ferroptosis subtype, GSVA was performed to determine pathway enrichment (Subtype A vs. Subtype B; Subtype B vs. Subtype C). As shown in Figure 3A, ferroptosis subtype A was significantly enriched in genes involved with RNA polymerase, pyrimidine metabolism, purine metabolism, cell cycle, and mismatch repair. Ferroptosis subtype B was markedly enriched in genes related to fatty acid and drug metabolism, including bile acid biosynthesis, alpha-linolenic acid metabolism, arachidonic acid metabolism, linoleic acid metabolism, and cytochrome P450 drug metabolism (Figure 3B). In addition, subtype B genes are also involved in cell adhesion and molecules (CAMS), intestinal immune network, and JAK-STAT signaling pathway (Figure 3A). Subtype C genes were predominantly related to pathogenic Escherichia coli infection, the p53 signaling pathway, ubiquitin-mediated proteolysis, and basal transcription factors (Figure 3B). Subsequent evaluation of infiltrated immune cells revealed conspicuous differences in the immune microenvironment among subtypes (Figure 4A). Particularly, as the subtype with the best prognosis, the ferroptosis subtype B genes were prominently enriched in Tfh, NKT, NK, CD4⁺ T, and CD8⁺ T cells (Figure 4A).

To construct a prognostic model related to the defined ferroptosis subtypes, the intersecting DEGs identified in the pairwise comparisons were identified, and a total of 308 genes were selected as initial features (Figure 4B). After the dimension reduction by the Boruta algorithm, PCA was performed on the residual 174 genes, and the results indicated that the expression of these genes had great discriminatory power for the different subtypes (Figure 4C). Thus, these 174 genes constituted a representative collection of genes in the defined ferroptosis subtypes. These demonstrated that the dimensionality reduction process successfully removed unuseful features (n = 134) and retained relevant genes (n = 174). By applying univariate Cox regression, 56 genes with prognostic significance (p < 0.05) were retained for repeated LASSO Cox regression to generate the best model (Supplementary Table 2). A total of 1,000 iterations were conducted, and the distribution of model occurrences is shown in Figure 4D. As presented in Figure 4D, a prognostic signature with 6 genes (EIF5A, CACYBP, CYCS, ANLN, ARNTL2, and PPM1M) had the highest frequencies, with 881, which greatly exceeded those of the other four models and was, therefore, chosen as the FRP signature. The risk score of the FRP signature was calculated as follows: risk score = (0.0413 × EXP_EIF5A) + (0.0068 × EXP_CACYBP) + (0.0254 × EXP_CYCS) + (0.1300 × EXP_ANLN) + (0.1503 × EXP_ARNTL2) – (0.0848 × EXP_PPM1M). After calculating the risk score for each sample in the training cohort, samples were stratified into low-risk and high-risk groups based on the median risk score, and the survival analysis demonstrated that patients in the high-risk group had a significantly shorter overall survival time (log-rank p < 0.001, HR = 2.33, 95% CI = 1.72–3.17) than patients in the low-risk group (Figure 4E). The FPI was applied to quantify the activity of ferroptosis in each group, and the results suggested that the FPI in the high-risk group was significantly (p = 0.044) higher than that in the low-risk group (Figure 4F), which indicated a higher ferroptosis activity in samples that with higher risk. The association of attributes including ferroptosis subtype, risk group, FPI, and survival status of individual patients is presented in an alluvial diagram (Figure 4G).

To decrease the effects of other confounding factors on the relationship between the FRP risk score and OS time, multivariate Cox regression analysis including age, sex, stage, and smoking history was conducted with the TCGA-LUAD cohort. The results proved that age, stage, and risk score (p < 0.001, HR = 2.17, 95% CI = 1.54–3.04) were independent prognostic factors of the OS time for LUAD patients (Figure 5).

The basic information related to the five validation cohorts is described in Table 2. Briefly, a total of 668 cases made up our validation datasets. Consistent and significant differences were observed in validation cohort 1 (log-rank p = 0.001, HR = 2.62, 95% CI = 1.44–4.76), validation cohort 2 (log-rank p = 0.002, HR = 2.62, 95% CI = 1.40–4.89), validation cohort 3 (log-rank p = 0.005, HR = 2.75, 95% CI = 1.32–5.74), validation cohort 4 (log-rank p = 0.014, HR = 2.39, 95% CI = 1.17–4.92), and validation cohort 5 (log-rank p < 0.001, HR = 2.59, 95% CI = 1.57–4.26; Figures 6A–E). In addition, the C-index was calculated to assess predictive performance and the C-index was 0.663 for the training cohort, 0.641 for validation cohort 1, 0.694 for validation cohort 2, 0.666 for validation cohort 3, 0.647 for validation cohort 4, and 0.685 for validation cohort 5 (Figure 6F). Statistical significance (p < 0.05) was observed between the actual value of the C-index and the null hypothesis (C-index = 0.5) in all datasets (Figure 6F).

The abundance of infiltrated immune cells was estimated using CIBERSORT and ImmuCellAI for each sample, and a heat map was generated to illustrate the association between immune cell infiltration and the FRP risk score (Figure 7A). The top five immune cell types with significant changes in abundance between risk groups were CD4⁺ T cells (p < 0.001), nTregs (p < 0.001), Tfhs (p < 0.001), Tems (p < 0.001), and monocytes (p < 0.001). The detailed results of these differential analyses are presented in Supplementary Table 3. Furthermore, an association analysis of the FRP risk score and potential immunotherapeutic biomarkers, including TMB, PD-L1 protein, and neoantigen, was performed. The results suggested that TMB in the high-risk group was significantly (p < 0.001) higher than that in the low-risk group (Figure 7B). A correlation analysis revealed that the FRP risk score was significantly associated (R = 0.365, p < 0.0001) with TMB (Figure 7C). In addition, the expression level of the PD-L1 protein (p < 0.001) and the numbers of subclonal neoantigens (p = 0.027), clonal neoantigens (p < 0.001), and total neoantigens (p < 0.001) were all significantly higher in the high-risk group than in the low-risk group (Figures 7D–G). These results implied that patients in the high-risk group may potentially benefit from immunotherapy.

Considering that chemotherapy constitutes an important component of LAUD treatment, drug susceptibility data including 138 compounds from the GDSC were applied to the training group to predict drug sensitivities for each sample in the TCGA-LUAD cohort. Based on the significance level, the top 10 compounds that might benefit patients in the high-risk group were identified: A.443654, obatoclax mesylate, epothilone B, GW843682X, BI.2536, TW.37, thapsigargin, ZM.447439, PF.562271, and S-Trityl-L-cysteine (Figures 8A–J). These results provide a novel perspective for exploring new drug candidates. A detailed list showing the significance value for all 138 compounds is available in Supplementary Table 4.

To further investigate the expression pattern in vitro, six genes that constituted the FRP signature were quantified by qPCR and those in normal bronchial epithelial cells (BEAS-2B cells) and LUAD cells (HCC827 cells) were compared. Of these genes, EIF5A and ARNTL2 were significantly overexpressed in tumor cells (p < 0.05; Figures 9A,E). In contrast, the expression levels of CYCS, ANLN, and PPM1M were significantly elevated in normal lung cells (p < 0.05; Figures 9C,D,F). No significant difference in CACYBP expression was observed between the two cell lines (Figure 9B).