A Nomogram Integrating Ferroptosis- and Immune-Related Biomarkers for Prediction of Overall Survival in Lung Adenocarcinoma

Ferroptosis plays a dual role in cancer, which is known to be affected to antitumor immune responses. However, the association between ferroptosis and antitumor immune responses is uncertain in lung adenocarcinoma (LUAD). In this work, 38 ferroptosis-related genes (FRGs) and 429 immune-related genes (IRGs) were identified as being differentially expressed between tumor and normal samples. Two risk score formulas consisting of seven FRGs and four IRGs, respectively, were developed by Lasso-penalized Cox regression and verified in the GSE13213 dataset. The CIBERSORT algorithm was used to estimate the relative abundance of immune cells in tumors. The correlation between FRGs and immune cells was evaluated using the TIMER database. The results indicated that the development of ferroptosis was synergistic with that of anti-tumor immunity in LUAD. The concordance index and calibration curves showed that the performance of a nomogram that combines clinical staging and risk scores is superior to that of models using a single prognostic factor. In conclusion, ferroptosis might be synergistic with anti-tumor immunity in LUAD. The combined nomogram could reliably predict the probability of overall survival of LUAD patients. These findings may be useful for future investigation of prognostic value and therapeutic potential related to ferroptosis and tumor immunity in LUAD.


INTRODUCTION
Lung cancer is the most common type of cancer and the leading cause of cancer-related death worldwide. An estimated 1.8 million deaths from lung cancer occurred in 2020 (Sung et al., 2021). Lung adenocarcinoma (LUAD) is the most prevalent subtype of lung cancer, accounting for about 40% (Ma et al., 2013). Although there are many therapeutic options for LUAD, overall prognosis remains poor (Herbst et al., 2018). Therefore, it is necessary to further explore factors related to survival time that may contribute to the development of more effective treatment methods for lung adenocarcinoma.
Ferroptosis is an iron-dependent cell death modality marked by the oxidative modification of the phospholipid membrane (Stockwell et al., 2017). The effect of ferroptosis on tumor growth is not clear and requires further investigation. On the one hand, ferroptosis inhibits tumor growth in mice (Badgley et al., 2020). For example, although pancreatic cancer cells are prone to resisting chemotherapy, they are highly sensitive to artemisinin-induced ferroptosis (Efferth, 2017). On the other hand, ferroptosis can promote tumor growth (Dai et al., 2020). Therefore, examining the association between ferroptosis marker genes and survival may increase understanding of the role played by ferroptosis in LUAD.
Moreover, there seems to be an inextricable link between ferroptosis and immunity. Levels of CD8 + and CD4 + T cells cannot increase if they lack glutathione peroxidase 4 (GPX4), a key regulator of ferroptosis (Matsushita et al., 2015). Induction of ferroptosis is related to release of PGE2 (Yang et al., 2014), which attenuates antitumor immunity by affecting cDC1s and NK cells (Wang and DuBois, 2015). The synergistic effect of ferroptosis and immune regulation might not only inhibit the primary tumor but also stimulate immune responses in combination with immune checkpoint blockade (Li and Rong, 2020). In a word, FRGs and immunity have a profound impact on each other, and it is necessary to do more research to understand their relationship.
In this work, we constructed two formulas for calculating ferroptosis-related and immune-related risk scores using The Cancer Genome Atlas (TCGA) database. LUAD patients from the Gene Expression Omnibus (GEO) database were used for validation. Some immune cells had the same infiltration trend as risk score increased in two prognostic multigene signatures. In addition, a nomogram combining tumor stage, ferroptosisrelated risk score, and immune-related risk score was developed for more accurate prediction of patient survival. Finally, the correlation between ferroptosis-related genes (FRGs) and immune-related genes (IRGs) was assessed.

Data Acquisition
Gene expression data, including 497 tumor samples and 54 normal samples, in the HTSeq-FPKM format were obtained from the TCGA database using the GDC tool. 1 We also downloaded clinical data at the GDC portal. Meanwhile, the microarray dataset GSE13213 consisting of 117 patients with LUAD was downloaded from the GEO database 2 as an external validation dataset.
Probes were annotated based on annotation files. In total, 123 FRGs and 2,483 IRGs were obtained from the FerrDb database 3 and ImmPort Resources 4 (Supplementary Tables S1, S2), respectively. There was no ethical conflict since all data was obtained from public databases.

