Dysregulated ferroptosis-related genes indicate potential clinical benefits for anti-PD-1/PD-L1 immunotherapy in lung adenocarcinoma.

BACKGROUND: Ferroptosis is an iron-dependent programmed cell death mechanism that influences the development of malignancy. Lung adenocarcinoma (LUAD) is the most common type of lung cancer with no known cure. Anti-PD-1/PD-L immunotherapy is effective for patients with partial LUAD. Therefore, there is an immediate requirement of novel markers to predict the individualised benefits of immunotherapy.
METHODS: We manually collected the ferroptosis-related gene (FERG) set and employed the Wilcoxon rank-sum test to identify the differentially expressed FERGs. Subsequently, we constructed a recursive partitioning and regression tree (RPART) model to predict the benefits of anti-PD-1/PD-L1 immunotherapy. Subsequently, the ROC curve and AUC were used to evaluate the model efficiency in an independent dataset.
RESULTS: In this study, we found that the dysregulated FERGs were closely associated with multiple metabolic processes in LUAD. Furthermore, we identified three ferroptosis-related tumour subtypes (F1, F3 and F3). The F3 subtype exhibited higher immunoactivity and lower tumour purity, mutation count and aneuploidy and had better survival outcomes compared with the other two subtypes, implying that FERGs played an important role in intertumoral immune heterogeneity. We further explored the role of FERGs in the anti-PD-1/PD-L1 immunotherapy. We identified a set of three-FERGs signature (CD44, G6PD and ZEB1) that acted as a promising indicator (AUC = 0.697) for the prediction of the benefits of anti-PD-1/PD-L1 immunotherapy.
CONCLUSION: Ferroptosis, as emerging programmed cell death mechanism, was associated with cancer development. We used ferroptosis-related genes to predict the immunotherapy benefits that may facilitate the development of individualised anti-cancer treatment strategies.

INTRODUCTION

Lung adenocarcinoma (LUAD) is a type of non‐small cell lung cancer that accounts for approximately 40% of all lung cancer cases. , Immunotherapy has achieved promising clinical benefits as an emerging cancer treatment strategy for patients with partial LUAD. Anti–PD‐1/PD‐L1 immunotherapy is one of the most common immunotherapy strategies, which has demonstrated its activity as a single agent in the treatment of lung cancer. However, the individualised differences in immunotherapy remain a pressing issue. Therefore, there is an immediate requirement for novel biomarkers to predict the response of PD‐1‐guided immunotherapies.

Ferroptosis is a form of programmed cell death that is iron‐dependent and distinct from apoptosis, cell necrosis and autophagy. It influences many aspects of cell proliferation, developmental defects, malignancy development and immune regulation. Previous studies have revealed that ferroptosis participated in the development of cancer. Foetal metabolic imbalances resulting from glutathione (GSH) depletion or GSH peroxidase 4 inactivation are the enforcers of iron toxicity within cancer cells. , In addition, small molecule‐induced ferroptosis may act as a strong inhibitor of tumour growth and further improve the sensitivity of chemotherapeutic drugs. , Evidence suggests that ferroptosis plays an important role in cancer treatment. However, the role of ferroptosis in LUAD remains unclear.

In the present study, we systemically characterised dysregulated ferroptosis‐related genes (FERGs) between LUAD tumours and adjacent normal tissues. FEGRs were capable of participating in multiple metabolism‐related processes. Furthermore, we demonstrated a somatic mutation landscape of FERGs. Somatic mutation was an important influence factor of the expression levels of FERGs, wherein TP53, KEAP1 and ZEB1 were the most affected FERGs. Furthermore, we identified the ferroptosis‐related tumour subtypes (F1, F3 and F3). The F3 subtype exhibited higher immunoactivity and lower tumour purity, mutation count and aneuploidy than those of F1 and F2. Notably, patients in the F3 subtype also showed better survival outcomes and radiotherapy benefits. In addition, we also explored the role of FERGs in anti–PD‐1/PD‐L1 immunotherapy by employing a recursive partitioning and regression tree (RPART) model (Figure S1). The FERGs CD44, G6PD and ZEB1 could act as immunotherapy benefit indicators (area under the receiver operating characteristic (ROC) curve [AUC] = 0.697). These FERGs were generally involved in cancer‐related and metabolic processes.

