Development and Validation of a Robust Ferroptosis-Related Prognostic Signature in Lung Adenocarcinoma

Background Lung adenocarcinoma (LUAD) is the most common subtype of non-small cell lung cancer. Ferroptosis is a newly recognized process of cell death, which is different from other forms of cell death in morphology, biochemistry, and genetics, and has played a vital role in cancer biology. This study aimed to identify a ferroptosis-related gene signature associated with LUAD prognosis. Methods Dataset TCGA-LUAD which came from the TCGA portal was taken as the training cohort. GSE72094 and GSE68465 from the GEO database were treated as validation cohorts. Two hundred fifty-nine ferroptosis-related genes were retrieved from the FerrDb database. In the training cohort, Kaplan–Meier and univariate Cox analyses were conducted for preliminary screening of ferroptosis-related genes with potential prognostic capacity. These genes then entered into the LASSO Cox regression model, constructing a gene signature. The latter was then evaluated in the training and validation cohorts via Kaplan–Meier, Cox, and ROC analyses. In addition, the correlations between risk score and autophagy were examined by Pearson correlation coefficient. The analyses of GSEA and immune infiltrating were performed for better studying the function annotation of the gene signature and the character of each kind of immune cells played in the tumor microenvironment. Results A 15-gene signature was found from the training cohort and validated by Kaplan–Meier and Cox regression analyses, revealing its independent prognosis value in LUAD. Moreover, the ROC analysis was conducted, confirming a strong predictive ability that this signature owned for LUAD prognosis. One hundred fifty-one of 222 (68.01%) autophagy-related genes were discovered significantly correlated with risk scores. Analyses of GSEA and immune infiltration exhibited in detail the specific pathways that associate with the 15-gene signature and identified the crucial roles of resting mast cells and resting dendritic cells owned in the prognosis of the 15-gene signature. Conclusion In this present study, a novel ferroptosis-related 15-gene signature (RELA, ACSL3, YWHAE, EIF2S1, CISD1, DDIT4, RRM2, PANX1, TLR4, ARNTL, LPIN1, HERPUD1, NCOA4, PEBP1, and GLS2) was built. It could accurately predict the prognosis of LUAD and was related to resting mast cells and resting dendritic cells, which provide potential for the personalized outcome prediction and the development of new therapies in LUAD population.