Construction of A Predictive Nomogram
The R/Bioconductor package "limma," which provides an integrated solution for analyzing gene expression data and is a popular choice for gene discovery through differential expression analyses, was used to identify differentially expressed FRGs and IRGs. The results were visualized in a volcano map using the popular drawing packages "ggplot2" and "ggrepel." |log2 foldchange| > 1 and FDR < 0.05 were considered significant. Prognosis-related genes were identified using univariate Cox proportional hazards regression analysis (p < 0.05). "Venn" R package is a common way to compare different datasets and identify and visualize intersections between sets. It was used to identify intersections between differentially expressed genes (DEGs) and prognosis-related genes as candidates for risk scoring.
To reduce the risk of overfitting, we applied LASSO-Cox regression to construct survival-predicting models. "Glmnet" is the most widely used package for LASSO analysis, and can generate a variety of models, including binary and multinomial logistic regression models, Poisson models, Cox proportional hazards models and SVM models. The risk score was calculated using "glmnet" and "survival" packages. The Cox coefficient and expression levels of the prognosis-related DEGs were extracted to calculate risk scores. The formula was as follows: Coei Expi , where N represents gene number, Coei represents coefficient value and Expi represents gene expression level. Based on the median risk score, LUAD patients were dichotomized into low-and high-risk groups. Considering the significant batch effect between the TCGA and GEO dataset, we used the median value of each to divide the high and low risk groups. The predictive ability of the model was evaluated using the log-rank test and receiver operating characteristic (ROC) analysis. The t-distributed stochastic neighbor embedding (t-SNE) test was implemented in "Rtsne" package to visualize clustering. Principal component analysis (PCA) was completed using the prcomp function. A nomogram was constructed to predict overall survival (OS) based on the results of multivariate Cox regression. The concordance index (C-index) was calculated to assess the stability of the nomogram by 1,000 bootstrap replicates. The performance of the prognostic nomogram was evaluated by plotting calibration curves.

Functional Annotation Analysis and Evaluation of Immune Cell Infiltration
In order to identify different pathways between the two groups, DEGs among groups were analyzed using R package "clusterProfiler" for Gene Ontology (GO). The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed by the same procedure. CIBERSORT is a gene-based deconvolution algorithm to estimate the abundance of any of 22 human immune cell types. We loaded "e1071" package to execute this algorithm to quantify the distribution of types of infiltrating immune cells in lung adenocarcinoma samples. The Wilcox test was conducted to identify disparities in the infiltration levels of immune cells between different risk groups, including B cells, plasma cells, T cells, natural killer cells, monocytes, macrophages, dendritic cells, mast cells, eosinophils, and neutrophils. Inter-group differences were identified by "limma" package and displayed in Violin plots using "Vioplot" package. The TIMER web server 5 is a comprehensive resource for analysis of immune infiltrates in multiple cancers. We analyzed the association between gene expression and abundance of infiltrating immune cells using gene modules.
Total RNA was extracted using the RNeasy mini kit (Qiagen, United States) and reverse-transcribed with the iScript cDNA Synthesis Kit. qRT-PCR was performed using the CFX96 Real-Time System (Bio-Rad, Hercules, CA, United States). All samples were tested in triplicate. Primers were purchased from Ribobio (Guangzhou, China) and are listed in Supplementary Table S3.

