We collected transcriptome profiles of lung adenocarcinoma (LUAD) available in the TCGA database (https://portal.gdc.cancer.gov/) on 6 July 2021.

Preprocessed and aligned RNA-Seq samples were downloaded by selecting HTSeq-Counts as the workflow type on the portal. Clinical and pathological information related to the TCGA–LUAD cohort was retrieved from the cBioportal website (https://www.cbioportal.org) with the “cdgsr” package (Cerami et al., 2012; Jacobsen, 2015).

Meanwhile, we collected two microarray datasets and corresponding clinical information, including GSE68465 and GSE72094, from Gene Expression Omnibus using the “GEOquery” package (https://www.ncbi.nlm.nih.gov/geo/) (Davis and Meltzer, 2007). These two affymetrix microarrays were preprocessed using the RMA method (R package “affy”) (Gautier et al., 2004).

Removing samples without overall survival (OS) information or with an OS time of 0 and converting the TNM stage to AJCC staging groups, the TCGA–LUAD cohort (N = 306) was used as the training cohort, whereas GSE68465 (N = 441) and GSE72094 (N = 398) were used as the validation cohorts (Supplementary Table S1).

The specific data processing and research flow are shown in Figure 1.

One ferroptosis-related gene set was downloaded from FerrDb (http://www.zhounan.org/ferrdb/index.html), the world’s first database of ferroptosis regulators, markers, and associations (Zhou and Bao, 2020). There were 108 drivers annotated as genes that promoted ferroptosis, 69 suppressors that prevented ferroptosis, and 111 markers that indicated the occurrence of ferroptosis. Removing multi-annotated genes and selecting ferroptosis-related genes measured by using all the three cohorts, 209 genes were included in the ferroptosis-related gene set (Supplementary Figure S1).

The mean expression value of replicated genes was calculated. We defined the combination of two ferroptosis-related genes (FRG-1 and FRG-2) as a ferroptosis-related gene pair (FRGP). In every specific sample, 209 FRGs were compared by the gene expression level with each other to build 21736 FRGPs and generate a score for each FRGP. An FRGP score of 1 was assigned if FRG-1 was less than FRG-2. Otherwise, the FRGP score was 0. This gene-pair-based approach calculated the FRGP score based totally on the gene expression value of each individual sample and could be applied without normalization.

To reduce biases and be meaningful for subsequent analyses, some FRGPs counting 0s or 1s in more than 80% samples were filtered out. Finally, we got 6374 FRGPs as columns and 306 samples as rows to form the training matrix.

The least absolute shrinkage and selection operator (LASSO) is a statistical method to reduce data dimensionality. We applied the LASSO regression operation with the R package “glmnet” and “survival” in the TCGA FRGP matrix to construct a prognostic ferroptosis-related gene-pair index (FRGPI)(Simon et al., 2011). Insignificant variables whose coefficients became zero and any collinear variables were removed. 10-fold cross validation (CV) divided data into ten equal parts, nine parts as the training set and the remaining one part as the validation part. When partial likelihood deviance was the smallest, we got the minimum of lambda and nine gene pairs as our best FRGPI model. The prognostic index signature is expressed as FRGPI risk score = ∑n_i (FRGP_i *coef_i) (i = 1,2,3 … … 9, n = 9, where n is the number of FRGPs, FRGPi is the score (0 or 1) of the ith FRGP, and coefi is the regression coefficient of the ith FRGP).

First, we used the R package “survivalROC” to draw ROC curves and calculated the AUC values in the training and validation cohorts. AUC values greater than 0.5 and closer to 1 indicated the prognostic ability of FRGPI.

Second, we used the Kaplan–Meier method to compare survival outcomes between high- and low-FRGPI risk score groups in training and validation cohorts. The optimal cut-off value, determined based on the best balance of sensitivity and specificity to achieve the best AUC in the training cohort, was investigated using the ROC curves with the R package “survivalROC,” and “survminer” (Weiss et al., 2003).

After verification of the FRGPI significantly stratifying patients into low- and high-risk groups, we performed univariate and multivariate Cox hazard analyses to validate FRGPI as an independent prognostic factor. The hazard ratio (HR) in survival analyses less than 1 meant that the presence of the factor was protective, whereas the hazard ratio more than 1 was harmful.

After constructing and validating FRGPI, we compared the expression of ferroptosis-related genes between high- and low-FRGPI groups in the TCGA cohort and GEO validation cohorts, using the Wilcoxon rank-sum test and the reshape2 package (Wickham, 2007). Meanwhile, we drew the boxplots for visualization using the ggplot2 package (Wickham, 2009).In addition, we compared MKI67 expression between high- and low-FRGPI groups in the three cohorts. MKI67 encodes a nuclear protein Ki-67, which is a commonly used marker for cell proliferation.

Gene set enrichment analysis (GSEA) determines whether the gene sets, not single genes, present differences between different biological status groups and verifies that the gene sets are enriched in one specific clinical group (Subramanian et al., 2005). The gene sets are predefined by previous experiments and function annotations.

We chose gene sets from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway which is a collection of pathway maps representing molecular interactions, reactions, and relation networks (Kanehisa et al., 2016).

We performed GSEA–KEGG analyses and drew enrichment plots in the TCGA cohort and GEO validation cohorts, using “c2.cp.kegg.v7.4.entrez.gmt” and the R package “clusterProfiler,” and “ggplot2,” defining p value < 0.05 and q value < 0.05 as the filtering criteria (Yu et al., 2012).

Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) is a deconvolution method for characterizing the cell composition of complex tissues from gene expression profiles (Newman et al., 2015). We made 22 kinds of immune cells as the target characterizing composition. We operated the “CIBERSORT” algorithm with the “Leukocyte signature matrix” (Newman et al., 2015). Based on the composition of 22 kinds of immune cells, we screened out what kind of immune cells differently infiltrated between high- and low-FRGPI groups.

All statistical analyses were based on R Programming Language software (Rx64 3.3.3). The online website “www.genome.jp/kegg/” was used for GSEA–KEGG gene sets. The website “cibersortx” offered thoughts for CIBERSORT analyses in R software.

First, we downloaded one ferroptosis-related gene set including 209 genes. Two genes as one couple, any couple of the 209 genes could be created. So, 209 genes formed 21736 FRGPs in every individual sample. To calculate each FRGP in every sample, if the FRG-1 value was less than the FRG-2 value, the FRGP score was 1. Otherwise, the FRGP score was 0. Then, we filtered out FRGPs with constant values (0 or 1) in TCGA or GEO datasets. Finally, we got 6374 FRGPs in every individual sample. In TCGA training cohort, we used 306 samples as rows and 6374 FRGP scores as columns to make up the training FRGP matrix.

The training FRGP matrix was used for evaluating the relationship between FRGPs and overall survival rates applied using the LASSO regression operation. The lambda and coef diagrams of FRGPs (Figure 2A) were plotted using the LASSO algorithm. With the increase in the lambda value, the coefficients of some FRGPs decreased to zero, which meant that the scores of these FRGPs did not affect the model. We then used 10-fold CV to calculate the partial likelihood deviance of the model (Figure 2B). The minimum deviance exported the best model. The best model included nine gene pairs. Finally, we constructed a prognostic FRGPI signature with nine FRGPs and the corresponding coef values (Figure 2C).

The time-dependent ROC curves with AUC values of the training and validation cohorts are all presented in Figures 2D,F,G. All the AUC values were more than 0.5 and even greater than 0.7 in training and validation cohorts, indicating that FRGPI had a favorable prognostic ability.

In the training cohort, the AUC value reached 0.773 for 1606 days. We defined the optimal cut-off 0.081 in the curve with the best AUC and “1606” as the time point (Figure 2E). 0.081 was used as a cut-off for FRGPI to stratify patients into the low- or high-FRGPI risk score group.

We then performed Kaplan–Meier curves between high- and low-FRGPI groups. All three curves showed that the low-FRGPI group had a significantly better prognosis than the high-FRGPI group (p < 0.01, Figures 3A–C).

To further validate FRGPI as a prognostic biomarker, the low-FRGPI group also had significantly better prognoses than the high-FRGPI group for 2 and 5 years (p < 0.05, Supplementary Figures S2A–F). In the early-stage LUAD, all the Kaplan–Meier curves showed that the low-FRGPI group had a significantly better prognosis (p < 0.01, Supplementary Figures S2G–I). Additionally, ever-smokers with low FRGPI scores owned a better prognosis (p < 0.05, Supplementary Figure S3). Overall, FRGPI successfully stratified the risk of LUAD patients in all the training and validation sets.

Univariate regression results showed that FRGPI was statistically significant in the training and validation cohorts (p < 0.01, Figures 3D–F). The multivariate regression results showed that FRGPI was an independent prognostic factor in all three cohorts (p < 0.05, Figures 3G–I). All univariate and multivariate results revealed that the high-FRGPI group matched with a worse prognosis, with the HR and 95% confidence interval HR of FRGPI more than 1.

In the multivariate regression, it was shown that our FRGPI could be as good as stage III vs. stage I in stratifying patients. Patients whose cancer has progressed to stage III, especially stage IIIB, could hardly get radical surgical therapy. The median PFS of these patients is about 10 months (Ryan et al., 2019).

Our FRGPI included 16 ferroptosis-related genes. In addition, SLC7A11 as the main target of FINs was added. We compared the expression of 17 ferroptosis-related genes and MKI67 between high- and low-FRGPI groups (p < 0.05, Figures 4A–C). In training and validation cohorts at the same time, solute carrier family 2 member 1 (SLC2A1), gelsolin-like actin-capping protein (CAPG), ribonucleotide reductase regulatory subunit M2 (RRM2), SLC7A11, and MKI67 were significantly up-regulated in the high-FRGPI group. On the contrary, we found GLS2 and phosphatidylethanolamine-binding protein 1 (PEBP1) were down-regulated in the high-FRGPI group in all three cohorts. The high-FRGPI group was marked with significantly higher MKI67 expression and exhibited higher cancer proliferation potential.

FRGPI separated patients into high and low groups. The differences of enriched GSEA–KEGG pathways between the two groups are shown in Figures 4D–F. There was no enriched KEGG pathway in the low-FRGPI group in the training and validation cohorts at the same time. Conversely, we found that “KEGG CELL CYCLE,” “KEGG DNA REPLICATION,” “KEGG PROTEASOME,” and “KEGG P53 SIGNALING PATHWAY” were significantly enriched in the high-FRGPI group in all three cohorts (p < 0.05, Figures 4D–F). We inferred that the consistently enriched pathways in the high-FRGPI group in part played important roles in the worse prognosis.

The differences of infiltration levels of 22 kinds of immune cells are shown in Figures 4G–I. We found that M1 macrophages were significantly up-regulated in the high-FRGPI group in all three cohorts; however, activated CD4⁺ cells, CD8⁺ cells, dendritic cells, plasma cells, and natural killer cells did not present different infiltration levels. There were few differences in adaptive immune responses between high- and low-FRGPI groups.