INTRODUCTION
Lung cancer is the leading cause of death from cancer. In the United States, there will be approximately 228,820 newly diagnosed cases and 135,720 deaths in 2020 (Siegel et al., 2020). 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 is made up of two major types, lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (Herbst et al., 2018). LUAD is the predominant histology, and the rates are still increasing (Siegel et al., 2020). In recent years, several treatment advances have been made, especially the advancement of targeted therapy and the emergence of immunotherapy (Li et al., 2018;Low et al., 2019). However, these two methods can only benefit a limited number of subtypes, and the overall survival rate of LUAD is still very low . Therefore, further studies on therapeutic monitoring and prognostic assessment of LUAD are crucial for clinicians and patients.
Ferroptosis is a novel type of cell death that was discovered in recent years and is usually accompanied by a large amount of iron accumulation and lipid peroxidation during the cell death process . It is characterized by increased mitochondrial membrane density and cell volume contraction, which is different from other morphological, biochemical, and genetically regulated cell deaths Tao et al., 2020). Recent research has shown that ferroptosis is closely associated with the pathophysiological process of many diseases, such as tumors, neurological disorders, ischemia-reperfusion injury, kidney injury, and blood diseases . For the past few years, the induction of ferroptosis has emerged as a promising therapeutic alternative to trigger cancer cell death, especially for malignancies that are resistant to traditional treatments (Hassannia et al., 2019;Liang et al., 2019;Bebber et al., 2020). Apart from ferroptosis-inducing agents, numerous genes have also been identified as modulators or markers of ferroptosis (Stockwell et al., 2017;Bersuker et al., 2019;Doll et al., 2019;Hassannia et al., 2019;Liu et al., 2020). The fast-growing studies of ferroptosis in cancer have boosted a perspective for its usage in cancer therapeutics (Mou et al., 2019). The expression of FSP1 is correlated with the ferroptosis resistance of lung cancer cell lines, indicating that the upregulation of FSP1 is a strategy of ferroptosis escape in lung cancer (Doll et al., 2019).
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 (Poursaitidis et al., 2017). Interestingly, LSH inhibits ferroptosis and promotes lung tumorigenesis by affecting metabolic genes through chromatin modification (Jiang et al., 2017). RNA sequencing in NSCLC cells showed that SLC7A11, a key gene associated with ferroptosis through its role in controlling iron concentrations, can be downregulated by XAV939 (an inhibitor of NSCLC), as the target genes of lncRNAs, and suppress the development of NSCLC via ferroptosis-mediated pathways (Yu et al., 2019). LSH promotes the expression of LINC00336 by upregulating ELAVL1 through the p53 signaling pathway in lung cancer Wu et al., 2020). LINC00336 acts as a crucial inhibitor of ferroptosis in carcinogenesis by decreasing intracellular levels of iron and lipid ROS through interacting with ELAVL1 (ELAV-like RNA-binding protein 1), which has been recognized as a novel regulator of ferroptosis Wu et al., 2020). The above evidence has highlighted the importance of ferroptosis in lung cancer therapeutics, but the roles of ferroptosis in tumorigenesis and development remain unclear.
Autophagy is the natural, regulated mechanism of the cell that removes unnecessary or dysfunctional components. It allows the orderly degradation and recycling of cellular components (Mizushima and Komatsu, 2011). The original study shows that ferroptosis is morphologically, biochemically, and genetically distinct from autophagy and other types of cell death (Kang and Tang, 2017). However, recent studies demonstrate that activation of ferroptosis is indeed dependent on the induction of autophagy (Kang and Tang, 2017). In addition, accumulating studies have revealed cross talk between autophagy and ferroptosis at the molecular level .
Currently, several studies were mining prognostic gene signatures related to ferroptosis in tumors from public databases (Liang et al., 2020;Liu et al., 2020). Liu confirmed that the ferroptosis-related 19-gene signature could predict glioma patient survival . Liang et al. discovered a novel ferroptosis-related prognostic gene signature for hepatocellular carcinoma (Liang et al., 2020). However, there is still no study to declare whether there is a ferroptosis-related prognostic gene signature able to predict LUAD outcomes. In order to fill this blank and widen the potential of diagnosis and therapy of LUAD, the present study performed comprehensive analyses utilizing TCGA and GEO, along with ferroptosis-related genes identified in previous studies to determine and validate the minimum number of potentially robust prognostic genes of LUAD. In addition, we determined the correlations between the signature and autophagy and studied the characteristics of the gene signature in the tumor microenvironment through GSEA and immune infiltration analysis.
FIGURE 1 | Flowchart of the study. LASSO, the least absolute shrinkage and selection operator Cox regression model; ROC, receiver operating characteristic; UM, uveal melanoma; TICs, tumor-infiltrating immune cells.
Frontiers in Cell and Developmental Biology | www.frontiersin.org

Mining From TCGA and GEO Databases
We selected the TCGA-LUAD dataset as the training cohort and only included cases that meet the following criteria: (1) gene expression data is available; (2) survival data is available; (3) and follow-up time is greater than 0 days. Finally, 500 LUADs were included, and their gene expression profile, survival data, and clinical characteristics were downloaded from GDC Xena Hub 1 . Besides, datasets of GSE72094 (n = 442) and GSE68465 (n = 443) were downloaded from the GEO database and taken as validation cohorts to examine the gene signature we trained.

Ferroptosis-Related Genes
FerrDb 2 is the world's first manually curated database for regulators and markers of ferroptosis and ferroptosis-disease associations (Zhou and Bao, 2020). In this database, genes of drivers, suppressors, and markers of ferroptosis were listed based on the knowledge of previous studies (Supplementary Table 1). The unique 259 genes (Supplementary Table 2) were identified as the ferroptosis-related genes and entered to the next procedure.