Statistical Analysis
The Wilcoxon rank-sum test was used to compare the gene expression in tumor and normal tissues. Risk scores were validated as independent prognostic factors by the application of Cox regression. Pearson correlation analysis was implemented to assess the correlation of gene expression levels. The statistical analysis tool used in this study was R software 4.0.2. All statistical tests were two-tailed. A value of p < 0.05 was considered statistically significant.
We calculated the risk score of FRGs according to the following formula: risk score = (0.048216 × SLC7A11) + (0.161996 × DDIT4) + (0.004447 × SLC7A5) + (−0.068488 × GDF15) + (−0. 061102 × IL33) + (0.027458 × SLC2A1) + (0.140388 × RRM2). A total of 464 patients were divided into high-risk (N = 232) and low-risk (N = 232) groups with the median as the cutoff value. The patients in the low-risk group had significantly higher survival rates than those in the high-risk group (left panel of Figure 3A). The DEGs between the two groups were mainly related to cell cycle transition, mitosis, protein synthesis, and chemotaxis (Supplementary Table S8). The signaling pathways obtained by KEGG analysis were not only correlated with cell division and maturation, but also correlated with immunity (Supplementary Table S9).
In the ROC analysis, the AUCs in the first, second and third year were 0.681, 0.658, and 0.684, respectively (left panel of Figure 3B). PCA and t-SNE analysis showed that patients were clustered into distinct groups (Figures 3C,D). In the univariate Cox regression analyses, the risk score was significantly associated with OS in the FRG prognostic risk model (HR = 3.617, 95% CI = 2.346-5.577, p < 0.001; upper panel of Figure 3E). Multivariate Cox regression revealed that the classifier was an independent prognostic factor (HR = 3.150, 95% CI = 2.037-4.872, p < 0.001; upper panel of Figure 3F). After applying the risk score formula to the GSE13213 dataset, significant differences in OS between the two groups still existed (right panel of Figures 3A,B). Moreover, the risk score remained an independent prognostic factor in the GSE13213 dataset (bottom panel of Figures 3E,F).
The prognostic risk score of IRGs = (−0.043721 × HLA− DRB5) + (−0.055446 × SFTPD) + (−0.053414 × PTGDS) + (−0. 003929 × S100B). As shown in the left panel of Figure 4A, patients with high-risk scores tended to have shorter survival times than those with low-risk scores. AUCs in the first, second, and third-year reached 0.679, 0.603, and 0.610, respectively (left panel of Figure 4B). Patients were also distributed in 2 regions (Figures 4C,D). In addition, in the univariate Cox regression analysis, the IRGs prognostic risk score was a variable closely related to prognosis (HR = 3.859, 95% CI = 2.002-7.440, p < 0.001; upper panel of Figure 4E). After multivariate Cox regression analysis of multiple factors, the IRGs prognostic risk model remained a reliable independent prognostic factor (HR = 3.818, 95% CI = 1.928-7.563, p < 0.001; upper panel of Figure 4F). The validity of the model was also verified in the GSE13213 dataset (right panel of Figures 4A,B; bottom panel of Figures 4E,F).

Correlation of Immune Cell Infiltration With Risk Scores and Gene Expression
To better understand the relationship between risk score and immune response, we calculated the proportions of 22 immune cells ( Figure 5A). In addition, the correlation matrix of immune cells showed that the infiltration level of activated memory CD8 + T cells was highly correlated with CD4 + T cells (Figure 5B). In addition, the level of resting mast cells highly correlated with that of monocytes was also high. Subsequently, immune cell infiltration was compared between high and low-risk groups. There were similar trends of differences in infiltration levels for memory B cells, activated memory CD4 + T cells, resting memory CD4 + T cells, monocytes, M0 macrophages, resting mast cells, and resting dendritic cells in the 2 models ( Figure 5C). Next, the relationship between gene expression and immune cell subtype infiltration was further analyzed in the TIMER database (Figures 5D,E). SLC7A11 and IL33 showed a strong correlation with immune cell infiltration. In contrast to SLC7A11, IL33 promoted the infiltration of nearly all types of immune cell. All IRGs were positively correlated with the level of immune cell infiltration (Supplementary Figure S1).

Construction of the Nomogram
The above results above showed that tumor stage, FRGs risk score, and IRGs risk score were independent prognostic factors in LUAD. The c-indices of the stage, ferroptosis, immune and combined models were 0.663, 0.644, 0.616, and 0.712, respectively ( Figure 6A). Therefore, the combined model was selected for prediction of 1-, 3-, and 5-year OS rates ( Figure 6B). Calibration plots showed that the combined model performed well in predicting 1-and 3-year survival but not 5-year survival ( Figure 6C). Taken together, compared with models established using a single prognostic factor, the combined model was superior for short-term survival prediction, which might be beneficial to diagnosis and treatment.