MATERIALS AND METHODS

Data collection

The Cancer Genome Atlas (TCGA) level 3 gene expression and somatic mutation data of patients with LUAD were downloaded from the Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/). The clinical data of corresponding samples were obtained from the cBioPortal database (https://www.cbioportal.org/). Two anti–PD‐1/PD‐L1 immunotherapy datasets of lung cancer GSE136961 and GSE119144 were obtained from the Gene Expression Omnibus (GEO). The datasets used in the study are listed in Table 1.

TABLE 1

Clinical characteristics of the patients in this study

	TCGA	GSE136961	GSE135222
Sample size	516	21	27
Age (mean ± SD)	65.22 ± 9.97		62.07 ± 8.79
Gender
Male	240	15	22
Female	276	6	5
Radiotherapy	61
Anti‐PD‐1/PD‐L1 immunotherapy
DCB		9	8
NDB		12	19

Collection of ferroptosis‐related genes

We manually collected 60 FERGs from recently published papers, , , , including ABCC1, ACACA, ACO1, ACSF2, ACSL3, ACSL4, AIFM2, AKR1C1, AKR1C2, AKR1C3, ALOX12, ALOX15, ALOX5, ATP5MC3, CARS, CBS, CD44, CHAC1, CISD1, CRYAB, CS, DPP4, EMC2, FADS2, FANCD2, FDFT1, FTH1, G6PD, GCLC, GCLM, GLS2, GOT1, GPX4, GSS, HMGCR, HMOX1, HSBP1, HSPB1, IREB2, KEAP1, LPCAT3, MT1G, NCOA4, NFE2L2, NFS1, NOX1, NQO1, PEBP1, PGD, PHKG2, PTGS2, RPL8, SAT1, SLC1A5, SLC7A11, SQLE, STEAP3, TFRC, TP53 and ZEB1.

Differential expression analysis

For the gene expression data obtained, first, we removed the genes with fragments per kilobase of exon per million mapped fragments (FPKM) of 0 in ≥70% of the tumour or normal samples. Subsequently, we identified the differentially expressed genes (DEGs) between LUAD tumours and adjacent normal tissues using the Wilcoxon rank‐sum test. p‐values were adjusted using the Benjamini and Hochberg procedure. Genes with adjusted p‐value <0.01 were defined as DEGs in LUAD.

Mutation analysis

We first integrated the somatic mutation data from four platforms (MuSE, MuTect2, SomaticSniper and VarScan2). Mutations in FERGs were considered for the following analysis. The R package ‘maftools’ was used to perform mutation‐related analyses.

Identification of ferroptosis‐related tumour subtypes

We performed consensus clustering on LUAD samples based on the expression levels of deferentially expressed FERGs using the R package ‘ConsensusClusterPlus’. When the number of clusters (k) increased from 2 to 10, there was a minimal increase in the area under the cumulative distribution function curve in consensus clustering. When k reached 3, there was an optimal number of clusters (F1, F2 and F3).

Evaluation of tumour immune cell infiltration levels

We characterised the immune cell infiltration abundance (B cells, CD4‐T cells, CD8‐T cells, macrophages, neutrophils and dendritic cells) using TIMER. In addition, the tumour purity and immune scores were evaluated using the R package ‘ESTIMATE’. We also compared the immune cell infiltration levels among three ferroptosis‐related tumour subtypes using the Wilcoxon rank‐sum test. p‐Values were adjusted using the Benjamini and Hochberg procedure.

Construction of decision model for predicting immunotherapy benefits

We trained an RPART model to predict the immunotherapy benefits in LUAD based on the expression levels of deferentially expressed FERGs in the dataset GSE136961. A three‐FERG signature (FER signature) was identified as the immunotherapy benefit indicator. To confirm the effectiveness of the signature, we extracted FER features among samples based on the principal component analysis of the FER signature in an independent test set GSE119144. Subsequently, the first and second principal components were used to construct the FES score; FES = PC1 + PC2. Finally, the ROC curve and AUC were used to evaluate the efficiency of the model.

Functional enrichment analysis

To annotate the potential biological functions of the FER signature, we performed functional enrichment analysis using Metascape. We first obtained the associated genes with each gene in the FER signature using the protein–protein interaction data from STRING (https://string‐db.org/). For these genes, the pathway and process enrichment analyses were performed using the following ontology sources: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO) biological processes, Reactome gene sets and canonical pathways. Terms with p < 0.01, a minimum count of 3 and an enrichment score >1.5 were collected and grouped into clusters based on their membership similarities. p‐Values were calculated based on the accumulative hypergeometric distribution and adjusted using the Benjamini and Hochberg procedure.

Statistical analysis

All analyses were conducted using the R (v3.6.3) software.

RESULTS

Dysregulated ferroptosis‐related genes were associated with metabolic processes in LUAD

To explore the expression levels of 60 FERGs, we performed a differential expression analysis between LUAD tumours and adjacent normal tissues (as mentioned in Materials and Methods). In total, 45 significantly and differentially expressed FERGs (adjusted p < 0.01) were identified, which included 29 upregulated and 16 downregulated FERGs (Figure 1A). Furthermore, we annotated their biological functions based on Metascape. We found that both upregulated and downregulated FERGs were mainly enriched in ferroptosis (WP4313) and various multiple metabolism‐related processes. For example, upregulated FERGs exhibited associations with response to oxidative stress (GO:0006979), NADP metabolic process (GO:0006739) and ion homeostasis (GO:0050801) (Figure 1B). However, downregulated FERGs were enriched in lipid metabolism (R‐HSA‐556833) and cholesterol metabolism with the Bloch and Kandutsch–Russell pathways (WP4718) and nucleoside bisphosphate metabolism (GO:0033865) (Figure 1C). These results indicated that dysregulated FERGs were associated with metabolic processes and might play an important role in LUAD.

FIGURE 1

Differentially expressed ferroptosis‐related genes in lung adenocarcinoma. (A) CircosPlot showing the expression levels of ferroptosis‐related genes in LUAD and adjacent normal tissues. From inner to outer, the mean expression levels of normal and tumour tissues are demonstrated. Differentially expressed FERGs are labelled in red (upregulated) and blue (downregulated). (B–C) The functional enrichment of upregulated (B) and downregulated (C) FERGs. *p < 0.05, **p < 0.001 and ***p < 0.0001 as calculated by (A) the Wilcoxon rank‐sum test

Tumour somatic mutation landscape of ferroptosis‐related genes

We also investigated the somatic mutation of FERGs using the TCGA mutation data (as mentioned in Materials and Methods). We found that the FERGs TP53 (69%), KEAP1 (25%) and ZEB1 (9%) had the highest mutation frequency (Figure 2A). Moreover, the missense mutation was most common across these FERGs (Figure 2B). Furthermore, comparing the expression levels of FERGs, we observed that there were general differences in expression levels between mutated and wild tumours (Figure 2C). It is interesting to note that the somatic mutations of TP53, KEAP1 and ZEB1 were also important factors affecting their expression levels. These results revealed that tumour somatic mutation was generally distributed among the FERGs, and the mutated FERGs were capable of affecting their expression, which might be associated with LUAD development.

FIGURE 2

Mutation‐perturbed ferroptosis‐related genes in lung adenocarcinoma. (A) Mutation landscape of FERGs. The top 20 genes in mutation levels are demonstrated. (B) Bar plot depicting the SNV class. (C) Violin plot depicting the expression levels of FERGs among mutated, wild tumour and normal samples. *p < 0.05, **p < 0.001 and ***p < 0.0001 as calculated by (C) the Wilcoxon rank‐sum test

Ferroptosis‐related subtypes revealed immune heterogeneity and prognostic values

To further assess the role of FERGs in LUAD, we systemically identified ferroptosis‐related tumour subtypes based on the expression levels of differentially expressed FERGs using the consensus clustering analysis (as mentioned in Materials and Methods). The analysis resulted in three ferroptosis‐related tumour subtypes (F1, F2 and F3) (Figure 3A,B). We discovered that there were distinct expression patterns of FERGs among the ferroptosis‐related subtypes, implying the differences among them at the molecular and clinical levels (Figure 3C). Furthermore, the aneuploidy score (Figure 3D) and mutation count (Figure 3E) in the F3 subtypes were significantly lower than those in F1 and F2 indicating that F3 had higher genome stability. Interestingly, we also observed that the tumour immune score of F3 was significantly greater (Figure 3F) and the tumour purity of F3 was significantly lower than those of F1 and F2 (Figure 3G), respectively, implying that F3 was a tumour subtype with high immune activity. Therefore, we compared six immune cell infiltration levels among these subtypes. B cells, CD4‐T cells, CD8‐T cells, macrophages, neutrophils and dendritic cells had higher infiltration abundances in F3 than in F1 and F2, thus highlighting that the F3 subtype might have a greater ability to regulate and kill tumour cells (Figure 3H). Eventually, we also found that patients in F3 had better overall survival outcomes than those of F1 and F2 (Figure 3I). These results suggested that ferroptosis‐related subtypes played an important role in revealing intertumoral immune heterogeneity and had potential prognostic value in LUAD.

FIGURE 3

Ferroptosis‐related subtypes revealed immune heterogeneity and distinct survival outcomes. (A) Heat map of consensus scores when k is 3. (B) Relative change in area under the CDF curve when k ranges from 2 to 10. (C) Heatmap showing the expression levels of FERGs in different ferroptosis‐related subtypes. (D–H) Violin plots showing the (D) aneuploidy score, (E) mutation count, (F) immune score, (G) ESTIMATE score and (H) immune cell infiltration abundance among ferroptosis‐related subtypes. (I) Kaplan–Meier curves of patients stratified by the different ferroptosis‐related subtypes with log‐rank test. *p < 0.05, **p < 0.001 and ***p < 0.0001 as calculated by (D–H) the Wilcoxon rank‐sum test

Ferroptosis‐related subtypes prompt treatment benefits of radiation therapy

We investigated the potential clinical treatment benefits of different ferroptosis‐related subtypes in patients who received radiation therapy in the TCGA cohort. We found that patients in the F3 subtype had better overall (Figure 4A) and disease‐free (Figure 4B) survival outcomes than those of F1 and F2 subtypes. The results indicated that ferroptosis‐related subtypes could prompt the clinical treatment benefits of radiation therapy.

FIGURE 4

Ferroptosis‐related subtypes prompt the clinical treatment benefits of radiation therapy. Kaplan–Meier curves of patients who received radiotherapy stratified by the different ferroptosis‐related subtypes with log‐rank test p‐values provided based on the (A) overall and (B) disease‐specific survival

Exploration of ferroptosis‐related genes in clinical benefits of anti–PD‐1/PD‐L1 immunotherapy

The abovementioned results revealed that ferroptosis‐related subtypes have the potential for indicating intertumoral immune heterogeneity and clinical treatment benefits. Therefore, we investigated the role of FERGs in anti–PD‐1/PD‐L1 immunotherapy. First, we characterised the correlation between the expressions of each FERG and PD‐1/PD‐L1, which suggested that they were generally associated (Figure 5A). We discovered that glucose‐6‐phosphate dehydrogenase (G6PD) was most related to PD‐1 and PD‐L1, implying that it might be a good effector in immunotherapy. Furthermore, we established an RPART classification model to predict the durable or non‐durable clinical treatment benefits in an immunotherapy dataset (as mentioned in Materials and Methods). We identified a FER signature containing CD44, G6PD and ZEB1 that acted as the immunotherapy benefit indicator. To further test the efficiency of the FER signature, we constructed an FES score for extracting the sample features based on the expression levels of the signature in another independent immunotherapy dataset (as mentioned in Materials and Methods). The analysis indicated that the FES score had a promising potential for predicting the immunotherapy benefits of patients (AUC = 0.697) (Figure 5C). These results revealed that FER signature was a good biomarker for prompting the immunotherapy benefits of patients.

FIGURE 5

Ferroptosis‐related subtypes prompt the clinical treatment benefits of immunotherapy. (A) Bubble plot showing the Spearman's correlation between the expression of PD‐1 or PD‐L1 and FERGs in the TCGA LUAD cohort (B) RPART decision tree for predicting immunotherapy benefits (C) ROC curves showing the decision efficiency based on the FES score, CD44, G6PD and ZEB1

Immunotherapy‐related ferroptosis signature was involved in complex biological processes

For a better understanding of the effect of FER signature on immunotherapy benefits, we constructed a protein–protein interaction (PPI) network for each FERG in the signature (as mentioned in Materials and Methods). In these PPI networks, we observed that the FER signature could interact with some well‐known cancer‐related genes, such as G6PD–TP53 (Figure 6A), CD44–EGFR, CD44–ERBB2 (Figure 6B) and ZEB1–MYC (Figure 6C). The results revealed that the FER signature had an important role in cancer development. Furthermore, we systemically annotated the biological functions of each FERG in FER signature based on their corresponding PPI interaction network. The results suggested that CD44 and ZEB1 were mainly enriched in the cancer‐related (hsa05200, ko05206) and cell cycle‐related (GO:0030155, GO:0045596 and GO:2000147) terms. It is interesting to note that ZEB1 specifically engaged in some signalling pathways, such as the Wnt signalling pathway (WP399), interleukin signalling pathway (R‐HSA‐449147) and intracellular receptor signalling pathway (GO:0030522). Alternatively, the biological functions of G6PD are presented in multiple metabolic processes, including the monocarboxylic acid metabolic process (GO:0032787), carbon metabolism (hsa01200) and glucose 6‐phosphate metabolic process (GO:0051156). These results further highlighted the importance of the FER signature in LUAD development.

FIGURE 6

Functional enrichment of the FER signature. (A–C) The PPI network for (A) G6PD, (B) CD44 and (C) ZEB1 (D) heat map showing the biological functions of each FER signature involved

DISCUSSION

Emerging evidence has revealed that ferroptosis is closely associated with lung cancer. , Herein, we reported that the expression of most FERGs was dysregulated in LUAD. These FERGs were involved in many metabolic processes that might affect cancer development. Furthermore, we constructed the somatic mutation landscape of FERGs and found that somatic mutations influenced the expression levels of FERGs. TP53, KEAP1 and ZEB1 had the highest mutation frequency, and their expression levels were also impacted between mutated and wild samples. The results indicated that the genetic alterations of FERGs were essential signals in LUAD progression.

To better understand the role of altered FERGs in LUAD, we identified ferroptosis‐related subtypes based on the consensus clustering analysis. The analysis resulted in three robust subtypes: F1, F2 and F3. Comparing the molecular features, we found that patients in the F3 subtype usually had lower mutation counts and aneuploidy scores than those of the other two subtypes, implying that F3 had higher genome stability. Notably, we also observed that the F3 subtype showed lower tumour purity and higher immunoactivity. B cells, CD4‐T cells, CD8‐T cells, macrophages, neutrophils and dendritic cells exhibited higher infiltration abundances in F3 than those in F1 and F2. These results suggested that F3 was an immune‐related subtype. In addition, we found that patients in F3 typically experienced longer overall survival and had a better prognosis. Moreover, we observed that the F2 subtype exhibited the highest tumour purity, lowest tumour immunoactivity and worst overall survival outcomes among the three subtypes. In addition, patients with radiotherapy in the F2 subtype also had the poorest prognostic benefits. The poor prognosis outcomes might be associated with the immune‐depleted tumour immune microenvironment, which implied the crosstalk between ferroptosis and tumour immune activity. These results also highlighted that the ferroptosis‐related subtypes could contribute to the clinical performance and intertumoral immune heterogeneity of patients.

First, to enhance the application of FERGs in the clinical management of patients, we explored the relationship between radiation therapy and ferroptosis‐related subtypes. Consistent with the abovementioned results, the F3 subtype showed better survival outcomes in both overall and disease‐free survival in patients who received radiotherapy. In addition, we assessed the role of FERGs in anti–PD‐1/PD‐L1 immunotherapy. By constructing an RPART decision model, we identified a FER signature for predicting the clinical benefits of immunotherapy. The FER signature consisted of CD44, G6PD and ZEB1. The robustness of the FER signature was verified in an independent immunotherapy dataset.

In conclusion, our study provided a novel insight for understanding ferroptosis in LUAD. FRGs were also associated with immune heterogeneity and had a promising potential for predicting immunotherapy benefits.

CONFLICT OF INTEREST

The authors declared no potential conflicts of interest in terms of the research, authorship and/or publication of this article.

AUTHOR CONTRIBUTIONS

Y.C. and M.Z. designed the study and performed analysis. X.Z. collected data. T.L. proofread the manuscript. Y.C. supervised the research. All authors read and approved the final manuscript.

Figure S1

Click here for additional data file.