Construction and Validation of the Prognostic Ferroptosis-Related Gene Signature in LUAD
Kaplan-Meier and univariate Cox analyses were applied in the training cohort to determine the association between the gene expression and patients' overall survival for potential prognostic genes, and P-value < 0.05 was considered to be statistically significant. The genes in the overlapped part of potential prognostic genes and ferroptosis-related genes were identified as potential prognostic ferroptosis-related genes, which then entered into an overall survival-based LASSO Cox regression model with penalty parameter tuning conducted by 10−fold cross-validation to detect the best penalty parameter lambda (Tibshirani, 1997;Sauerbrei et al., 2007;Friedman et al., 2010;Goeman, 2010). Based on the best lambda value, a list of prognostic genes with coefficients was harvested. As shown in the below formula, the risk score of each LUAD case could be obtained based on the expression level of each prognostic gene and its corresponding coefficient. In the formula, n, Expi, and βi indicate the number of hub genes, gene expression level, and regression coefficient value, respectively.
In the training cohort, patients were divided into low-and high-risk groups by using the median risk score as a cutoff point, and the survival difference of the two groups was measured by Kaplan-Meier analysis. Also, Cox and ROC analyses were conducted for further assessing of the gene signature prognostic ability. Moreover, in the validation cohorts, the same formula and statistical methods was adopted to validate the prognostic capacity of the gene signature.

Relationships Between Gene Signature and Autophagy in LUAD
To explore the relationship between autophagy and our gene signature, we identified 222 autophagy-associated genes from the Human Autophagy Database (HADb 3 ), which contains an exhaustive, up-to-date list of human autophagy-related genes (Abdul Rahim et al., 2017) (Supplementary Table 3). The Pearson correlation coefficient was applied to assess the correlation between autophagy and gene signature risk score. In addition, the GEPIA2 online tool 4 was applied to plot survival heatmaps of top 10 correlated genes . P-value < 0.05 was considered statistically significant.

Gene-Set Enrichment Analysis (GSEA)
Gene-set enrichment analysis was performed based on Hallmark (v7.1 5 ) gene-set collections using GSEA software (v4.1.0 6 ) in the training cohort to uncover the functions and pathways in the differentially expressed genes (DEGs) between high-and low-risk groups. Gene sets with | NES | > 1, NOM p-value < 0.05, and FDR q-value < 0.25 were considered significant.

Identification of the Relationship Between Gene Signature and 22
Tumor-Infiltrating Immune Cells (TIC) The relative proportion of 22 TICs in the training cohort was calculated using the CIBERSORT algorithm (Thorsson et al., 2018;Newman et al., 2019). After the quality filtering (Pvalue < 0.05), 391 LUADs were selected for the next analysis. The correlations between 22 kinds of TICs were examined by the Pearson coefficient.
For identifying the relationship between the 22-TIC proportion and risk score, an integrated analysis of consisting of the Spearman coefficient and Wilcoxon rank-sum was applied. Based on the 22-TIC proportion and survival data, Cox and Kaplan-Meier analysis were deployed to screen 22 TICs with prognostic meaning. P-value < 0.05 was a statistically significant threshold.

Statistical Analysis
The LASSO analysis was conducted applying the "glmnet" R package. Kaplan-Meier analysis was conducted using "survival" and "survminer" R packages. Univariate and multivariate Cox proportional hazard regression analyses were conducted using "survival" R package. The ROC analysis was applied using "survivalROC" R package. P-value < 0.05 was a statistically significant threshold.

Cohort Characteristics
The flowchart of the present research is shown in Figure 1. Five hundred LUAD cases that came from the TCGA-LUAD cohort were taken as the training cohort. The datasets GSE72094 and GSE68465, containing 442 and 443 LUAD cases, respectively, were selected as the validation cohorts. The detailed clinical characteristics of cohorts included in this study are summarized in Table 1.