Relationship Between the FRGs and the IRGs in the Prognostic Risk Models
To identify interactions between FRGs and IRGs, we performed a correlation analysis of gene expression ( Figure 7A). Interleukin 33 (IL33) was synergistically co-expressed with HLA-DRB5, SFTPD, PTGDS, and S100B, and their high expression levels helped to prolong patient survival ( Figure 7B). The association between IL33 and prostaglandin D2 synthase (PTGDS) was the strongest among all genes, followed by IL33 and surfactant protein D (SFTPD). All the other FRGs were negatively correlated with IRGs. This result was consistent with their opposite effects on survival time. Similar results were observed in the TIMER databases ( Figure 7C).

Validation of Gene Expression by qRT-PCR
We detected the levels of 7 FRGs and 4 IRGs in the BEAS-2B, PC9, and H1299 cell lines by using qRT-PCR (Figure 8). SLC7A11, DDIT4, SLC7A5, GDF15 SLC2A1, and HLA-DRB5 were highly expressed in LUAD cell lines. The mRNA levels of IL33, RRM2, and SFTPD were downregulated in LUAD cell lines. PTGDS was highly expressed in the H1299 cell line but low in PC9 cell line. The above results were generally consistent with the TCGA results, which indicates that our bioinformatics analysis was credible.

