Construction and validation of a robust ferroptosis-associated gene signature predictive of prognosis in lung adenocarcinoma.

ABSTRACT: To construct and validate a ferroptosis-associated signature predictive of prognosis in lung adenocarcinoma (LUAD), and systematically evaluate the underlying molecular connections in cancer biology.We retrieved mRNAs sequencing profiles of LUAD from the cancer genome atlas (TCGA) data portal and clinical information from the cBio Cancer Genomics Portal. The differentially expressed ferroptosis-associated genes (DEFAGs) were screened between normal samples and LUAD by packages "limma" in R. Then the total TCGA cohort was randomly divided into training set and testing set. Based on the training set, a DEFAG signature was built and further validated in the test set, the total TCGA cohort and other independent cohorts from the gene expression omnibus data portal. A nomogram was constructed and validated, and the correlation between high-risk group and cancer biology was further evaluated.We initially identified 68 DEFAGs from TCGA cohort. A 6 DEFAG signature was built and further validated in the test set, the total TCGA cohort and other 2 independent cohorts including GSE31210 and GSE72094 from gene expression omnibus data portal. Further exploration indicated that high-risk group combined with TP53 mutation harbored the most unfavorable prognosis while low-risk group with TP53 wild-type status had the most favorable survival advantage over other groups. Moreover, high-risk group was associated with higher cancer stemness, tumor mutation burden, and CD274 (programmed cell death 1 ligand 1) expression.We constructed a robust ferroptosis-associated gene signature and a nomogram predictive of prognosis in LUAD, and provided a new perspective on associations between ferroptosis and cancer.

Introduction