Identification of Prognostic Ferroptosis-Related Gene Signature From the Training Cohort
Kaplan-Meier and univariate Cox analyses were conducted on 500 LUAD cases in the training cohort to examine the prognostic relationship between each gene and overall survival.
Only the gene showing P-values < 0.05 in both Kaplan-Meier and univariate Cox analyses was considered having potential prognostic capacity. Finally, 4,138 genes were identified as potential prognostic genes (Supplementary Table 4). The 4,138 potential prognostic genes and 259 ferroptosis-related genes were intersected to obtain a list containing 29 ferroptosisrelated potential prognostic genes ( Table 2). An overall survivalbased LASSO Cox regression model was built using the 29 ferroptosis-related potential prognostic genes (Figure 2A). When 15 genes were included, the model achieved the best performance ( Figure 2B). The regression coefficient of each gene was calculated and shown in Table 3.

Prognostic Value of the 15-Gene Signature in the Training and Validation Cohorts
The risk score of each LUAD case was a linear combination of each 15-gene signature expression level and its risk coefficient.
Patients were sorted to high-and low-risk groups based on their median. The distribution of risk scores, outcome status, and gene profiles of the 15-gene signature in training and validation cohorts are shown in Figure 3. As demonstrated in the graph (Figures 3A,D), more deaths or events happened in the high-risk groups than those in their corresponding low-risk groups. We checked the performance of this 15gene signature in 5-year survival (Figures 3E,H) and found consistent patterns. The boxplots in Supplementary Figure 1 show the 15-gene expression distributions in the high-and low-risk groups. The analysis in the training cohort witnesses that ACSL3, YWHAE, DDIT4, PANX1, RELA, CISD1, EIF2S1, and RRM2 were overexpressed, while GLS2, PEBP1, ARNTL, NCOA4, LPIN1, HERPUD1, and TLR4 were downregulated in high-risk groups.
Univariate and multivariate Cox analyses were applied in the training and two validation cohorts based on overall survival or progression-free survival, using the available co-variables including risk score, gender, 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 15-gene signature among other clinicpathologic characteristics ( Table 4). In the overall survival-based Cox regression model of the training cohort, both univariate and multivariate results suggested that the 15-gene signature was a powerful player (HR = 3.77, 95% CI = 2.77-5.14, P-value = 4.59E-17, and HR = 6.52, 95% CI = 2.74-15.5, P-value = 2.20E-05, respectively). Consistent with that in the training cohort, in either  univariate or multivariate analysis, the 15-gene signature showed excellent ability in the other two independent validation cohorts in predicting overall survival (P-value < = 1.04E-04). Also, we utilized progression-free survival data in the training cohort to perform the Cox analysis and found that the 15-gene signature had the ability not only in univariate but also in multivariate models to predict the outcomes (P-value < = 2.84E-03). These evidences demonstrated that the 15-gene signature was to be an independent and strong variable of prognosis. We built ROC curves and time-dependent dynamic AUC comparisons with other variates that showed an independent prognostic value in the multivariate Cox analysis in Table 4 to evaluate how the 15-gene signature could behave in predicting outcome and whether risk score is superior to other variates. As shown in Figure 5A, the AUCs of the 15-gene risk score model performed in the training cohort were 0.724, 0.755, and 0.681 at 1, 3, and 5 years, respectively, and the time-dependent dynamic AUC plot exhibited that the risk score displayed the best performance among the independent variates at all time points within the 5-year period. Consistently, the AUCs of the risk score in the GSE72094 ( Figure 5C) and GSE68465 (Figure 5D) cohorts were > = 0.681 at 1, 3, and 5 years and were greater than other vital variates at any time within 5 years. Besides, we used progression-free survival data in the training cohort to evaluate the predicting ability of risk score and found a similar pattern that in the other three cohorts, risk score had the best capacity ( Figure 5B).