DISCUSSION
Non-small cell lung cancer (NSCLC) is one of the deadliest form of cancers. LUAD and lung squamous cell carcinoma (LUSC) are the predominant histological phenotypes of NSCLC, and they exhibit significant differences in morphologic differentiation, underlying drivers, and response to various therapies (Wilkerson et al., 2010;Cancer Genome Atlas Research Network, 2014). In recent years, the incidence of LUAD has significantly increased compared with LUSC (Kinoshita et al., 2016). The possible benefits and potential harm of inducing ferroptosis during treatment are receiving increasing attention. Cell death in LUAD in response to treatment with siramesine and lapatinib has been reported to be mediated by ferroptosis (Villalpando-Rodriguez et al., 2019). Additionally, triggering ferroptosis increased sensitivity to radiotherapy in human patient-derived models of LUAD (Ye et al., 2020). This suggests that new therapies that combine immunotherapy with regulation of ferroptosis may lead to improved outcomes. Therefore, studies exploring the interaction between ferroptosis and immunity in LUAD may have far-reaching implications. The ferroptosis-related risk score calculation formulas used here contained seven FRGs. The qRT-PCR results showed that five genes were upregulated and two genes were downregulated. Except for RRM2, the other results were consistent with TCGA results. SLC7A11 is overexpressed in multiple types of cancer (Shi et al., 2019), and suppressing transcription and protein expression of SLC7A11 can effectively induce ferroptosis in cancer cells (Chang et al., 2018). DDIT4 has potential not only as a prognostic biomarker (Ho et al., 2020) but also as a therapeutic target . Paradoxically, upregulation of DDIT4 may promote cellular ferroptosis (Dixon et al., 2014), and therefore its specific role in cancer merits further exploration. Inhibition of SLC7A5 can regulate amino acid transport and affect ferroptosis (Dixon et al., 2014). Frontiers in Genetics | www.frontiersin.org 7 September 2021 | Volume 12 | Article 706814 Ansari et al. (2020) found that there is a significant correlation between overexpression of SLC7A5 and specific immune cell subtypes. Their study not only provided clinical evidence that SLC7A5 in breast cancer can aid the personalization of anti-PD1/ PDL1 inhibition therapies but also suggested that targeting SLC7A5 may enhance the efficacy of anti-PDL1 immunotherapy. GDF15 may promote ferroptosis (Dixon et al., 2014), and some studies indicate that GDF15 plays an anti-cancer role, but other data suggest that it may promote tumor progression and metastasis (Tsui et al., 2015;Wang et al., 2019b). In a study of acute kidney injury, Martin-Sanchez et al. (2017) found that IL-33 release was associated with ferroptosis, leading them to hypothesize that ferroptosis may regulate inflammation by activating IL-33 in vivo. Moreover, tumor-derived IL33 enhanced the antitumor effects of checkpoint inhibitors . SLC2A1 can be activated by lymphoid-specific helicase (LSH) to inhibit ferroptosis (Jiang et al., 2017), and high SLC2A1 expression usually correlates with poorer patient outcomes (Kawamura et al., 2001;Tohma et al., 2005). RRM2 was reported to be highly expressed in liver cancer tissues and to prevent ferroptosis (Zhang et al., 2019), and upregulation of RRM2 in NSCLC cells promoted proliferation and chemotherapeutic resistance (Huang et al., 2019).
The immune-related risk score calculation formulas used here contained four IRGs. The HLA-DRB5 alleles were associated with cervical neoplasia through a linkage disequilibrium with amino acid variations and HLA-DRB1 alleles (Bao et al., 2018). HLA-DRB5 has also been considered a possible prognostic factor for gastric cancer (Hang et al., 2018). High expression of SFTPD might function to prevent progression of lung cancer (Yamaguchi et al., 2011). Increased SFTPD levels have been shown to be associated with fewer distant metastases and progression-free survival in LUAD that harbors EGFR mutations (Umeda et al., 2017). The low expression of SFTPD in H1299 cells and TCGA samples indicates that SFTPD may be a protective factor in lung adenocarcinoma patients but is inhibited in tumor tissue. Prostaglandin H2 is converted to prostaglandin D2 (PGD2) under the catalytic action of PTGDS (Fukuhara et al., 2012). PGD2 has previously been shown to inhibit migration of cancer cells (Shyu et al., 2013). This effect can be achieved by influencing immune responses (Fagerberg et al., 2014). S100B is significantly downregulated in esophageal squamous cell carcinoma, and may cause cell growth stagnation and apoptosis through synergistic action with p53 (Ji et al., 2004). CacyBP/SIP is a target protein of S100B and an inhibitor of gastric cancer (Ning et al., 2007). Ferroptosis seems to be involved in the immune response to tumors. We observed that expression of the FRGs that inhibit ferroptosis negatively correlated with expression of IRGs. IL33 expression was consistent with the presence of ferroptosis and positively correlated with expression of IRGs. These results might suggest that ferroptosis had a positive relationship with anti-tumor immunity in LUAD. It has been shown that reduced uptake of cystine results in increased ferroptosis of tumor cells after tumor immunotherapy, which helps improve anti-tumor efficacy (Wang et al., 2019a). This means that the immune system can drive ferroptosis to mediate inhibition of tumor growth. By calculating the abundance of immune cells, it can be shown that the two systems overlap in differences in the levels of multiple types of immune cells. Xu et al. (2020) found that IFN-γ secreted in tumor tissues by infiltrating lymphocytes can help downregulate expression of SLC7A11, thereby promoting tumor cell ferroptosis and decreasing tumor volume. Efimova et al. (2020) observed that dendritic cells can be induced to mature and active by phagocytosis of early ferroptotic cancer cells. Further exploration of the clinical relevance of ferroptosis-and immune-related marker genes, and potential connection between them, would help to identify more efficient diagnostic and therapeutic approaches to LUAD.
Although many studies have attempted to elucidate an association between ferroptosis and immunity, there is no report of combining two variables for predicting OS of LUAD patients. At present, prognosis in LUAD is based on tumor stage. In the present study, the nomogram composed of tumor staging combined with FRGs and IRGs-related risk scores provided more precise prediction of patient outcomes. However, there are several limitations. Firstly, the data used in this work were all downloaded from public databases. The results should also be externally validated with other primary data. Secondly, although we have found some similarity in immune function, further understanding of functional connections between the seven FRGs and four IRGs identified here requires additional experimental study.
In conclusion, our study identified seven FRGs and four IRGs that are differentially expressed in LUAD and significantly associated with prognosis. A nomogram that combines these sets of genes is more beneficial to individualized prognosis in LUAD than the current standard that uses a single prognostic factor. Using Pearson correlation analysis, we inferred that immune response is positively correlated with ferroptosis in LUAD. These results provide valuable information for development of new therapies that combine immunotherapy with ferroptosis-related drugs.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: TCGA database: https://portal.gdc.cancer. gov/repository. GSE13213 was from GEO database: https:// www.ncbi.nlm.nih.gov/geo/.

AUTHOR CONTRIBUTIONS
MC and XL designed the study, executed the bioinformatics analysis, and prepared figures and tables. YZ and YT collected the data. MC, XL, YZ, PS, and JL wrote parts of the manuscript. XH, LW, and KS revised the manuscript.