Prognostic and immune implications of a novel ferroptosis-related ten-gene signature in lung adenocarcinoma.

BACKGROUND: Lung cancer has been the focus of attention for many researchers in recent years due to its leading contribution to cancer-related death worldwide, with lung adenocarcinoma (LUAD) being the most common histological type. Ferroptosis, a novel iron-dependent form of regulated cell death, can be induced by sorafenib. Emerging evidence shows that triggering ferroptosis has potential as a cancer therapy. This work aimed to build a ferroptosis-related gene signature for predicting the outcome of LUAD.
METHODS: The TCGA-LUAD dataset was set as the training cohort, and the GSE72094 and GSE68465 datasets were set as the validation cohorts. Sixty-two ferroptosis-related genes were retrieved from the literature. A univariate Cox regression model was constructed for the training cohort to preliminarily screen for potential prognostic ferroptosis-related genes. A gene signature was generated from a LASSO Cox regression model and assessed with the training and validation cohorts through Kaplan-Meier, Cox, and ROC analyses. In addition, the correlation between the risk score and autophagy-related genes was determined by the Pearson test. Finally, GSEA and immune infiltrating analyses were performed to better study the functional annotation of the signature and the role of each kind of immune cell.
RESULTS: A ten-gene signature was constructed from the training cohort and validated in three cohorts by Kaplan-Meier and Cox regression analyses, revealing its independent prognostic value in LUAD. Moreover, a ROC analysis conducted with all cohort data confirmed the predictive ability of the ten-gene signature for LUAD prognosis. A total of 62.85% (308/490) of autophagy-related genes were found to be significantly correlated with risk scores. GSEA detailed the exact pathways related to the gene signature, and immune-infiltrating analyses identified crucial roles for resting mast cells and resting dendritic cells in the prognosis of LUAD.
CONCLUSIONS: We identified a novel ferroptosis-related ten-gene signature (PHKG2, PGD, PEBP1, NCOA4, GLS2, CISD1, ATP5G3, ALOX15, ALOX12B, and ACSL3) that can accurately predict LUAD prognosis and is closely linked to resting mast cells and resting dendritic cells. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer is the most common malignant tumor in humans and causes more cancer-related mortality worldwide than any other disease (1). Lung cancer mainly consists of two subtypes, non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC). NSCLC accounts for almost 80% of lung cancer cases and comprises two major types, lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (2). LUAD is the predominant histology, and the incidence rate is still increasing (1). Several treatment advances have been made in recent years, in particular, the advancement of targeted therapy and the emergence of immunotherapy (3,4). However, these two methods can only benefit a limited number of subtypes, and the overall survival rate of patients with LUAD is still very low (5). Therefore, it is important to continue efforts to find specific prognostic methods for the survival of specific patients so that the most suitable treatment and management programs can be designed for different subgroups of patients with LUAD (6).

Ferroptosis is a type of programmed cell death dependent on iron and characterized by the accumulation of lipid peroxides. It is genetically and biochemically distinct from other forms of regulated cell death, such as apoptosis (7). The induction of ferroptosis has become a promising treatment option for cancer cell death over the past few years, especially for malignant tumors resistant to traditional therapies (8-10). In addition to ferroptosis-inducing agents, many genes have also been identified as regulators or hypertrophy markers (10-14).

Recently, several studies have mined the prognostic gene signature related to ferroptosis in tumors from public databases (11,15). Liu confirmed that the ferroptosis-related nineteen-gene signature could predict survival in patients with glioma (11). Liang et al. discovered a novel ferroptosis-related prognostic gene signature for hepatocellular carcinoma (15). Several studies have identified prognostic gene signatures of LUAD involving glycolysis (16), metabolism (17), or immunity (18). These studies did not explain the role of ferroptosis in their models. To date, no studies have determined whether a ferroptosis-related prognostic gene signature can predict LUAD outcomes. To fill this gap and discover potential diagnostic methods potentially targeting ferroptosis, this study used ferroptosis-related genes identified in previous studies and public databases to identify a ferroptosis-related prognostic gene signature, which was also validated in independent datasets. Gene Set Enrichment Analysis (GSEA) and immune infiltration analysis were used to better understand the gene signature. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/atm-20-7936).

Methods

Cohorts and ferroptosis-related genes