Identification of the Autophagy Correlation With the Gene Signature Risk Score
Moreover, we conducted Pearson correlation to evaluate the relationship between autophagy-related genes and the 15-gene signature risk score. Of the 222 autophagy-related genes, 151 (68.01%) were significantly correlated with risk scores, of which 74 were positively correlated and 77 were negatively correlated (Supplementary Table 5). As shown in Figure 6, GAPDH, BIRC5, ERO1L, EIF2S1, SPHK1, ATIC, GNAI3, NAMPT, EIF4EBP1, and FADD are the top 10 autophagy-related genes that positively corrected with the risk score; 8/10 showed a significant elevated hazard ratio in LUAD (Supplementary Figure 2A). ERN1, ATG16L2, CCR2, IKBKB, HSPB8, PRKCD, DAPK1, DRAM1, DLC1, and DAPK2 are the leading 10 that have negative relationships with the 15-gene signature risk score; three of them exhibited a decreased hazard ratio (Supplementary Figure 2B).

GSEA With the 15-Gene Signature
In view of the negative correlation between the 15-gene signature risk score level and the LUAD outcomes, GSEA was performed between the high-risk and low-risk groups. As displayed in Figure 7A and Supplementary Table 6, enriched gene sets of HALLMARK collection in the high-risk group were mainly involved in pathways related to glycolysis, unfolded protein response, mTORC1, MYC, G2/M checkpoint, E2F, DNA repair, mitotic spindle assembly, ultraviolet radiation, hypoxia, cholesterol homeostasis, and reactive oxygen species, whereas the gene set concerned with metabolism of bile acids and salts was primary enriched in the low-risk group (Figure 7B and Supplementary Table 6).

Identification of the Relationship Between the Fifteen-Gene Signature and 22 TICs
To better study how the 15-gene signature and the immune microenvironment interact, we used the CIBERSORT algorithm to evaluate the proportion of tumor-infiltrating immune subpopulations and made comprehensive comparisons with the risk score. The relative content distribution of 22 TICs in the TCGA-LUAD cohort and the correlation between 22 TICs are shown in Supplementary Figure 3.
Incorporating the results of difference analysis ( Figure 8A) and correlation analysis (Figure 8B and Supplementary Table 7), 10 TICs, including resting mast cells, activated mast cells, M0 macrophages, activated CD4 memory T cells, resting dendritic cells, resting memory CD4 T cells, neutrophils, B cell memory, monocytes, and regulatory T cells (Tregs), were identified associating with the 15-gene signature risk score ( Figure 8C). Among them, mast cells activated, macrophages M0, and T cells CD4 memory activated were positively correlated with risk score, the remaining negatively correlated.
Furthermore, the prognostic capacity of each TIC was examined using univariate Cox and Kaplan-Meier analyses. As shown in Figure 9, the univariate Cox regression model ( Figure 9A) highlighted that resting mast cells, resting dendritic cells, and M0 macrophages impacted the prognosis; besides, Kaplan-Meier analysis (Figure 9B and Supplementary Table 8) indicated that resting mast cells and resting dendritic cells can predict the survival of LUAD. It can be seen that resting mast cells and resting dendritic cells showed significance in both analyses and may have potential prognostic ability in LUAD. Also, resting mast cells and resting dendritic cells not only had prognostic value but also owned significant correlations with the risk score, as mentioned earlier. Therefore, the significant infiltration with resting mast cells and resting dendritic cells may play a vital role in contributing to the prognosis value of the 15-gene signature in LUAD.