Lung cancer remains a major public health problem worldwide and is the leading cause of cancer-related mortality in human, accounting for 228,820 estimated new cases and 135,720 estimated deaths in 2020 in the United States nationally.[ Lung adenocarcinoma (LUAD) is one of the most common nonsmall cell lung carcinoma types, which account for almost 85% of all lung cancer cases.[ Although various therapeutic options are administrated including surgery, chemotherapy, radiation, molecular targeted therapy and immunotherapy, prognosis of lung cancer is still not satisfactory, with a 5-year survival rate of 57% for stage I disease and only 4% for those with stage IV disease.[ Unfortunately, approximately 70% of lung cancer patients are found to be locally advanced or metastatic at the time of diagnosis.[ With the rapid development of molecular biology, there still requires efforts in exploring predictive models based on molecular biomarkers for early diagnosis, predictive therapeutic responses to specific therapies and prognosis in lung cancer.

Ferroptosis, an iron-dependent form of regulated cell death, is characterized by accumulation of lipid peroxides and plays an important role in cancer biology.[ Ferroptosis is triggered by the inactivation of an essential metabolic process recognized as glutathione depletion and glutathione peroxidase 4 inactivation, leading to the accumulation of lipid peroxides and reactive oxygen species.[ Regarding inherited resistance to apoptosis of cancer, investigation targeting nonapoptotic regulated cell death process might provide clinical significance for cancer therapy.[

Although some publications investigated the correlations between ferroptosis and LUAD,[ this study aimed to explore multi-omics data on LUAD including mRNA expression, protein expression, mutational status as well as clinical information, construct and validate a ferroptosis-associated signature predictive of prognosis, and systematically evaluate the underlying molecular connections in cancer biology.

Materials and methods

Data collection

This study was not involved with new participants, so ethics committee approval was not necessary. The level 3 mRNA transcriptome profiling and mutation data of LUAD with corresponding clinical information were downloaded from the cancer genome atlas (TCGA) data portal (https://tcga-data.nci.nih.gov/tcga/) and cBioPortal (https://www.cbioportal.org/) in December, 2020. Other cohorts included GSE31210 with 226 LUADs cases and GSE72094 with 398 LUAD cases were retrieved from the gene expression omnibus portal (https://www.ncbi.nlm.nih.gov/geo/).[ A total of 91 ferroptosis-associated genes (Table S1, Supplemental Digital Content) were collected from previous literature and “WP_FERROPTOSIS” in the gene set enrichment analysis (GSEA) database (https://www.gsea-msigdb.org/gsea/index.jsp).[

Identification of DEFAGs and prognostic signature

A total of 68 differentially expressed ferroptosis-associated genes (DEFAGs) were identified between 54 normal tissues and 497 LUAD tissues from TCGA cohort at the threshold of |log2 fold change| > 0 and false discovery rate < 0.05 by packages “limma” in R in the whole TCGA LUAD cohort. A total of 457 LUAD cases with detailed follow-up information were randomly divided into training set (n = 229) and testing set (n = 228) by the “caret” package in R. Univariate and multivariate Cox regression analysis were used to explore overall survival (OS)-associated DEFAGs and construct a prognostic signature based on OS-associated DEFAGs by package “survival” in R in the training set. The risk scores were calculated based on the gene expressions and corresponding coefficients by the following formula:

The prognostic signature was further validated in the testing set and the whole TCGA LUAD set (n = 457), and other 2 independent cohorts (GSE31210 and GSE72094). Moreover, the mRNA and protein expression differences of the selected genes in the signature were compared between normal tissue and LUAD based on transcriptome profiling from TCGA portal and the human protein atlas database (https://www.proteinatlas.org/).[

Nomogram and clinical characteristics

The patients of the whole TCGA LUAD cohort were stratified into high- and low-risk groups based on the median value of risk scores. The OS differences were compared between high- and low-risk groups by stratification of various clinical characteristics including age, gender, T status, N status, M status, American joint committee on cancer (AJCC) stage. A nomogram consisting of clinical variables and the risk scores were constructed to predict the 1-, 3-, and 5-year OS of LUAD cases, while calibration plots and time-dependent receiver operating characteristic (ROC) curves were applied to evaluate the nomogram.

Enrichment analysis between high- and low- risk groups

The differently expressed genes (DEGs) (|log2 fold change| ≥ 1, false discovery rate < 0.05) between high- and low- risk groups based on the prognostic signature were identified by using the R package “limma”. The “clusterProfiler” R package was utilized to perform gene ontology analyses based on the DEGs between the high-risk and low-risk groups. The Kyoto encyclopedia of genes and genomes and potential activated pathways of DEGs in the high-risk group were analyzed using GSEA method.[

Prognostic signature and cancer biology

Tumor mutation burden (TMB) was defined as the total amount of mutations per million bases of tumor tissue from TCGA database via the GDC data portal (https://portal.gdc.cancer.gov/). The stemness index based on mRNA expression (mRNAsi) calculated using a one-class logistic regression machine learning algorithm was applied to LUAD samples from the TCGA cohort, and the mRNAsi scores were represented using β values ranging from zero (no gene expression) to one (complete gene expression) as a marker of cancer stemness.[

Statistical analysis

Wilcoxon rank-sum test was used to compare mRNA levels between tumor tissues and LUAD tissues, and programmed cell death 1 [PDCD1 (PD1)] expression, CD274 (programmed cell death 1 ligand 1 [PDL1]) expression, TMB, and mRNAsi between the high- and low- risk groups. Univariate and multivariate Cox regression analyses were used to select OS-associated factors and construct predictive models for prognosis. The OS between high- and low- risk groups was compared by Kaplan–Meier analysis with the log-rank test. All statistical analyses were performed with R software (Version 4.0.2, R Foundation for Statistical Computing, Vienna, Austria) or GraphPad Prism v7.00 (GraphPad Software Inc., USA).

Results

Construction and validation of a prognostic signature

The total flowchart was demonstrated in Figure 1. Among the 91 ferroptosis-associated genes, 68 DEFAGs were identified. As demonstrated in the heatmap (Fig. 2A) and volcano plot (Fig. 2B), 40 ferroptosis-associated genes were upregulated while 28 ferroptosis-associated genes were downregulated in LUAD tissues, when compared with normal tissues. Univariate Cox regression was conducted to analyze the correlations between the 68 DEFAGs and OS, and 9 OS-associated DEFAGs were obtained in the training set, as demonstrated in the forest plot (Fig. 2C). Based on the 9 OS-associated DEFAGs selected from the univariate Cox regression analysis, a multivariate Cox proportional hazard model consisting of 6 OS-associated DEFAGs, was established using stepwise regression. Risk scores of the 6 DEFAG signature were evaluated based on gene expressions and coefficients as follows: Risk score = (0.509∗CISD1) + (0.287∗ACSL3) + (−0.045∗PEBP1) + (−2.192∗ NOX1) + (0.549∗CHAC1) + (−12.221∗ GLS2). Then the median risk score was set as the cutoff value, the LUAD patients were stratified into the low-risk group and high-risk group in the training set, testing set and the whole LUAD set, respectively. The survival curves indicated significant OS advantage of low-risk group over high-risk group in the training set (Fig. 3A, P < .001), testing set (Fig. 3C, P = .011), and the whole LUAD set (Fig. 3E, P < .001). The time-dependent ROC curves showed that area under curve of 1 year, 2 years, and 3 years with moderate predictive performance (Fig. 3 B, D, and F). To evaluate whether the 6 DEFAG signature harbored similar predictive value in other cohorts, the same formula was used to generate risk scores of cases in GSE31210 and GSE 72094 cohorts. Consistent with the results in the TCGA cohort, patients with high-risk scores also had poor OS than those with low-risk scores (Figure S1A and C, Supplemental Digital Content), and predictive efficacy were also satisfactory (Figure S1B and D, Supplemental Digital Content).

Figure 1

The total flowchart of this study. KEGG = Kyoto encyclopedia of genes and genomes, GO = gene ontology, GSEA = gene set enrichment analysis, LUAD = lung adenocarcinoma, TCGA = the cancer genome atlas.

Figure 2

Exploration of prognostic DEFAGs in the TCGA cohort. (A) Heatmap indicated 68 DEFAGs between 54 normal tissues and 497 lung cancer tissues. (B) Volcano plot represented the upregulated and down-regulated DEFAGs in the TCGA cohort. (C) Forest plot of univariate Cox regression analysis of prognostic DEFAGs in training set. DEFAGs = differentially expressed ferroptosis-associated genes, FC = fold change, TCGA = the cancer genome atlas.

Figure 3

Predictive value of the 6 DEFAGs signature in the TCGA cohort. (A), (C), and (E) Kaplan–Meier curves for high- and low-risk groups in the training set, testing set, and the whole LUAD set, respectively. (B), (D), and (F) AUC of time-dependent ROC curves verified the prognostic performance of the risk score in the training set, testing set and the whole LUAD set, respectively. DEFAGs = differentially expressed ferroptosis-associated genes, LUAD = lung adenocarcinoma, ROC = receiver operator characteristic, TCGA = the cancer genome atlas.

Gene expressions and mutational status

The different mRNA and protein expressions of genes in the signature between normal and LUAD tissues were demonstrated in Figures 4 and 5. Compared with normal tissues, only PEBP1 were downregulated while other selected genes were upregulated in LUAD tissues. Moreover, the gene mutational status was shown in Figure 6, including the mutational proportion and hot spots, among which GLS2 was the most common mutational gene. Based on the mutational data from the dataset of LUAD (TCGA, PanCancer Atlas), the top 10 mutational genes were identified including TP53, TTN, MUC16, CSMD3, RYR2, LRP1B, ZFHX4, USH2A, KRAS, and XIRP2, as demonstrated in Figure S2, Supplemental Digital Content.

Figure 4

Different mRNA expressions of genes in the signature between normal and LUAD tissues. LUAD = lung adenocarcinoma.

Figure 6

Gene mutational status of genes in the signature, (A) mutational proportion, (B) hot spots.

Different protein expressions of genes in the signature between normal and LUAD tissues from the HPA database. HPA = human protein atlas, LUAD = lung adenocarcinoma.

Independent prognostic factors and construction of a nomogram

After stratified by clinical characteristics, the low-risk group still harbored OS advantage over high-risk group in patients with age <65 years, >65 years, female, male, T2, T3/4, N0, N1/2/3 status (Figure S3, Supplemental Digital Content), and M0, unknown M status, stage II, stage III–IV (Figure S4, Supplemental Digital Content). Univariate and multivariate Cox regression analysis were performed to explore the predictive factors for OS in the TCGA cohort, including clinical variables and risk scores based on the 6 DEFAG signature. The results indicated that the risk scores based on the signature was significantly associated with OS of LUAD as an independent factor (Figure S5A and B, Supplemental Digital Content). The ROC curve analysis showed that the risk score based on the signature was 0.695, only secondary to the stage and better than other clinical parameters (Figure S5C, Supplemental Digital Content), which indicated that our signature was an independent prognostic factor to predict OS of LUAD patients. Moreover, a nomogram consisted of risk score and other clinical parameters including age, gender, T stage, M stage, N stage, AJCC stage was constructed (Fig. 7A). The calibration plots demonstrated the predictive ability for 1-, 3-, and 5-year OS as shown in Figure 7 B–D. The area under curve values of the nomogram at 1-, 3-, and 5- year were 0.73, 0.723, and 0.68, respectively (Fig. 7E).

Figure 7

(A) Nomogram consisted of risk score and other clinical parameters including age, gender, T stage, M stage, N stage, AJCC stage. (B)–(D) Calibration plots of the predictive ability for 1-, 3-, and 5-year OS. (E) AUC values of the nomogram at 1-, 3-, and 5- year. AJCC = American joint committee on cancer, AUC = area under curve.

Functional enrichment analysis and GSEA

Gene ontology functional enrichment indicated that the DEGs between high- and low- risk groups in the TCGA cohort were mostly associated with humoral immune response, chromatin silencing, chromatin organization involved in negative regulation of transcription, negative regulation of gene expression, epigenetic (Fig. 8A). The Kyoto encyclopedia of genes and genomes functional enrichment indicated that the DEGs were enriched in cancer related pathways including IL-17 signaling pathway, necroptosis and Wnt signaling pathway. (Fig. 8B). The GSEA analysis demonstrated that the high-risk group were enriched in the process of cell cycle, homologous recombination, P53 signaling pathway, pancreatic cancer, pentose phosphate pathway, small cell lung cancer, ubiquitin mediated proteolysis (Fig. 9A), while the low-risk group were enriched in basal cell carcinoma, calcium signaling pathway, fatty acid metabolism (Fig. 9B).

Figure 8

Functional analyses between high- and low- risk groups in the TCGA cohort. (A) The GO enrichment of DEGs between the high-risk and low-risk groups. (B) The KEGG functional enrichment of DEGs between the high-risk and low-risk groups. DEGs = differently expressed genes, GO = gene ontology, KEGG = Kyoto encyclopedia of genes and genomes, TCGA = the cancer genome atlas.

Figure 9

The GSEA analysis demonstrated that the high-risk group were enriched in the process of cell cycle, homologous recombination, P53 signaling pathway, pancreatic cancer, pentose phosphate pathway, small cell lung cancer, ubiquitin mediated proteolysis (A), while the low-risk group were enriched in basal cell carcinoma, calcium signaling pathway, fatty acid metabolism (B). GSEA = gene set enrichment analysis, KEGG = Kyoto encyclopedia of genes and genomes.

TP53 mutation, TMB, cancer stemness between high- and low-risk groups

Considering the most frequent mutation of TP53 in the total TCGA cohort and the P53 signaling pathway enriched in the high-risk group in the GSEA analysis, we explored TP53 mutation frequency differences between high- and low-risk groups in the training set, testing set and the whole LUAD set, and compared the survival differences under TP53 mutation status stratification, and found that both in low-risk and high-risk groups,TP53 mutation was an unfavorable factor affecting OS of LUAD patients. Moreover, TP53 mutation in high-risk group showed the most unfavorable survival curve while TP53 wild-type in low-risk group indicated the obvious survival advantage over other groups (Fig. 10). Compared with low-risk group, high risk group was associated with higher CD274 (PDL1) (Fig. 11A, P < .001), mRNAsi (Fig. 11D, P < .001), not significantly associated with PDCD1 (PD1) (Fig. 11B, P = .319), CTLA4 (Fig. 11C, P = .269), TMB (Fig. 11E, P = .0825).

Figure 10

Survival differences under TP53 mutation status stratification in low-risk and high-risk groups. TCGA = the cancer genome atlas.

Figure 11

The relationship of the signature with cancer biology. (A) CD274 (PDL1), (B) PDCD1 (PD1), (C) CTLA4, (D) cancer stemness (mRNAsi), and (E) TMB. PDCD1 (PD1) = programmed cell death 1, mRNAsi = stemness index based on mRNA expression, TMB = tumor mutation burden.

Discussions

Ferroptosis, as a form of lipid peroxidation-induced cell death, can be regulated in many ways, from altering the activity of antioxidant enzymes to the level of transcription factors.[ Induction of ferroptosis by directly targeting GPX4 and its compensatory members provided a promising option for cancer treatment.[ This study systematically investigated the expression differences and potential roles of the 91 ferroptosis-associated genes in LUAD. A total of 68 DEFAGs were identified in LUAD compared to normal tissues based on the RNA-seq data from the TCGA cohort. A robust ferroptosis-associated gene signature was constructed in the training set, and validated in multiple cohorts with satisfactory predictive efficacy.

The selected 6 ferroptosis-related genes could be categorized into iron metabolism (CISD1), lipid metabolism (PEBP1, ACSL3), oxidant metabolism (NOX1, CHAC1), and energy metabolism (GLS2).[ CISD1, an outer mitochondrial membrane protein, played a significant role in regulation of iron and reactive oxygen species (ROS) homeostasis and inhibiting mitochondrial iron uptake, lipid peroxidation as well as subsequent ferroptosis.[ PEBP1 was shown to contribute to ferroptosis by complexing with lipoxygenases and allowing them to produce lipid peroxides.[ ACSL3 activated monounsaturated fatty acids, which promoted a ferroptosis-resistant cell state by suppressing ROS accumulation at the plasma membrane and decreasing levels of phospholipids containing oxidizable polyunsaturated fatty acids.[ Moreover, a recent study found that lymph protected metastasizing melanoma cells from ferroptosis due to higher levels of glutathione and oleic acid and less free iron in lymph when compared with blood plasma. In an ACSL3-dependent manner, oleic acid protected melanoma cells from ferroptosis and increased their capacity to form metastatic tumors.[ As a core member of the mitochondrial glutaminases, GLS2 was identified as a transcriptional target of the tumor suppressor protein P53 and was responsible for P53-mediated oxygen consumption, increased cellular antioxidant function in cancer cells.[ Interestingly, among the 6 genes, only GLS2 expression was risky for OS, while other 5 genes were protective against mortality. However, coefficients in the signature were not consistent with gene roles and vary significantly; the conflict might be associated with different expression levels of included genes and the method of Cox regression, which was used to establish the best-fit signature without considering gene roles.

TP53 was one of the 91 ferroptosis-associated genes, and its expression levels were not associated with prognosis in LUAD. Further analysis of mutational data indicated that TP53 was the most frequent mutated gene in LUAD, with 52.1% mutational frequency. Moreover, P53 signaling pathway was related to the high-risk group based on the ferroptosis-associated gene signature, which suggested the tight connection between TP53 mutation and risk scores. Our study indicated that high-risk group was characterized with a greater proportion of TP53 mutation, and survival analysis demonstrated that TP53 mutation or not exerted a profound role in both low-risk and high-risk groups. Regarding the tumor suppressing character of TP53, the high-risk group with TP53 mutation harbored the most unfavorable prognosis, while TP53 wild-type in low-risk group was predictive of most favorable survival. Multi-omics analysis in this study provided a promising predictive model for cancer prognosis. Compared with previous literature,[ this study constructed a new predictive model consisting of 6 ferroptosis-associated genes, innovatively investigated the correlation between TP53 mutations and this model, and explored the value of this model combined with TP53 mutational status in prognosis prediction.

Clinicopathological characteristics were significant predictors of cancer treatment. After stratified by various clinical factors, the risk scores based on the ferroptosis-associated gene signature were still tightly connected with LUAD prognosis. Univariate and multivariate Cox regression analysis indicated that the risk scores based on the signature was significantly associated with OS of LUAD as an independent factor. A nomogram based on risk scores and clinical parameters including age, gender, T stage, M stage, N stage, AJCC stage was built with excellent predictive values in 1-, 3-, and 5-year OS, which might have great potential to guide therapeutic decision-making in LUAD.

Cancer immunotherapy based on immune checkpoint inhibitors (ICIs) has achieved promising however limited success, the optimal method to select potentially beneficial responders to ICIs is significant in the clinic. High TMB was found to be beneficial in patients treated with ICIs.[ Compared with low-risk group, high-risk group based on the signature was associated with higher cancer stemness, TMB, PDL1 expression, which might suggest better immunotherapeutic response to ICIs in the high-risk group. Moreover, the induction of ferroptosis combined with ICIs showed synergistically enhanced antitumor activity, which was a novel strategy to overcome ICI-resistance in cancer treatment.[

Conclusions

In this study, we constructed a robust ferroptosis-associated gene signature and a nomogram predictive of prognosis in LUAD, which was closely related to prognosis, TMB, cancer stemness, and provided a new perspective on associations between ferroptosis and LUAD.

Acknowledgments

The mRNA expressions and corresponding clinical data were obtained from TCGA data portal (https://tcga-data.nci.nih.gov/tcga/), the cBio Cancer Genomics Portal (http://cbioportal.org), and gene expression omnibus (GSE72094, GSE31210). The protein expression data was obtained from the human protein atlas database.

Author contributions

Xin Zhu designed this study; Mi Zhou and Xin Zhu collected the data; Xin Zhu and Mi Zhou performed the statistical analysis; Mi Zhou and Xin Zhu wrote the paper. All authors read and approved the final version of the paper.

Conceptualization: Xin Zhu.

Data curation: Mi Zhou, Xin Zhu.

Formal analysis: Mi Zhou.

Methodology: Mi Zhou.

Software: Mi Zhou, Xin Zhu.

Supervision: Xin Zhu.

Validation: Mi Zhou, Xin Zhu.

Writing – original draft: Mi Zhou.

Writing – review & editing: Xin Zhu.