The TCGA-LUAD dataset was chosen as the training cohort and downloaded from the GDC Xena Hub (https://gdc.xenahubs.net). Five hundred patients with LUAD in the TCGA-LUAD dataset who had both gene expression and survival data available were included in our research. Two independent cohorts, GSE72094 and GSE68465, which consisted of 442 and 443 LUAD cases, respectively, were selected for the study as validation cohorts to validate the developed prognostic model. The detailed survival data of the cohorts are shown in https://cdn.amegroups.cn/static/public/atm-20-7936-1.pdf. A comprehensive list consisting of a total of 62 ferroptosis-related genes was retrieved from the literature (10-14) and is provided in https://cdn.amegroups.cn/static/public/atm-20-7936-2.pdf. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification and validation of the prognostic ferroptosis-related gene signature

A univariate Cox regression analysis was performed based on overall survival in the training cohort to identify potential prognostic genes (P value <0.05). The potential prognostic genes that overlapped with ferroptosis-related genes were identified as potential prognostic ferroptosis-related genes, which were then entered into an overall survival-based LASSO Cox regression model in the training cohort. The “glmnet” R package with a setting of 10-fold cross-validations was applied for LASSO analysis (19-22). R software generated a list of prognostic genes with coefficients according to the selected optimal lambda value. Each patient’s risk score can be obtained based on the expression level of each prognostic gene and its corresponding coefficient, as shown in the following formula:

where n, Expi, and βi represent the number of hub genes, gene expression level, and regression coefficient value, respectively.

Kaplan-Meier analysis, univariate and multivariate Cox analyses and ROC analyses were conducted for the training and validation cohorts to assess the prognostic value of the gene signature risk score.

Relationship between gene signature and autophagy

Autophagy is a conserved intracellular degradation system that plays a dual role in cell death; thus, therapies targeting autophagy in cancer are somewhat controversial (23). Accumulating studies have revealed crosstalk between autophagy and ferroptosis at the molecular level (23). We explored the relationship between autophagy and our gene signature and first identified 232 autophagy-associated genes from the Human Autophagy Database (HADb; http://www.autophagylu/index.html), which consists of an exhaustive, up-to-date list of human autophagy-related genes (24). Another 363 autophagy-related genes were retrieved from the Molecular Signatures Database (version 7.1, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). A list of 490 autophagy-related genes was obtained by merging these lists (https://cdn.amegroups.cn/static/public/atm-20-7936-3.pdf). The Pearson correlation coefficient was used to evaluate the correlation between autophagy and the risk score. A P value<0.05 was considered statistically significant.

GSEA

GSEA was performed based on Hallmark gene set collections (v7.1, https://www.gsea-msigdb.org/gsea/downloads.jsp) to uncover the pathways and functions of the genes that were differentially expressed in the high-risk and low-risk groups. The significant gene set selection criteria were as follows: | NES | >1, NOM p-val <0.05, and FDR q-val <0.25.

Relationship between the gene signature and 22 tumor-infiltrating immune cells (TICs)

The relative proportion of 22 TICs in the training group was downloaded from the National Cancer Institute web portal (https://gdc.cancer.gov/about-data/publications/panimmune); it had been calculated using the CIBERSORT algorithm (25,26). Three-hundred-ninety LUAD cases were eligible to be entered into the subsequent analysis after quality filtering (P value <0.05), and the Pearson coefficient was used to evaluate the correlations between the 22 TICs and LUAD. The Spearman coefficient and Wilcoxon rank-sum tests were applied to identify the relationship between the proportions of the 22 TICs and the risk score. Additionally, univariate Cox and Kaplan-Meier analyses were used to assess the prognostic value of the 22 TICs. The potential candidate TIC that plays a role in the gene signature’s prognostic capacity was found by combining the results of these analyses.

Statistical analysis

Univariate and multivariate Cox proportional hazard regression analyses were conducted using the “survival” R package. Kaplan-Meier analysis was performed using the “survival” and “survminer” R packages. ROC analysis was applied using the “survivalROC” R package. The area under the ROC curve (AUC) served as a marker of prognostic accuracy. A P value <0.05 was considered statistically significant.

Results

Cohorts’ characteristics

This research flowchart is shown in . Five hundred LUAD cases from the TCGA-LUAD cohort were taken as the training cohort. The GSE72094 and GSE68465 cohorts from the GEO database, consisting of 442 and 443 LUAD cases, respectively, were selected as the validation cohorts. summarizes the detailed clinical characteristics of the cohorts in this study.

Figure 1

Flow chart of the study. LASSO, least absolute shrinkage and selection operator Cox regression model; ROC, receiver operating characteristic; LUAD, lung adenocarcinoma; TICs, tumor-infiltrating immune cells.

Table 1

Clinical characteristics of patients involved in the study

Characteristics	Training cohort (TCGA-LUAD, 500 cases)	Validation cohort (GSE72094, 442 cases)	Validation cohort (GSE68465, 443 cases)
Age
<65	219 (43.8%)	115 (26.02%)	214 (48.31%)
≥65	271 (54.2%)	306 (69.23%)	229 (51.69%)
Unknown	10 (2%)	21 (4.75%)	0
Gender
Female	270 (54%)	240 (54.3%)	220 (49.66%)
Male	230 (46%)	202 (45.7%)	223 (50.34%)
T classification
T1	167 (33.4%)	NA	150 (33.86%)
T2	267 (53.4%)	NA	251 (56.66%)
T3	45 (9%)	NA	28 (6.32%)
T4	18 (3.6%)	NA	12 (2.71%)
Unknown	3 (0.6%)	NA	2 (0.45%)
N classification
N0	324 (64.8%)	NA	299 (67.49%)
N1	94 (18.8%)	NA	88 (19.86%)
N2	69 (13.8%)	NA	53 (11.96%)
N3	2 (0.4%)	NA	0
Unknown	11 (2.2%)	NA	3 (0.68%)
M classification
M0	332 (66.4%)	NA	NA
M1	24 (4.8%)	NA	NA
Unknown	144 (28.8%)	NA	NA
Tumor stage
Stage I	268 (53.6%)	265 (59.95%)	NA
Stage II	119 (23.8%)	69 (15.61%)	NA
Stage III	80 (16%)	63 (14.25%)	NA
Stage IV	25 (5%)	17 (3.85%)	NA
Unknown	8 (1.6%)	28 (6.33%)	NA
Race
White	386 (77.2%)	399 (90.27%)	295 (66.59%)
Black or African American	52 (10.4%)	13 (2.94%)	12 (2.71%)
American Indian or Alaska Native	1 (0.2%)	0	1 (0.23%)
Asian	7 (1.4%)	3 (0.68%)	6 (1.35%)
Unknown	54 (10.8%)	27 (6.11%)	129 (29.12%)
Ethnicity
Hispanic or Latino	7 (1.4%)	10 (2.26%)	NA
Not Hispanic or Latino	381 (76.2%)	402 (90.95%)	NA
Unknown	112 (22.4%)	30 (6.79%)	NA
Tobacco smoking history
Ever	415 (83%)	335 (75.79%)	300 (67.72%)
Never	71 (14.2%)	33 (7.47%)	49 (11.06%)
Unknown	14 (2.8%)	74 (16.74%)	94 (21.22%)
Number pack-years smoked
<30	115 (23%)	NA	NA
≥30	227 (45.4%)	NA	NA
Unknown	158 (31.6%)	NA	NA
Radiation therapy
Yes	58 (11.6%)	NA	65 (14.67%)
No	361 (72.2%)	NA	364 (82.17%)
Unknown	81 (16.2%)	NA	14 (3.16%)
Chemotherapy
Yes	NA	NA	89 (20.09%)
No	NA	NA	341 (76.98%)
Unknown	NA	NA	13 (2.93%)
Additional pharmaceutical therapy
Yes	61 (12.2%)	NA	NA
No	68 (13.6%)	NA	NA
Unknown	371 (74.2%)	NA	NA
Additional radiation therapy
Yes	61 (12.2%)	NA	NA
No	71 (14.2%)	NA	NA
Unknown	368 (73.6%)	NA	NA
EGFR mutation
Yes	80 (16%)	47 (10.63%)	NA
No	191 (38.2%)	395 (89.37%)	NA
Unknown	229 (45.8%)	0	NA
KRAS mutation
Yes	23 (4.6%)	154 (34.84%)	NA
No	39 (7.8%)	288 (65.16%)	NA
Unknown	438 (87.6%)	0	NA
STK11 mutation
Yes	NA	68 (15.38%)	NA
No	NA	374 (84.62%)	NA
TP53 mutation
Yes	NA	111 (25.11%)	NA
No	NA	331 (74.89%)	NA
EML4-ALK mutation
Yes	34 (6.8%)	NA	NA
No	207 (41.4%)	NA	NA
Unknown	259 (51.8%)	NA	NA
Tumor intermediate dimension
<1	287 (57.4%)	NA	NA
≥1	88 (17.6%)	NA	NA
Unknown	125 (25%)	NA	NA
Location in lung parenchyma
Central lung	62 (12.4%)	NA	NA
Peripheral lung	118 (23.6%)	NA	NA
Unknown	320 (64%)	NA	NA
Person neoplasm cancer status
With tumor	134 (26.8%)	NA	NA
Tumor free	243 (48.6%)	NA	NA
Unknown	123 (24.6%)	NA	NA
Vital status
Alive	318 (63.6%)	298 (67.42%)	207 (46.73%)
Dead	182 (36.4%)	122 (27.6%)	236 (53.27%)
Unknown	0	22 (4.98%)	NA

Identification of a prognostic ferroptosis-related gene signature from the training cohort

A total of 7,596 genes were identified from the overall survival-based univariate Cox regression analysis in the training cohort (https://cdn.amegroups.cn/static/public/atm-20-7936-4.pdf). The 7,596 potential prognostic genes and 62 ferroptosis-related genes were intersected to obtain a list consisting of 14 ferroptosis-related potential prognostic genes (). The 14 ferroptosis-related potential prognostic genes were then subjected to an overall survival-based LASSO Cox regression model (). The model achieved the best performance when 10 genes were included (). The regression coefficient of each gene was calculated and is shown in .

Table 2

14 ferroptosis-related potential prognostic genes generated from the training cohort

Gene symbol	Description	HR	95% confidence interval		P value
			Lower	Upper
PEBP1	Phosphatidylethanolamine-Binding Protein 1	0.603658936	0.464865289	0.78389185	0.000152746
ACSL3	Acyl-CoA Synthetase Long-Chain Family Member 3	1.435108616	1.157210163	1.779743045	0.001003251
GLS2	Glutaminase 2	0.356469987	0.176890046	0.71836067	0.003911814
GCLC	Glutamate-Cysteine Ligase Catalytic Subunit	1.135847381	1.034367787	1.24728292	0.007639188
SLC7A11	Solute Carrier Family 7 Member 11	1.137868621	1.023951703	1.264459052	0.016406631
PHKG2	Phosphorylase Kinase Catalytic Subunit Gamma 2	0.69374619	0.511654438	0.940642238	0.018577862
PGD	Phosphogluconate Dehydrogenase	1.177966423	1.02132216	1.358635842	0.02446422
CISD1	CDGSH Iron Sulfur Domain 1	1.403301444	1.041497232	1.890792296	0.025931782
ALOX12B	Arachidonate 12-Lipoxygenase, 12R Type	1.528527805	1.03748207	2.251988076	0.031866869
DPP4	Dipeptidyl Peptidase 4	0.911077238	0.836744384	0.992013511	0.031983077
ALOX15	Arachidonate 15-Lipoxygenase	0.857033797	0.740142144	0.992386306	0.039196308
ATP5G3	ATP Synthase Membrane Subunit C Locus 3	1.311553769	1.01163069	1.700396506	0.040633369
ABCC1	ATP Binding Cassette Subfamily C Member 1	1.23437377	1.007318781	1.512608156	0.042328557
NCOA4	Nuclear Receptor Coactivator 4	0.776056577	0.603024235	0.998738982	0.048865024

Figure 2

Construction of the prognostic gene signature using LASSO regression analysis. (A) LASSO coefficient profiles of 14 ferroptosis-related potential prognostic genes. Each curve corresponds to a gene. (B) Tenfold cross-validation for tuning parameter selection in the LASSO model. The partial likelihood deviance is plotted against log (Lambda), where Lambda is the tuning parameter. Partial likelihood deviance values are shown, with error bars representing SE. The dotted vertical lines are drawn at the optimal values by minimum criteria and 1-SE criteria. LASSO: the least absolute shrinkage and selection operator Cox regression model.

Table 3

10 ferroptosis-related prognostic genes obtained from LASSO Cox regression model

Gene symbol	Description	Risk coefficient
PHKG2	Phosphorylase Kinase Catalytic Subunit Gamma 2	−0.376487238
PGD	Phosphogluconate Dehydrogenase	0.027750809
PEBP1	Phosphatidylethanolamine Binding Protein 1	−0.369553219
NCOA4	Nuclear Receptor Coactivator 4	−0.498545244
GLS2	Glutaminase 2	−0.399525822
CISD1	CDGSH Iron Sulfur Domain 1	0.313326653
ATP5G3	ATP Synthase Membrane Subunit C Locus 3	0.019018653
ALOX15	Arachidonate 15-Lipoxygenase	−0.049973524
ALOX12B	Arachidonate 12-Lipoxygenase, 12R Type	0.323899465
ACSL3	Acyl-CoA Synthetase Long Chain Family Member 3	0.246495026

Prognostic ability of the ten-gene signature

Each patient’s risk score was a linear combination of each ten-gene signature expression and its risk coefficient. LUAD cases were assigned to high- and low-risk groups based on the median risk score. The distribution plots of risk scores, outcome status, and expression profiles of the gene signature in the three cohorts are shown in . As demonstrated in , more events occurred in the high-risk groups than in their corresponding low-risk groups. Additionally, the patients in the high-risk groups had a shorter survival time than those in the respective low-risk groups. The heat maps show that in high-risk cases, ACSL3, PGD, ATP5G3, CISD1, and ALOX12B were overexpressed, while the expression of NCOA4, GLS2, ALOX15, PEBP1, and PHKG2 was downregulated. In addition, we checked the five-year survival prediction capacity of the ten-gene signature () and found that more events and less survival time existed in the high-risk groups than in the low-risk groups, which were consistent with those shown in .

Figure 3

The overall performance of the ten-gene signature in the training (A, B, E, and F) and validation (C, D, G, and H) cohorts. The upper and middle parts of each plot indicate the distribution of the risk score and patient survival time, respectively. The bottom part of each graph shows a heatmap of ten gene expression profiles.

Kaplan-Meier curves showed that the high-risk patients had poor survival rates in the TCGA-LUAD (overall survival: P value <0.0001, ; progression-free survival: P value =0.00064, ), GSE72094 (P value <0.0001, ), and GSE68465 (P value =0.00043, ) cohorts compared to the survival of specific low-risk patients among the cohorts. The Kaplan-Meier curves of five-year survival showed the same pattern: high-risk score groups had significantly more unfavorable outcomes than their corresponding low-risk groups ().

Figure 4

Kaplan-Meier curves of the ten-gene signature risk score in the training (A, B, E, and F) and validation (C, D, G, and H) cohorts. The middle part of each graph indicates the number of patients at risk. The two-sided log-rank test measured the differences between the high- and low-risk groups with a P value <0.05.

Univariate and multivariate Cox analyses were applied in the three cohorts based on overall or progression-free survival using available covariables, including ten-gene signature risk score, sex, age, race, tumor stage, tobacco smoking history, radiation therapy, KRAS mutation, and EML4-ALK mutation to validate the prognostic capacity and the independence of the ten-gene signature among other clinicopathological characteristics (). In the overall survival-based Cox regression model of the training cohort, both univariate Cox and multivariate Cox models suggested that the ten-gene signature was a powerful player (HR =3.39, 95% CI, 2.49–4.62, P value =8.78E-15, and HR =7.32, 95% CI, 2.88–18.6, P value =2.94E-05, respectively). Consistent with that in the training cohort, in both univariate and multivariate analyses, the ten-gene signature showed excellent ability in the other two independent validation cohorts to predict overall survival (P value ≤1.22E-04). Additionally, we utilized progression-free survival data in the training cohort to perform the Cox analysis. We found that the ten-gene signature had the ability to predict the outcomes in both the univariate and multivariate models (P value ≤7.35E-03). The above evidence shows that the ten-gene signature is a powerful independent prognostic factor.

Table 4

Univariate analysis and multivariate analysis of the correlation of gene signature risk score with outcomes among lung adenocarcinoma patients in three cohorts

Variable	Univariate Cox analysis					Multivariate Cox analysis
	coef	HR (95% CI)	z	P value		coef	HR (95% CI)	z	P value
TCGA-LUAD (overall survival)^
Gender (male vs. female)	0.048	1.05 (0.784–1.4)	0.322	7.47E-01
Age	0.00816	1.01 (0.993–1.02)	1.05	2.94E-01
Tumor stage II (vs. stage I)	0.905	2.47 (1.72–3.56)	4.87	1.09E-06*		1.32	3.75 (1.44–9.8)	2.7	6.93E-03
Tumor stage III (vs. stage I)	1.25	3.49 (2.38–5.12)	6.4	1.51E-10*		0.455	1.58 (0.445–5.59)	0.705	4.81E-01
Tumor stage IV (vs. stage I)	1.34	3.82 (2.2–6.62)	4.76	1.92E-06*		−0.107	0.899 (0.255–3.17)	−0.166	8.68E-01
Tobacco smoking history (ever vs. never)	−0.127	0.88 (0.583–1.33)	−0.605	5.45E-01
Number pack years smoked	0.00327	1 (0.996–1.01)	0.938	3.48E-01
Race (White vs. non-White)	0.36	1.43 (0.875–2.34)	1.43	1.53E-01
Ethnicity (Hispanic or Latino vs. Non-Hispanic or Latino)	0.387	1.47 (0.466–4.65)	0.658	5.10E-01
Tumor intermediate dimension	0.441	1.55 (0.922–2.62)	1.65	9.80E-02
Radiation therapy (yes vs. no)	0.763	2.14 (1.44–3.19)	3.76	1.73E-04*		0.446	1.56 (0.415–5.88)	0.66	5.09E-01
Additional pharmaceutical therapy (yes vs. no)	−0.502	0.605 (0.382–0.959)	−2.14	3.26E-02*		−1.03	0.356 (0.158–0.8)	−2.5	1.25E-02*
Additional radiation therapy (yes vs. no)	−0.0153	0.985 (0.625–1.55)	−0.0661	9.47E-01
Location in lung parenchyma (central lung vs. peripheral lung)	0.0908	1.09 (0.684–1.75)	0.378	7.06E-01
KRAS mutation (yes vs. no)	0.492	1.63 (0.672–3.98)	1.08	2.78E-01
EGFR mutation (yes vs. no)	0.268	1.31 (0.828–2.06)	1.15	2.50E-01
EML4-ALK mutation (yes vs. no)	0.592	1.81 (1.01–3.24)	1.98	4.73E-02*		1.18	3.25 (1.22–8.63)	2.36	1.84E-02*
Risk score	1.22	3.39 (2.49–4.62)	7.76	8.78E-15*		1.99	7.32 (2.88–18.6)	4.18	2.94E-05*
TCGA-LUAD (progression-free survival)#
Gender (male vs. female)	0.07	1.07 (0.815–1.41)	0.499	6.17E-01
Age	−0.000998	0.999 (0.985–1.01)	−0.141	8.88E-01
Tumor stage II (vs. stage I)	0.798	2.22 (1.61–3.07)	4.85	1.22E-06*		1.2	3.33 (1.92–5.78)	4.28	1.88E-05*
Tumor stage III (vs. stage I)	0.709	2.03 (1.38–2.99)	3.61	3.10E-04*		0.95	2.59 (1.31–5.11)	2.73	6.25E-03*
Tumor stage IV (vs. stage I)	0.849	2.34 (1.32–4.13)	2.93	3.40E-03*		0.828	2.29 (0.979–5.35)	1.91	5.60E-02
Tobacco smoking history (ever vs. never)	−0.0557	0.946 (0.637–1.4)	−0.276	7.82E-01
Number pack years smoked	−0.000552	0.999 (0.993–1.01)	−0.165	8.69E-01
Race (White vs. non-White)	0.0557	1.06 (0.703–1.59)	0.268	7.89E-01
Ethnicity (Hispanic or Latino vs. non-Hispanic or Latino)	0.291	1.34 (0.495–3.62)	0.574	5.66E-01
Tumor intermediate dimension	0.328	1.39 (0.84–2.29)	1.28	2.01E-01
Radiation therapy (yes vs. no)	0.707	2.03 (1.42–2.89)	3.92	9.01E-05*		0.432	1.54 (0.829–2.86)	1.37	1.72E-01
Additional pharmaceutical therapy (yes vs. no)	0.27	1.31 (0.914–1.88)	1.47	1.41E-01
Additional radiation therapy (yes vs. no)	0.16	1.17 (0.825–1.67)	0.89	3.73E-01
Location in lung parenchyma (central lung vs. peripheral lung)	−0.0709	0.932 (0.584–1.49)	−0.298	7.66E-01
KRAS mutation (yes vs. no)	−0.275	0.759 (0.316–1.82)	−0.616	5.38E-01
EGFR mutation (yes vs. no)	0.483	1.62 (1.06–2.49)	2.21	2.71E-02*		0.277	1.32 (0.812–2.14)	1.12	2.64E-01
EML4-ALK mutation (yes vs. no)	0.509	1.66 (0.965–2.87)	1.83	6.72E-02
Risk score	0.853	2.35 (1.74–3.16)	5.63	1.80E-08*		0.727	2.07 (1.22–3.52)	2.68	7.35E-03*
GSE72094 (overall survival)&
Gender (male vs. female)	0.44	1.55 (1.07–2.25)	2.33	1.98E-02*		0.464	1.59 (1.09–2.33)	2.39	1.67E-02*
Age	0.00696	1.01 (0.988–1.03)	0.702	4.83E-01
Race (White vs. non-White)	−0.0694	0.933 (0.38–2.29)	−0.151	8.80E-01
Ethnicity (Hispanic or Latino vs. non-Hispanic or Latino)	−0.642	0.526 (0.0733–3.78)	−0.638	5.24E-01
Tobacco smoking history (ever vs. never)	0.314	1.37 (0.597–3.14)	0.741	4.59E-01
Tumor stage II (vs. stage I)	0.758	2.13 (1.32–3.44)	3.11	1.85E-03*		0.708	2.03 (1.25–3.28)	2.88	3.92E-03*
Tumor stage III (vs. stage I)	1.13	3.09 (1.93–4.97)	4.67	3.00E-06*		1.17	3.22 (1.99–5.22)	4.76	1.91E-06*
Tumor stage IV (vs. stage I)	1.21	3.35 (1.59–7.06)	3.18	1.48E-03*		1.14	3.13 (1.48–6.62)	2.99	2.80E-03*
KRAS mutation (yes vs. no)	0.376	1.46 (1–2.12)	1.97	4.92E-02*		0.157	1.17 (0.796–1.72)	0.799	4.24E-01
EGFR mutation (yes vs. no)	−1.34	0.262 (0.0965–0.71)	−2.63	8.49E-03*		−0.926	0.396 (0.143–1.1)	−1.78	7.52E-02
STK11 mutation (yes vs. no)	−0.0393	0.961 (0.58–1.59)	−0.153	8.79E-01
TP53 mutation (yes vs. no)	0.211	1.23 (0.82–1.86)	1.01	3.13E-01
Risk score	0.613	1.85 (1.47–2.31)	5.33	1.00E-07*		0.508	1.66 (1.31–2.11)	4.14	3.54E-05*
GSE68465 (overall survival)$
Gender (male vs. female)	0.33	1.39 (1.07–1.8)	2.5	1.25E-02*		0.324	1.38 (1.05–1.83)	2.27	2.32E-02*
Age	0.0266	1.03 (1.01–1.04)	3.97	7.16E-05*		0.0357	1.04 (1.02–1.05)	4.93	8.10E-07*
Race (White vs. non-White)	0.32	1.38 (0.702–2.71)	0.931	3.52E-01
Chemotherapy (yes vs. no)	0.488	1.63 (1.22–2.18)	3.29	9.99E-04*		0.243	1.28 (0.861–1.89)	1.21	2.25E-01
Radiation therapy (yes vs. no)	0.698	2.01 (1.47–2.75)	4.37	1.26E-05*		0.26	1.3 (0.838–2)	1.17	2.43E-01
T classification (T2 vs. T1)	0.381	1.46 (1.09–1.97)	2.52	1.18E-02*		0.126	1.13 (0.829–1.55)	0.786	4.32E-01
T classification (T3 vs. T1)	1.2	3.31 (2.06–5.34)	4.92	8.78E-07*		0.76	2.14 (1.29–3.55)	2.94	3.27E-03*
T classification (T4 vs. T1)	1.45	4.26 (2.17–8.34)	4.22	2.41E-05*		1.06	2.89 (1.36–6.13)	2.77	5.59E-03*
N classification (N1 vs. N0)	0.821	2.27 (1.67–3.09)	5.23	1.74E-07*		0.775	2.17 (1.57–3)	4.69	2.74E-06*
N classification (N2 vs. N0)	1.36	3.9 (2.77–5.48)	7.82	5.47E-15		1.24	3.45 (2.33–5.09)	6.22	5.08E-10*
Tobacco smoking history (ever vs. never)	0.22	1.25 (0.79–1.96)	0.946	3.44E-01
Risk score	0.331	1.39 (1.18–1.64)	3.94	8.17E-05*		0.356	1.43 (1.19–1.71)	3.84	1.22E-04*

^, multivariate Cox analysis: Concordance =0.759 (se =0.042), Likelihood ratio test =32.44 on 7 df, P=3e-05, Wald test =29.03 on 7 df, P=1e-04, Score (logrank) test =30.94 on 7 df, P=6e-05; #, multivariate Cox analysis: Concordance = 0.725 (se =0.031), Likelihood ratio test = 43.39 on 6 df, P=1e-07, Wald test =41.16 on 6 df, P=3e-07, Score (logrank) test =46.29 on 6 df, P=3e-08; &, multivariate Cox analysis: Concordance =0.715 (se =0.024), Likelihood ratio test =61.28 on 7 df, P=8e-11, Wald test =59.97 on 7 df, P=2e-10, Score (logrank) test =63.51 on 7 df, P=3e-11; $, multivariate Cox analysis: Concordance =0.716 (se =0.018), Likelihood ratio test =125 on 10 df, P=<2e-16, Wald test =129.4 on 10 df, P=<2e-16, Score (logrank) test =143.4 on 10 df, P=<2e-16. *, indicates P value <0.05.

Then, we established a ROC curve to evaluate the performance of the ten-gene signature in predicting LUAD overall survival rates. As shown in , the training cohort’s ROC analysis demonstrated that the AUC of the ten-gene signature risk score was 0.728, the highest among all clinical characteristics. The AUC of the ten-gene signature in the GSE72094 cohort also reached 0.774, which was the leading variable compared to other factors (). In the GSE68465 cohort, the AUC of the gene signature was 0.663, which was the best and superior to that of the N classification (AUC =0.657). Interestingly, when combined with the risk score and N classification for analysis, the AUC reached 0.711 ().

Figure 5

ROC analysis of the ten-gene signature risk score. ROC curves show the sensitivity and specificity of the prognosis prediction by the ten-gene signature risk score, age, sex, race, T classification, N classification, M classification, tumor stage, race, KRAS status, etc., or risk score + N classification in the training cohort (A) and two validation cohorts (B,C) based on overall survival. ROC, receiver operating characteristic; AUC, area under the ROC curve.

Identification of the correlation between autophagy and the ten-gene signature

We conducted a correlation Pearson correlation coefficient to evaluate the relationship between autophagy-related genes and the ten-gene signature risk score. Of the 490 autophagy-related genes, 308 (62.85%) were significantly correlated with risk scores, of which 129 were positively correlated, and 179 were negatively correlated (https://cdn.amegroups.cn/static/public/atm-20-7936-5.pdf). As shown in , GAPDH, ERO1 L, BIRC5, SPHK1, VDAC1, ITGB1, NAMPT, HK2, CISD2, and GNAI3 were the top ten autophagy-related genes that were positively correlated with the risk score, while DAPK2, NOD1, LRRK2, DLC1, TMEM150A, PRKCD, TFEB, SESN1, DRAM1, and KAT8 were the ten most enriched genes that had negative relationships with the ten-gene signature risk score.

Figure 6

The correlations between the ten-gene signature and autophagy-related genes in the training (A,B) cohort. Each graph’s blue line fits a linear model that indicates the proportional trend of the expression level of each gene and the risk score. The shading around the blue line represents the 95% confidence interval. Pearson coefficient was applied for correlation testing.

GSEA with the ten-gene signature

Given the negative correlation between the risk score and the outcomes, GSEA was performed between the high-risk and low-risk groups. As shown in https://cdn.amegroups.cn/static/public/atm-20-7936-6.pdf and , the enriched gene sets in the high-risk group were mainly involved in pathways related to glycolysis, mTORC1, MYC, G2/M checkpoint, unfolded protein response, E2F, hypoxia, mitotic spindle assembly, epithelial-mesenchymal transition, and late response to estrogen. Gene sets concerned with the metabolism of bile acids and salts, metabolism of heme, and genes downregulated by KRAS activation were primarily enriched in the low-risk group ( and https://cdn.amegroups.cn/static/public/atm-20-7936-7.pdf).

Figure 7

Gene set enrichment analysis performed using the HALLMARK collection (A,B). Gene sets with | NES | >1, NOM p-val <0.05, and FDR q-val <0.25 were considered significantly enriched.

Identification of the relationship between the ten-gene signature and 22 TICs

The CIBERSORT algorithm was applied to evaluate the proportion of tumor-infiltrating immune subpopulations to better study how the ten-gene signature and the immune microenvironment interact. The relative content distribution of the 22 TICs in the TCGA-LUAD cohort and the correlation between the 22 TICs are shown in .

Figure 8

TIC profile in tumor samples in the training cohort and correlation analysis. (A) Bar plot showing the relative content distribution of 22 TICs in high- and low-risk LUAD samples in the training cohort. Each column represents a sample. (B) Heatmap showing the correlation between 22 TICs. The shade and numeric in each small color box indicate the coefficient between two TICs. The coefficients covered by the X shape are not significant. Pearson coefficient was applied for correlation tests. P value <0.05 is the cutoff. TIC, tumor-infiltrating immune cell; LUAD, lung adenocarcinoma.

Combining the results of the difference analysis () and the correlation analysis (, https://cdn.amegroups.cn/static/public/atm-20-7936-8.pdf) identified a total of eight TICs associated with the signature (). Among the eight TICs, M0 macrophages, M1 macrophages, activated memory CD4 T cells, and activated mast cells were positively correlated with the risk score. Monocytes, resting memory CD4 T cells, resting mast cells, and resting dendritic cells were negatively correlated with the risk score.

Figure 9

Relationship between TICs and ten-gene signature risk score. (A) The violin plot shows the ratio differentiation of each of 22 TICs between the high- and low-risk groups. Wilcoxon rank-sum was applied for the significance test. (B) The correlations between the TICs and ten-gene signature risk score (only correlations with significance were plotted). Each graph’s blue line fits a linear model that indicates the proportional trend of the TICs and the risk score. The shading around the blue line represents the 95% confidence interval. Spearman coefficient was applied for the correlation test. (C) The Venn diagram shows that ten TICs have a pronounced correlation with the risk score, which is determined by the results of the violin and the scatter plot. P value <0.05 is the cutoff. TIC, tumor-infiltrating immune cell.

Furthermore, to examine the prognostic capacity of each TIC, univariate Cox and Kaplan-Meier analyses based on overall survival were conducted. As shown in , Kaplan-Meier analysis ( and https://cdn.amegroups.cn/static/public/atm-20-7936-9.pdf) indicated that resting mast cells and resting dendritic cells can predict the survival rate of LUAD. Additionally, the univariate Cox regression model () highlighted that resting mast cells, resting dendritic cells, and M0 macrophages impacted prognosis. From the above survival analysis, resting mast cells and resting dendritic cells can be identified as having a potential prognostic capacity in LUAD.

Figure 10

Evaluation of the prognostic ability of 22 TICs. (A) Kaplan-Meier survival curves. Only graphs with a P value <0.05 in the log-rank test were plotted. (B) Univariate Cox regression model built for 22 TICs based on overall survival. The asterisks shown in plot B indicate a P value <0.05.

The analysis in this section found that resting mast cells and resting dendritic cells not only had pronounced correlations with the gene signature but also predicted the outcomes of LUAD. Therefore, the significant infiltration of resting mast cells and resting dendritic cells potentially plays a vital role in the prognostic ability of the ten-gene signature in LUAD.

Discussion

In this study, we mined the public databases of the TCGA and GEO and constructed a ten-gene signature related to ferroptosis that can be used to predict the prognosis of LUAD. The LASSO Cox regression model was applied after discovering the potential ferroptosis-related prognostic genes using univariate Cox analyses in the TCGA-LUAD cohort. A ten-gene signature associated with LUAD prognosis was generated. The ten-gene signature was then applied to the training cohort and two validation cohorts by performing univariate and multivariate Cox analyses, ROC analysis, and Kaplan-Meier analysis. It was found that the ten-gene signature had a powerful prognostic ability, which also showed the effectiveness and broadness of the gene signature in LUAD prediction. In the subsequent correlation analysis, the ten-gene signature was found to be correlated with most autophagy-related genes. GSEA revealed essential pathways related to the gene signature. The TIC analysis revealed that resting mast cells and resting dendritic cells may play crucial roles in the gene signature that affects prognosis in LUAD. We are the first to utilize ferroptosis-related genes for training and validation in two independent cohorts (more than 400 cases each) to evaluate LUAD prognostic gene signatures. We anticipate that this study will guide future LUAD research.

Ferroptosis is a type of iron-dependent programmed cell death that is characterized by the accumulation of lipid peroxides and is genetically and biochemically different from other forms of regulatory cell death (such as apoptosis) (7). Ferroptosis-inducing factors can directly or indirectly affect glutathione peroxidase through different pathways, resulting in a decrease in antioxidant capacity and accumulation of lipid reactive oxygen species (ROS) in cells, ultimately leading to oxidative cell death (7). Recent studies have shown that ferroptosis is closely linked to many diseases, such as cancers, blood diseases, nervous system diseases, kidney diseases, and local ischemia-reperfusion injury (7). The increasing number of studies investigating ferroptosis in cancer have revealed its potential as a therapeutic target (27). The expression of FSP1 is correlated with the ferroptosis resistance of lung cancer cell lines, indicating that the upregulation of FSP1 expression is a strategy of ferroptosis escape in lung cancer (13). Additionally, in non-small cell lung cancer (NSCLC) cell lines, it was shown that the level of MAPK pathway activity correlates with sensitivity to ferroptosis induced by cystine deprivation (28). Interestingly, LSH affects metabolic genes through chromatin modification, thereby inhibiting hypertrophy and promoting lung tumors (29). RNA sequencing in NSCLC cells showed that the expression of SLC7A11, an essential gene associated with ferroptosis through its role in controlling iron concentrations, can be downregulated by XAV939 (an inhibitor of NSCLC), a target gene of lncRNAs, which suppressed the development of NSCLC via ferroptosis-mediated pathways (30). LSH promotes the expression of LINC00336 by upregulating ELAVL1 expression through the p53 signaling pathway in lung cancer (31,32). LINC00336 interacts with ELAVL1 to reduce intracellular iron and lipid ROS levels, thereby acting as a critical inhibitor of ferroptosis in carcinogenesis (31,32). These findings have highlighted the importance of ferroptosis in lung cancer therapeutics, but the roles of ferroptosis in tumorigenesis and development remain unclear. Recently, many studies have begun to mine the prognostic gene signature related to ferroptosis in tumor information available from public databases. For example, Liu et al. investigated ferroptotic gene expression patterns in patients with glioma, identified their relationship to patient outcome, and confirmed that the ferroptosis-related nineteen-gene signature they discovered could predict glioma cell death and glioma disease progression (11). Liang et al. defined a novel prognostic model of 10 ferroptosis-related genes, which proved to be independently associated with overall survival in both the derivation and validation cohorts, providing insight into the prediction of hepatocellular carcinoma prognosis (15). However, a prognostic gene signature had yet to be determined in LUAD. Our study conducted in-depth research to fill this gap and discovered a ferroptosis-related ten-gene signature that is closely associated with LUAD prognosis.

The ferroptosis-related ten-gene signature that we discovered showed strong prognostic prediction capabilities in the training cohort and two independent validation cohorts after being tested by a variety of statistical methods. Our signature was composed of the following ten genes: PHKG2, PGD, PEBP1, NCOA4, GLS2, CISD1, ATP5G3, ALOX15, ALOX12B, and ACSL3. In the signature model, ATP5G3, PGD, ACSL3, CISD1, and ALOX12B were unfavorable genes for LUAD prognosis, while other genes showed a protective effect. ATP5G3 was reported to be involved in diverse processes in human cancer cells and suppressed erastin-induced ferroptosis upon knockdown (33). PGD is involved in the pentose phosphate pathway, and its knockdown suppresses elastin-induced ferroptosis in the NSCLC cell line Calu-1 (33). ACSL3-dependent MUFA metabolism is a critical regulator of ferroptotic cell death, and ACSL3 activity specifically promotes a ferroptosis-resistant cell state (34). CISD1, located in the mitochondrial outer membrane, inhibits the uptake of iron ions by mitochondria and blocks ferroptosis (35,36). Lipoxygenases, including ALOX15 and ALOX12B, are a family of non-heme iron enzymes involved in generating leukotrienes from arachidonic acid (31). Several lipoxygenases are engaged in promoting ferroptosis in human cells (31). The ALOX15 protein is involved in FIN‐induced ferroptosis in cancer cells (37). In addition, the P53-SAT1-ALOX15 pathway is also involved in the regulation of ferroptosis (38). SAT1 is a transcription target of P53 and a crucial rate-limiting enzyme of polyamine catabolism. The activation of SAT1 induces lipid peroxidation and ferroptosis caused by ROS, which is closely related to the expression level of ALOX15 (7). PHKG2 positively regulates ferroptosis through the modulation of available iron, and PHKG2 silencing may function as iron chelation. The detailed mechanism of PHKG2 in iron metabolism regulation is unknown and needs further investigation (39-41). Zhao et al. found that PEBP1 is pivotal for dynamic interactions between the ferroptosis cell death process and prosurvival autophagy. The concurrent activation of autophagy protects cells from ferroptosis death and mitochondrial DNA release (42). Furthermore, PEBP1 is a vital regulator of relevant intracellular signaling pathways, including signaling pathways controlled by EGFR, such as MAPK. It has been reported that PEBP1 can inhibit metastasis and has prognostic value in various solid tumors, including lung cancer (43). Recent studies show that NCOA4-mediated ferritinophagy modulates susceptibility to ferroptosis (44,45). The GLS2 gene is a transcriptional target of the tumor suppressor p53, and upregulation of GLS2 expression contributes to p53-dependent ferroptosis (14). A previous comprehensive bioinformatics analysis showed that the expression of GLS2 is decreased in lung cancer and is associated with a better prognosis (46), which is consistent with our research.

Autophagy is the natural, regulated mechanism that removes unnecessary or dysfunctional cellular components. It allows the orderly degradation and recycling of cellular components (47). The original study shows that ferroptosis is morphologically, biochemically, and genetically distinct from autophagy and other cell death types (48). However, more recent studies have demonstrated that the activation of ferroptosis is indeed dependent on the induction of autophagy (48). In addition, accumulating studies have revealed crosstalk between autophagy and ferroptosis at the molecular level (23). Autophagy is critical to the maintenance of glucose homeostasis and tumor growth in lung cancer (49). Numerous studies have demonstrated that autophagy is involved in epidermal growth factor receptor tyrosine kinase inhibitor (EGFR-TKI)-acquired resistance in NSCLC, partially due to the inhibition of the PI3K/AKT/mTOR signaling pathway (50-54). In this study, we found that the risk score correlated with most of the autophagy-related genes (62.85%, 308/40), which can be used to further measure the connection between the gene signature and LUAD and provide potential explanations and inspirations for further research and development of autophagy-targeted tumor therapy.

GSEA of the HALLMARK collection found that the glycolysis, mTORC1, MYC, and G2/M checkpoint gene sets were the most enriched in the high-risk group. Altered energy metabolism is a biochemical fingerprint of cancer cells representing one of the “hallmarks of cancer”. This metabolic phenotype is characterized by preferential dependence on glycolysis for energy production in an oxygen-independent manner (55,56). Mammalian target of rapamycin (mTOR) is a key signal node that integrates environmental cues to regulate cell survival, proliferation, and metabolism and is often deregulated in human cancers (57). mTOR kinase works in two functionally different complexes whose activity and substrate specificity are regulated by complex cofactors; namely, mTORC1 and mTORC2. The dysregulation of this central signaling pathway is related to a variety of human diseases, including cancer (57). Recent research in NSCLC cells showed that overactivated AKT-mTOR signaling directly increases the expression of PD-L1 (58). The MYC oncogene contributes to the genesis of many human cancers. Recent insights into its expression and function have led to new cancer therapeutic opportunities. MYC activation by bromodomain proteins could be inhibited by drug-like molecules, resulting in tumor inhibition in vivo (59). In addition, MYC is usually amplified in lung cancer cells, and MYC transfection can enhance the in vitro proliferation ability of human small cell lung cancer cells (60). The G2/M checkpoint is activated by DNA damage and by incompletely replicated DNA (61). This checkpoint prevents chromosome segregation when the chromosome is not intact (61). The signaling pathway leading to G2 arrest after DNA damage is frequently altered or mutated in human cancer (61,62). Deficient G2/M checkpoints increase the risk of lung cancer (63). The GSEA results detail the methods by which the genes in the signature participate in the progression of LUAD, which can benefit future targeted therapy research.

Moreover, the CIBERSORT algorithm-based TIC analysis revealed that resting mast cells and resting dendritic cells have strong prognostic capacity and a pronounced correlation with the gene signature, suggesting that the infiltration of resting mast cells and resting dendritic cells may play a key role in the gene signature prediction ability. Mast cells are immune cells that accumulate in tumors and their microenvironment during disease progression. They play a vital role in the early stages of the innate immune response to pathogens. Studies have shown that mast cells can trigger different mechanisms, affecting the functional homeostasis of the lung (64,65). Imada and colleagues showed that patients with LUAD with a higher number of mast cells have a worse prognosis than patients with a lower number of mast cells (66). Dendritic cells represent a heterogeneous group of innate immune cells that infiltrate tumors and process and present tumor-derived antigens to naïve T cells (67). Dendritic cells play a critical role in priming antitumor T cell immunity and thereby represent a major therapeutic target for cancer immunotherapy (67). Dendritic cells are the key factors providing protective immunity against lung tumors, and clinical trials have proven that dendritic cell function is reduced in lung cancer patients (68). It is evident that the immunoregulatory network may play a key role in the immune response’s failure to terminate tumors (68). Dendritic cell-based immunotherapy has become an important strategy for cancer treatment, mainly via the direct activation of cytotoxic T cell responses (68). The recent development of dendritic cell vaccines has furthered our understanding of the vital roles of dendritic cells in the control of tumor progression (68). Our research shows that resting mast cells and resting dendritic cells can potentially target the gene signature in LUAD treatments, suggesting that these cells should be given more attention in future research.

This study also has some limitations. The ten-gene signature came from retrospective data, and more prospective data are needed to prove its clinical utility. Due to the limited patient clinical characteristics included in the TCGA cohort and the other two independent cohorts, we could not perform specific clinical subgroup analyses. Additionally, there are currently no laboratory data to explain and support the relationship between these ten genes and the progression of LUAD. Therefore, further research is expected to reveal the connection between the ten-gene signature and LUAD development.

Conclusions

The present study defined a novel, robust ferroptosis-related ten-gene signature for LUAD. The signature was strongly associated with the outcome of LUAD and can precisely detect the prognostic risk. Additionally, we examined the reliability and applicability of the signature by applying it to two independent validation cohorts and identified the vital roles of resting mast cells and resting dendritic cells in the prognostic capacity of the signature, which potentially further advances the discovery of treatments for LUAD.

The article’s supplementary files as