DISCUSSION
In this study, we built a ferroptosis-related 15-gene signature for the prognosis of LUAD by comprehensively mining the TCGA and GEO databases. After discovering the potential ferroptosisrelated prognosis genes using Kaplan-Meier and univariate Cox analyses in the TCGA-LUAD cohort, the LASSO Cox regression model was applied, and a 15-gene signature was generated which was related to outcome of LUAD. By applying the 15-gene signature in the training and validation cohorts, pronounced statistical differences were seen in Kaplan-Meier analysis, univariate and multivariate Cox regression models, and ROC curves, demonstrating the effectiveness and broadness of the gene signature in predicting LUAD prognosis. In the following correlation analysis, the 15-gene signature was found FIGURE 4 | Kaplan-Meier curves of the 15-gene signature risk score in the training (A,B,E,F) and validation (C,D,G,H) cohorts. The middle part of each graph indicates the number of patients at risk. The differences between the high-and low-risk groups were measured by the two-side log-rank test with a P-value < 0.05. OS, overall survival; PFS, progression-free survival.
correlated with most autophagy-related genes. The GSEA and analysis of immune infiltration exhibited important pathways that relate to the 15-gene signature and vital roles that resting mast cells and resting dendritic cells may have played backing the 15-gene signature influencing the outcome of LUAD. Compared with previous studies discovering the prognostic gene signature in LUAD, we are the first to utilize ferroptosis-related genes for training and validation in two independent cohorts (more than 400 cases each). This work we have done aimed to present more hints in future LUAD research.
in tumor-stroma interactions and correlates strongly with the severity of tumor infiltration by inflammatory cells in NSCLC patients (Giopanou et al., 2015). ACSL3 activity was previously demonstrated specifically to promote a ferroptosis-resistant cell state . ACSL3 overexpression increased cell proliferation, migration, and invasion altering the metabolic properties of lung cancer cells and was associated with worse clinical outcomes in patients with high-grade NSCLC (Fernandez et al., 2020). YWHAE was upregulated in breast cancer , high-grade endometrial stromal sarcoma (Hemming et al., 2017), gastric cancer (Leal et al., 2016), and colorectal cancer (Bjeije et al., 2019) and associated with poor outcomes. However, the role of YWHAE in LUAD has not been adequately evaluated. The interaction between TIPRL and EIF2S1 leads to phosphorylation of EIF2S1 and activation of the EIF2S1-ATF4 pathway, thereby inducing autophagy and playing a role in affecting the prognosis of LUAD (Jeon et al., 2019). CISD1 inhibits ferroptosis by protection against mitochondrial lipid peroxidation (Yuan et al., 2016). A recent study revealed that the NEET protein CISD1 plays a critical role in promoting the proliferation of cancer cells, supporting tumor growth and metastasis (Mittler et al., 2019). The function of CISD1 in cancer cells was found to be dependent on the degree of lability of their 2Fe-2S clusters (Mittler et al., 2019), whereas the mechanism of CISD1 in LUAD remained unclear. In multiple malignancies, studies have shown that DDIT4 participates in tumorigenesis and impacts patient survival (Cheng et al., 2020). Mu et al. (2019) showed that inhibition of SIRT1/2 induces pro-survival autophagy via acetylation of HSPA5 and subsequent activation of ATF4 and DDIT4 to inhibit the mTOR signaling pathway in NSCLC cells. It is reported that RRM2 is involved in the progression of various cancers, including glioma, colorectal cancer, bladder cancer, and NSCLC . In a recent study, Su provided important evidence that PANX1 plays an important role in ferroptotic cell death (Su et al., 2019). In addition, Jalaleddine identified a role for PANX1 in breast cancer progression and that PANX1 mediates this tumorpromoting role by modification of the EMT pathway (Jalaleddine et al., 2019). However, the impact of PANX1 on LUAD and its prognosis remained unknown. The expression of TLR4 is highly observed in the cells of the immune system such as monocytes, lymphocytes, and splenocytes, but it is also expressed in lower levels in epithelial and endothelial cells as well as cancer cells (Shetab Boushehri and Lamprecht, 2018). TLR4 agonists have thus been, and are still being, wildly explored as highly potentially immunotherapeutic for the treatment of cancer (Shetab Boushehri and Lamprecht, 2018). Previous study has shown that ARNTL played a role in tumor suppression . Nevertheless, few studies have found that ARNTL has tumor promotion properties . LPIN1 plays an important role in cell proliferation and tumor development through the regulation of intracellular signaling pathways. However, the underlying molecular mechanisms and the specific signaling pathway of LPIN1 during tumor development remains to be elucidated (Kim et al., 2016). HERPUD1 protects against oxidative stress-induced apoptosis through downregulation of the inositol 1,4,5-trisphosphate receptor (Paredes et al., 2016). However, little is known about the links between HERPUD1 and cancers. Tang's research shows that overexpression of NCOA4 increases ferritin degradation and promotes ferroptosis (Hou et al., 2016). However, how NCOA4 affects the occurrence and progression of tumors is still largely unknown (Hou et al., 2016). It has been well-established that PEBP1 suppresses the metastatic spread of tumor cells; moreover, the downregulated expression of PEBP1 is observed in a number of human cancers (Xu et al., 2010). According to one recent study, downregulation of PEBP1 is associated with poor prognosis of hepatocellular carcinoma, while the relationship between PEBP1 and LUAD is still unknown (Xu et al., 2010). The GLS2 gene is a transcriptional target of p53, and in glioblastoma and liver cancer GLS2 has been described as a tumor suppressor (Lukey et al., 2019). The role of GLS2 in LUAD needs to be studied. Autophagy is the natural, regulated mechanism of the cell that removes unnecessary or dysfunctional components. It allows the orderly degradation and recycling of cellular components (Mizushima and Komatsu, 2011). The original study shows that ferroptosis is morphologically, biochemically, and genetically distinct from autophagy and other types of cell death (Kang and Tang, 2017). However, recent studies demonstrate that activation of ferroptosis is indeed dependent on the induction of autophagy (Kang and Tang, 2017). In addition, accumulating studies have revealed cross talk between autophagy and ferroptosis at the molecular level .
Autophagy has been shown to exert pleiotropy, exhibiting anticarcinogenic, pro-survival, and pro-apoptotic effects in different stages of lung cancers (Ryter and Choi, 2015). Recent efforts to develop lung cancer treatment strategies have focused on understanding the role of autophagy as a tumor-suppressing or tumor-promoting mechanism (Jeon et al., 2019). In this study, we discovered that the risk score correlated with more than half of the autophagy-related genes (68.01%, 151/222) and found that the top 10 correlated genes hold the same effect predicting the LUAD outcomes, which elaborated further the relation between the ferroptosis-related 15-gene signature and LUAD and also provided more possibilities for autophagytargeted strategies. The GSEA in HALLMARK collection found that gene sets about glycolysis, unfolded protein response, and mTORC1 were top enriched. Changes in energy metabolism are the biochemical fingerprints of cancer cells and represent one of the "hallmarks of cancer." This metabolic phenotype is characterized by preferentially relying on glycolysis to produce energy in FIGURE 8 | Relationship between TICs and 15-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 15-gene signature risk score (only correlations with significance were plotted). The blue line in each graph 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. The Spearman coefficient was applied for the correlation test. (C) The Venn diagram shows that 10 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. *, P-value < 0.05.
an oxygen-independent manner (Ganapathy-Kanniappan and Geschwind, 2013; Liberti and Locasale, 2016). Tumor cells (including lung cancer cells) ingesting a large amount of glucose will hinder the supply of nutrients to adjacent normal cells (Sotgia et al., 2011;Chang et al., 2020). Glycolysis can also induce deoxyribonucleic acid (DNA) mutations and peroxide production, both of which are conducive to the proliferation and metastasis of tumor cells (Sotgia et al., 2011;Chang et al., 2020). The unfolded protein response is a cellular stress response related to the endoplasmic reticulum stress (Madden et al., 2019). The internal pressure in the tumor (such as oncogenic activation) and the external pressure exerted by the tumor environment can increase the level of misfolded proteins in the endoplasmic reticulum, which triggers the activation of the unfolded protein response pathway (Madden et al., 2019). The tumor environment is an environment of hypoxia, acidity, and nutritional deficiency. The three arms of unfolded protein response are highly active in many types of cancers (including breast, lung, liver, and colorectal and glioma) (Madden et al., 2019). mTOR is a pathway found to be dysregulated in many disease states including cancer (lung cancer) (Ekman et al., 2012;Kim et al., 2017). Studies have shown that mTOR signal inhibition blocks tumor cell progression, disrupts angiogenesis, and induces apoptosis and autophagy (Ekman et al., 2012;Kim et al., 2017). Growth factors, stress, amino acids, energy, and oxygen give inputs to mTOR which activate the mTORC1 (Vicary and Roman, 2016;Murugan, 2019). Activated mTORC1 regulates stress/DNA damage and enhances nucleotide synthesis, protein synthesis, and metabolism that promote cell proliferation, cell survival metastasis, etc. (Vicary and Roman, 2016;Murugan, 2019). Targeting of mTOR is thus an attractive strategy in the development of therapeutic agents against lung cancer (Ekman et al., 2012). These GSEA results described in detail the ways and methods that the 15-gene signature participates in the progress of LUAD, which can benefit future targeted therapy research.
Moreover, the CIBERSORT algorithm-based TIC analysis discovered that resting mast cells and resting dendritic cells have strong prognostic capacity in LUAD and a significant correlation with the ferroptosis-related 15-gene signature risk score, revealing that the infiltration of resting mast cells and resting dendritic cells may play key roles affecting the prognostic ability of the 15-gene signature. Mast cells can be used as an important innate immune sentinel and have the ability to enhance the immune response mediated by T cells but in other cases also show the ability to suppress the immune response (Dudeck et al., 2019;Kaesler et al., 2019). Consistent with their functional plasticity, the number of mast cells in TME is reported to be associated with cancer progression and improvement in patient survival (Kaesler et al., 2019). In NSCLC, mast cell infiltration of tumor islets confers a survival advantage independently of tumor stage (Varricchi et al., 2017). In another study, it was found that only in stage I NSCLC were increased peritumoral mast cells associated with a better prognosis (Varricchi et al., 2017). Dendritic cells represent a heterogeneous group of innate immune cells that infiltrate and process tumors and present tumor-derived antigens to naive T cells (Wylie et al., 2019). Dendritic cells play a key role in triggering antitumor T cell immunity and therefore represent the main therapeutic target of cancer immunotherapy (Wylie et al., 2019). Dendritic cells are a key factor in providing protective immunity against lung tumors (Wang J.B. et al., 2019). Clinical trials have proved that the DC function of lung cancer patients is reduced (Wang J.B. et al., 2019). Based on our research, resting mast cells and resting dendritic cells have the potential to target the 15-gene signature for the treatment means of LUAD, and good efforts should be carried out to investigate these immune cells further.

