Construction and validation of a 15-gene ferroptosis signature in lung adenocarcinoma

Acquisition of ferroptosis-related genes

We downloaded the list of ferroptosis pathway genes (map04216) from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (https://www.genome.jp/kegg/pathway.html) and designated these genes as ferroptosis-related genes. Iron metabolism-associated genes were retrieved from the R-HAS-917937 pathway in the Reactome pathway database (https://reactome.org/) and the cellular iron ion homeostasis pathway in the AmiGO2 database (http://amigo.geneontology.org/amigo). After systematically searching and analyzing the original document, we discarded the genes that do not have a modulatory effect on ferroptosis. In addition, we collected and integrated newly reported ferroptosis-related genes for subsequent research.

Data collection

We obtained the level 3 mRNA expression profiles and corresponding clinical data of LUAD patients from the TCGA database (https://portal.gdc.cancer.gov/) and GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Rousseaux et al., 2013; Schabath et al., 2016) up to September 20, 2020. Samples with a follow-up time of less than 30 days or lack of prognostic data were excluded. Since the TCGA and GEO databases are public to researchers and we completely abided by the publication guidelines as well as the policies of access to the database, ethical review and approval were not required.

Identification of differentially expressed ferroptosis genes

The “limma” R package was used for the normalization of gene expression matrixes. Then, we matched the mRNA sequencing data with the ferroptosis-related gene list and performed differential expression analysis between the tumor tissues and normal tissues in the TCGA cohort by using false discovery rate (FDR) < 0.05 as the threshold. Thus, ferroptosis-related differentially expressed genes (DEGs) were identified. A heatmap and volcano plot to visualize the DEGs were generated by the “pheatmap” R package. In addition, the protein-protein interaction (PPI) network of the ferroptosis-related DEGs was analyzed in the Search Tool for the Retrieval of Interacting Genes (STRING) online database and visualized in Cytoscape 3.8.2 software. We also generated a correlation network of the DEGs using the “igraph” R package.

Establishment and validation of the prognostic model

DEGs with prognostic value were screened by univariate Cox regression analysis with a P value less than 0.05. Then, the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was conducted to establish the signature. In brief, we used the normalized expression data as the independent variable and the overall survival (OS) data of the patients in the TCGA cohort as the response variable to perform the LASSO algorithm for the shrinkage of variables by the “glmnet” R package. Tenfold cross-validation was used to narrow the number of candidate genes and identify the penalty parameter (λ), which corresponds to the lowest position of the likelihood deviance curve. We further evaluated the prognostic value of the LASSO genes by Kaplan–Meier (K–M) survival analysis. The K–M survival analysis was conducted using the optimal cutoff determined by the “surv_cutpoint” function in the “survival” R package.

The risk score of each patient was calculated with the following formula: risk score = Σ (ExpmRNAn × βmRNAn). The patients were stratified by the median risk score into high- and low-risk groups. After that, we performed principal component analysis (PCA) to evaluate the discriminatory ability by the “status” R package. We evaluated the predictive ability of the signature for OS through K–M survival analysis as well as time-dependent receiver operating characteristic (ROC) curve analysis conducted by the “survminer” and “survivalROC” R packages. Finally, we performed multivariate Cox regression analysis to identify independent risk factors for OS. In this study, the TCGA cohort was used as the derivation cohort, and the GSE72094 and GSE30219 cohorts were used as external validation cohorts.

Construction and evaluation of the predictive nomogram

We constructed a nomogram to predict the survival probabilities at 1, 2 and 3 years in the TCGA cohort by integrating all the independent risk factors using the “rms” R package. The patients were stratified into different risk groups, and K–M survival analysis was conducted to analyze the OS difference between the different risk groups. Then, we calculated the C index and performed time-dependent ROC analysis to further validate the prediction accuracy of the nomogram. Moreover, a calibration plot was used to evaluate the consistency between the predicted survival probability and the real observation. The GSE72094 cohort was used for the external validation of the nomogram.

Correlation with immune status

Using the “GSVA” R package, we calculated the enrichment score of 16 immune-related cells as well as 12 immune-related functions for each patient by single-sample gene set enrichment analysis (ssGSEA). In brief, using a set of genes that correspond to a particular immune cell or immune function, the enrichment score was calculated in the gene expression matrix through the ssGSEA algorithm, and the enrichment scores were normalized for subsequent analysis. Then, the enrichment scores for diverse immune cells and functions of patients in different risk groups were compared to illustrate the potential correlation between ferroptosis and immune status. The R script used for this part of the analysis is provided in the GitHub website (https://github.com/guangxu0109/ssGSEA.git).

Functional enrichment analysis

To gain insight into the molecular mechanisms of these ferroptosis-related genes, we performed Gene Ontology (GO) and KEGG enrichment analyses of the DEGs between the high-risk group and low-risk group in the TCGA cohort, which were screened by the thresholds of |log2 fold change (FC)| ≥ 1 and FDR < 0.05. GO enrichment analysis was performed by the “clusterProfiler” R package, and KEGG enrichment analysis was performed by gene set enrichment analysis (GSEA) using GSEA 4.0.2 software. The pathways with a P value less than 0.05 were considered to be significantly enriched.

Statistical analysis

All statistical analyses were conducted in R software version 4.0.2. The Mann–Whitney test was conducted to compare gene expression levels between tumor tissues and normal tissues as well as to compare the ssGSEA enrichment scores between different risk groups. A P value less than 0.05 (if not otherwise specified) was considered statistically significant.

Twenty-two genes correlated with prognosis were identified in the TCGA cohort

We generated a flowchart to describe the design of this study, as shown in Fig. 1. A total of 948 LUAD patients from the TCGA database (N = 477), GSE72094 cohort (N = 386) and GSE30219 cohort (N = 85) were used for subsequent analysis in this study. Table S1 shows the basic characteristics of these patients. A total of 125 ferroptosis-associated genes (as listed in Table S2) were identified to intersect with the mRNA expression matrix of the TCGA and GEO databases. Ninety-seven genes were confirmed to be differentially expressed between tumor tissues and nontumorous tissues (62 upregulated and 35 downregulated) in the TCGA cohort (Figs. 2A, 2B and Table S3). The PPI network and the correlation network of the DEGs are shown in Fig. S1. Twenty-two genes were identified to be correlated with prognosis through univariate Cox regression analysis (Fig. 3).

Establishment and assessment of the prognostic ferroptosis signature in the TCGA cohort

A signature consisting of 15 genes was established based on the minimum λ value identified by the LASSO algorithm in the TCGA cohort (Fig. S2). The risk score of each patient was calculated with the following formula: (−0.037 × expression level of AGER) + (0.261 × expression level of CISD1) + (−0.019 × expression level of DPP4) + (0.173 × expression level of EGLN1) + (0.077 × expression level of FANCD2) + (−0.253 × expression level of GLS2) + (−0.233 × expression level of ISCU) + (0.087 × expression level of ITGA6) + (0.024 × expression level of ITGB4) + (0.118 × expression level of KRAS) + (0.109 × expression level of NEDD4) + (−0.048 × expression level of PEBP1) + (−0.133 × expression level of SLC11A2) + (0.013 × expression level of TFAP2A) + (0.208 × expression level of VDAC1). We further performed K–M survival analysis for the 15 ferroptosis-related genes in the prognostic signature. The results confirmed that all these genes were significantly correlated with OS (Fig. S3). Among them, AGER, DPP4, GLS2, ISCU, PEBP1 and SLC11A2 were identified to be protective factors for OS, while the remaining factors (CISD1, EGLN1, FANCD2, ITGA6, ITGB4, KRAS, NEDD4, TFAP2A and VDAC1) were risk factors, which was consistent with the results of univariate Cox analysis, as shown in Fig. 3.

The patients in the TCGA cohort were divided into different risk groups by the median risk score (Table S4). The PCA plot in Fig. 4A shows that patients in different risk groups were clearly distributed in two directions, indicating that the gene signature had good discriminatory power. The K–M curve showed that the OS of patients in the high-risk group was significantly lower than that of patients in the low-risk group (Figs. 4B, 4C). In addition, OS was significantly decreased with increasing risk scores (Fig. 4D). The area under the curve (AUC) value of the risk score model to predict OS was 0.775, which was the highest among all the risk factors, indicating the superior predictive ability of the prognostic signature (Fig. 5E).

External validation of the ferroptosis-related gene signature

To further validate the robustness of the signature, we applied the 15-gene ferroptosis signature to the GSE72094 and GSE30219 cohorts. The patients in the two cohorts were stratified into high- and low-risk groups by the median risk score (Tables S5 and S6). In the GSE72094 cohort, the PCA plot showed that patients in different risk groups were distributed in two directions (Fig. 5A). The K–M curve indicated that patients in the high-risk group had significantly shorter OS than those in the low-risk group, and OS was significantly decreased with increasing risk scores (Figs. 5B–5D). The AUC value of the risk score model was 0.711, which was the highest among all the risk factors (Fig. 5E). In the GSE30219 cohort, the PCA plots showed that patients in different groups were distributed in two directions (Fig. 5F). The K–M survival curve confirmed a significantly worse OS for patients in the high-risk group than for those in the low-risk group (Fig. 5G). OS significantly decreased with increasing risk scores (Figs. 5H, 5I). The AUC value of the risk score model was 0.877, which was the highest among all the risk factors (Fig. 5J). All these results indicated that the 15-gene ferroptosis signature had a robust predictive performance.

Independent prognostic value of the ferroptosis-related gene signature

To further assess the independent prognostic value of the 15-gene ferroptosis signature, we performed univariate and multivariate Cox regression analyses on all characteristics, including age, sex and TNM stage, and the risk score model based on the ferroptosis signature in both the derivation cohort and two validation cohorts. The results are presented in Table 1. In the TCGA cohort, TNM stage and risk score model were identified to be independent prognostic factors for OS. In the GSE72094 cohort, sex, TNM stage and risk score model were independent prognostic factors for OS. The risk score model was the only independent prognostic factor for OS in the GSE30219 cohort.

Given that the risk score model was not the only independent prognostic factor in the TCGA and GSE72094 cohorts, we further performed subgroup survival analysis to verify whether the 15-gene ferroptosis signature could be independent of other risk factors to predict prognosis. In the TCGA cohort, patients in stage I–II and stage III–IV were stratified into high-risk and low-risk groups by the median risk score. K–M curves showed that the OS of the high-risk group was significantly worse than that of the low-risk group regardless of whether patients were stage I-II or stage III-IV, indicating that the 15-gene ferroptosis signature could predict OS independent of TNM stage (Figs. 6A, 6B). Similarly, in the GSE72094 cohort, the OS of the high-risk group was significantly worse than that of the low-risk group regardless of whether the patient was male or female, indicating that the signature could predict OS independent of sex (Figs. 6C, 6D). For the TNM stage subgroup, the OS of patients in the high-risk group was significantly shorter than that of patients in the low-risk group except for patients in stage III and stage IV, which was likely due to the small sample size of patients in these groups (Figs. 6E, 6F). Overall, the 15-gene ferroptosis signature could predict OS independent of other clinical characteristics.

Construction and assessment of the predictive nomogram

Based on the independent risk factors identified in the TCGA cohort, we constructed a nomogram model to predict the OS probabilities of patients at 1, 2 and 3 years for risk assessment and earlier intervention to improve patient survival time (as shown in Fig. 7A). As depicted in Fig. 7B, a significantly poorer prognosis was observed in the high-risk group. The C index for the nomogram model to predict OS was 0.770 (95% CI [0.724–0.816]), indicating that the nomogram model had a robust predictive accuracy. We further evaluated the discrimination and calibration of the nomogram by time-dependent ROC curves and calibration plots for 1-year and 3-year OS. The 1-year and 3-year calibration plots indicated that the predicted survival was highly consistent with the actual survival in the TCGA cohort (Figs. 7D, 7E). As shown in Figs. 7H and 7I, the AUC values of the nomogram for predicting 1-year and 3-year OS were 0.812 and 0.757, respectively. All the AUC values of the nomogram were superior to those of other independent risk factors, indicating that the nomogram had a better predictive performance.

The nomogram was further validated in the GSE72094 cohort, while TNM stage information was not given in the GSE30219 cohort. The results showed that the OS of patients in the high-risk group was significantly shorter than that of patients in the low-risk group (Fig. 7C). The C index of the nomogram model was 0.708 (95% CI [0.660–0.756]). The calibration plots indicated great consistency between the predicted survival rate and the real observation (Figs. 7F, 7G). The AUC values of the nomogram for predicting 1-year and 3-year OS were 0.735 and 0.758, respectively, which were also superior to those of all other independent risk factors (Figs. 7J, 7K).

Functional analysis of the 15-gene ferroptosis signature

GO and KEGG enrichment analyses were performed to elucidate the possible biological functions and pathways involved in the 15-gene ferroptosis signature. The results of GO enrichment analysis indicated that the DEGs between the high-risk group and the low-risk group were mainly enriched in pathways of the cell cycle and immune response, such as chromosome segregation, mitotic nuclear division, mitotic sister chromatid segregation, antimicrobial humoral response and humoral immune response (Fig. 8A and Table S7). As shown in Fig. 8B and Table 2, the KEGG enrichment analysis revealed that the top five pathways enriched in the high-risk group were the cell cycle, ubiquitin-mediated proteolysis, oocyte meiosis, homologous recombination and p53 signaling pathways. The top five pathways enriched in the low-risk group were the arachidonic acid metabolism, primary bile acid biosynthesis, alpha linolenic acid metabolism, asthma, and intestinal immune network for IgA production pathways.

Correlation with immune status

Given that pathways of the immune response were significantly enriched in the 15-gene ferroptosis signature, we further sought to investigate the correlation between immune status and risk score of the ferroptosis signature. The enrichment results are presented in Tables S8–10. Then, we analyzed the difference in the enrichment score of immune cells and immune functions between the high-risk and low-risk groups. By comparing the enrichment results of the three cohorts, we found that enrichment scores of activated dendritic cells (aDCs), dendritic cells (DCs), immature dendritic cells (iDCs), mast cells and neutrophils were significantly lower in the high-risk group than in the low-risk group (Figs. 9A, 9C, 9E). Regarding immune-related functions, higher enrichment score of antigen-presenting cell (APC) inhibition was found in the high-risk groups of all three cohorts, while a higher enrichment score of the type II IFN (IFN-γ) response was found in the low-risk group (Figs. 9B, 9D, 9F). All these results indicated a close correlation between risk score of the ferroptosis signature and immune cells as well as immune functions.