CONCLUSION
Our study found a novel robust ferroptosis-related 15-gene signature for LUAD. The signature is strongly associated with the prognosis of LUAD and can precisely detect the LUAD risk level. Remarkably, we validated the reliability and applicability of this signature by applying two independent validation cohorts and identified the vital role resting mast cells and resting dendritic cells may interplay in the prognostic capacity of the gene signature, which could potentially advance the discovery of new treatments for LUAD.

AUTHOR CONTRIBUTIONS
AZ, JY, and CM organized and wrote the manuscript. FL and HL produced the figures and visualized the data. CM revised the manuscript. All authors reviewed the manuscript and approved the manuscript for publication.

FUNDING
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.  validation (B,C) cohorts. The box plot shows the expression differentiation of each of the fifteen gene between the highand low-risk groups. Wilcoxon rank-sum was applied for the significance test. P-value < 0.05 was considered statistically significant; ns: P-value > 0.05; * P-value < 0.05; * * P-value < 0.01; * * * P-value < 0.001.

ACKNOWLEDGMENTS
Supplementary Figure 2 | Survival heatmaps of the top 10 autophagy-related genes positively (A) and negatively (B) correlated with risk score in LUAD (GEPIA2 web patrol, http://gepia2.cancer-pku.cn/). The survival heatmaps show the hazard ratios in the logarithmic scale (log10) for different genes. The red and blue blocks denote higher and lower risks, respectively. The rectangles with frames mean the significant unfavorable and favorable results in prognostic analyses. LUAD, lung adenocarcinoma.   Table 4 | 4138 genes were significantly predicting the prognosis of LUAD patients by both Kaplan-Meier and univariate Cox regression analyses (P-value < 0.05).
Supplementary Table 5 | The correlations of 222 autophagy-related genes with the fifteen-gene signature risk score examined by the Pearson correlation coefficient.