Comprehensive Analysis of Ferroptosis-Related Markers for the Clinical and Biological Value in Gastric Cancer

Gastric cancer is a highly malignant tumor with poor survival rate. Ferroptosis, a newly defined regulated cell death, is closely related to several tumors. Introduction of ferroptosis is promising for cancer treatments. However, the predictive role of ferroptosis in GC remains elusive. In this study, we screened the ferroptosis-related genes which were differentially expressed between normal and GC tissues. Then, based on these differentially expressed genes (DEGs), the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regressions were applied to construct the 10-gene prognostic signature (SP1, MYB, ALDH3A2, KEAP1, AIFM2, ITGB4, TGFBR1, MAP1LC3B, NOX4, and ZFP36) in TCGA training dataset. Based on the median risk score, all GC patients in TCGA training dataset and GSE84437 testing dataset were classified into a high- or low-risk group. GC patients in the low-risk group showed significantly higher survival possibilities than those in the high-risk group (P < 0.001). Combined with the clinical characteristics, the risk score was proven as an independent factor for predicting the OS of GC patients. Besides, the GC patients in the high- or low-risk group showed significantly different GO and KEGG functional enrichments, somatic mutation, fractions of immune cells, and immunotherapy response. Then, the expression levels of these genes in signature were further verified in the GC cell lines and our own GC samples (30-paired tumor/normal tissues). Furthermore, the effects of ferroptosis inducer Erastin on these 10 ferroptosis-related genes in GC cell lines were also explored in our study. In conclusion, our study constructed a prognostic signature of 10 ferroptosis-related genes, which could well predict the prognosis and immunotherapy for GC patients.


Introduction
Gastric cancer (GC) is one of the most common malignant tumors in the world, ranking fifth for incidence (1,089,103 cases) and fourth for mortality (768,793 cases) globally in 2020 [1]. Besides, the incidence rate of GC is the highest in digestive malignant tumors in China [2]. Due to the comprehensive treatments in the last few decades, including curative surgery, chemoradiotherapy, targeted therapy, and immunotherapy, the prognosis of GC patients has improved a lot. However, most patients are diagnosed at advanced stages; the overall survival (OS) rate of 5 years remains less than 40% [3]. Currently, the prognosis of GC patients was based on the TNM staging system; nevertheless, the patients at the same stage could show obviously different prognosis. Therefore, it is necessary to identify novel and reliable biomarkers to accurately predict prognosis, to find potential therapeutic targets, and finally to improve the outcomes of GC patients.
Distinct from apoptosis, ferroptosis is a newly defined form of programmed cell death characterized by irondependent peroxide lipid accumulation, inducing reactive oxygen species (ROS) production and subsequent cell death [4]. Emerging evidence has demonstrated that ferroptosis plays a critical role in the redox status, cell metabolism, and multiple diseases, such as ischemia-reperfusion injury, neurodegenerative and neuropsychiatric diseases, and diverse kidney diseases [4][5][6]. More importantly, ferroptosis dysfunction has been implicated in the process of various tumors, including glioma [7], lung cancer [8], breast cancer [9], renal cell carcinoma [10], colorectal cancer [11], and gastric cancer [12]. Ferroptosis is inhibited in various tumors, resulting in uncontrolled proliferation of tumor cells, and is also involved in immunotherapy and drug sensitivity. Therefore, ferroptosis can serve as a promising interventional target to induce tumor cell death. Recent studies also identified some key ferroptosis-related gene signature, such as glutathione peroxidase 4 (GPX4), solute carrier family 7 member 11 (SLC7A11), nuclear respiratory factor 2 (NRF2), and cysteine dioxygenase type 1 (CDO1), which are closely related to cancer progression and patients' prognosis [13][14][15][16]. However, only few studies to date focused on the role of ferroptosis in GC, and whether ferroptosis-related genes are related to patients' prognosis and clinical treatments still needs to be fully elucidated.
In this study, we constructed and validated a ferroptosisrelated gene signature, which could well predict the prognosis of GC patients. We further performed pathway and functional enrichment analysis to study the underlying mechanisms. The clinical value of the risk model based on this ferroptosis-related gene signature was also explored in immune microenvironment and tumor mutation burden of GC. In addition, the effects of ferroptosis inducer Erastin on these 10 ferroptosis-related genes in GC cell lines were also explored in our study. Similarly, the RNA-seq data of 174 normal human stomach samples in the Genotype-Tissue Expression (GTEx) database was downloaded from the University of California Santa Cruz (UCSC, https://xenabrowser.net/datapages/). Besides, the gene expression data and corresponding clinical information of the external validation cohort (GSE84437, n = 433; GSE29272, n = 268; normal: 134, tumor: 134) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/). Furthermore, the somatic mutation data of the TCGA-STAD was downloaded from the websites (https://portal.gdc.cancer.gov/). The 261 ferroptosis-related genes were downloaded from the FerrDb website (http://www.zhounan.org/ferrdb/), updating on 10 March 2021 [17].

2.2.
Screening of Candidate Gene. The RNA-seq data of TCGA and GTEx datasets was normalized into the transcripts per million (TPM) data. And the scale function in the dplyr R package was employed to further normalize the RNA-seq data (TPM normalized). Batch correction was performed using the sva R package. Then, the RNA-seq data of the 261 ferroptosis-related genes was extracted to perform subsequent difference analysis. Differentially expressed gene (DEG) analysis between the normal and tumor tissues was performed by limma R package, screening out the ferroptosis-related differentially expressed genes (FDEGs). And the results of the FDEGs were visualized by ggplot2 R package. Then, univariate Cox regression analysis was performed to screen out the overall survival-(OS-) associated FDEGs which were identified as the candidate genes for subsequent establishment of prognostic ferroptosis-related gene signature.

Establishment and Validation of a Prognostic
Ferroptosis-Related Gene Signature. In order to minimize the risk of overfitting, the least absolute shrinkage and selection operator (LASSO) regression was utilized to establish a gene prognostic signature and further screen the 10 potential hub genes from the FDEGs by the glmnet R package. Then, the protein-protein interaction (PPI) network was conducted to reveal the interaction of proteins among the protein coding between the 10 genes by the STRING database (http://www.string-db.org/). In order to explore the connection of the transcriptional level among these 10 candidate genes, the igraph and reshape2 R packages were utilized to construct correlation network of these 10 candidate genes. The multivariate Cox regression analysis based on these 10 genes was utilized to establish the prognostic ferroptosisrelated gene signature. The regression coefficients of genes and their corresponding mRNA expressions were utilized to calculate the risk scores of patients. The formula of risk score was established as follows: score = ∑ðcorresponding mRNA expressions × regression coefficientsÞ. The median value of risk scores was utilized as the cutoff value to divide GC patients into the high-and low-risk subgroups. To test the distribution of different groups, the principal component analysis (PCA) was performed by the Rtsne and ggplot2 R packages. The survival curves were performed to analyze the prognostic status between the high-and low-risk groups by the survminer R package. And the time-dependent receiver operating characteristic curve (ROC) was performed to evaluate the predictive value of the prognostic signature by the survival and timeROC R package. Besides, the univariate and multivariate Cox regression analysis was also performed to evaluate the independent prognostic value of the prognostic signature. Furthermore, nomograms of the training and testing groups were constructed to predict the survival probability of GC patients in 1, 2, and 3 years, and their corresponding nomogram calibration curves were also constructed based on the multivariate Cox regression analysis by the rms R package.

Functional Enrichment
Analysis. The OmicShare tools, a free online platform for data analysis (https://www.omicshare .com/tools), was employed to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis with P < 0:05 and normalized enrichment score > 1 based on the DEGs between the high-and low-risk groups. The Gene Set Enrichment Analysis (GSEA) software (https://www.gsea-msigdb.org/gsea/login.jsp/) was also utilized to further reveal the significantly enriched pathways 2 Oxidative Medicine and Cellular Longevity of these DEGs. Furthermore, the maftools R package was utilized to explore and visualize the MAF files of somatic mutation data and also calculate the tumor mutation burden (TMB) scores of patients in the training group.

Immunotherapy Targets and Immune Infiltration
Analysis. To explore the potential relationship between the 10-gene signature and the immune cell infiltration, Tumor Immune Estimation Resource (TIMER) which is a platform for analyzing the abundance of the six immune infiltration cells (macrophages, dendritic cells, neutrophils, CD4 + T cells, CD8 + T cells, and B cells) in malignant tumors was applied to evaluate the associations between the 10 hub genes and the infiltrating immune cells by Pearson correlation analysis and Student's t-test. And the estimate R package was utilized to explore the relationship between the risk scores and immune cell infiltration. Besides, considering the significant roles of more immune cells in the tumor microenvironment, the abundance of 22 infiltrating immune cell types in each tumor sample was calculated by CIBER-SORT. Furthermore, the expression of the target genes which have been reported to be related to the immunotherapy was compared between different risk groups. 2.12. Statistical Analysis. All the statistical analysis was conducted by the R software (version: 3.6.3) in this article. All P 3 Oxidative Medicine and Cellular Longevity values of statistical data were based on two-sided statistical tests. P < 0:05 was considered to be statistically significant.

Establishment of Ferroptosis-Related Prognostic
Signature. 36 samples without complete OS or OS time information in TCGA training group and 136 samples whose OS time over 9 years were eliminated. Then, 29 ferroptosis-related genes were correlated with OS by the univariate Cox analysis (P < 0:05) in the training group and 18 of them were differentially expressed between normal and tumor tissues (Figures 2(a) and 2(b)). LASSO Cox regression analysis was utilized to construct a prognostic signature using the expression value of the 18 prognostic FDEGs mentioned above. Then, a 10-gene signature (SP1, MYB, ALDH3A2, KEAP1, AIFM2, ITGB4, TGFBR1, MAP1LC3B, NOX4, and ZFP36) was filtered out by the minimum value of lambda (λ) (Figure 2(c)). The coefficients of these genes are shown in Figure 2(d). And the full names, function, and coefficients of these 10 genes are shown in Table S1. According to the value of coefficients and hazard ratio (HR), the genes TGFBR1, MAP1LC3B, NOX4, and ZFP36 were considered as the risk genes, while the genes SP1, MYB, ALDH3A2, KEAP1, AIFM2, and ITGB4 as the protective genes. The protein interaction network among these 10 genes indicated that NOX4, SP1, and KEAP1 were the hub genes ( Figure 2(e)). The gene correlation among them is shown in Figure 2(f). Besides, the risk scores of the signature were applied to predict prognosis in GC patients and median risk score was utilized to classify patients into the high-or low-risk groups, which was calculated as follows: risk score = ð−0:181Þ × expression of SP1 + ð−0:085Þ × expression of MYB + ð−0:076Þ × expression of ALDH3A2 + ð−0:075Þ × expression of KEAP1 + ð−0:031Þ × expression of AIFM2 + ð −0:026Þ × expression of ITGB4 + ð0:072Þ × expression of T GFBR1 + ð0:138Þ × expression of MAP1LC3B + ð0:148Þ × expression of NOX4 + ð0:345Þ × expression of ZFP36.

Evaluation and Validation of Ferroptosis-Related Gene
Signature. As shown in Figure S1, PCA of the training and testing groups revealed that the patients in different risk groups could be distributed in two discrete directions. The high-risk GC patients were more likely to die earlier than low-risk patients from the results of the scatterplots (Figure 3(a)) and heat maps (Figure 3(b)). Besides, the Kaplan-Meier survival curve (Figure 3(c)) indicated that patients with low-risk scores may have a better prognosis than patients with high-risk scores in the training group. The sensitivity and specificity of the risk scores to predict prognostic features were determined from the timedependent ROC curves by calculating the areas under the curve (AUC). And the risk scores presented the potential ability of predicting the OS status (1-year AUC = 0:722, 2year AUC = 0:704, and 3-year AUC = 0:680) in Figure 3(d).
In order to avoid the contingency of TCGA results, patients in the testing group were also classified into the low-and high-risk groups based on the median risk score. Similar to all the results above, patients with high-risk scores had a higher probability to encounter death earlier and had worse overall survival outcome than those with low-risk scores (Figures 3(e)-3(h)).

Construction and Validation of the Nomogram Prediction Model.
In order to predict the survival probability of GC patients at 1, 2, and 3 years, the clinicopathological characteristics, including grade, N stage, T stage, TNM stage, and risk score, were applied to construct the nomogram prediction model in the training group (Figure 4(e)). The corresponding calibration curves were shown to perform the good prediction in the observations in 1-3 years ( Figure 4(f)). Thus, the nomogram incorporating clinical features and the risk score was stable and accurate and may be applied in the clinical evaluation of GC patients.
3.6. Analysis of Functional Enrichment. In order to further explore the signature-related downstream molecular biological functions and pathways, the GO enrichment and KEGG pathway analyses were performed by the DEGs between the      Besides, tumor mutation burden (TMB) was calculated and analyzed in both groups, indicating that TMB level was significantly higher in the low-risk group (Figure 6(f)).

Analysis of Tumor Microenvironment and
Immunotherapy Response. According to the results of the functional enrichment analysis, immune process was significantly different between the high-and low-risk GC patients ( Figure S2). Thus, TIMER and CIBERSORT analysis was utilized to further explore the relationship between risk score and tumor microenvironment. Firstly, the results of TIMER analysis indicated that the 10 FDEGs were associated with all 6 immune infiltration cells (purity, B      11 Oxidative Medicine and Cellular Longevity cell, CD8 + T cell, CD4 + T cell, macrophage, neutrophil, and dendritic cell), especially for NOX4, AIFM2, and SP1 genes ( Figure S3). Meanwhile, CIBERSORT was also applied to estimate the different infiltration abundance of 22 immune cells between the high-and low-risk groups in the training group. The results showed that mast cells resting, B cells naive, dendritic cells resting, and monocytes were downregulated in the low-risk groups, while NK cells resting, macrophages M0, and T cells follicular helper were significantly upregulated (P < 0:05, Figures 7(a)-7(c)). Besides, the correlation analysis of risk score with common immune checkpoints (ICPs), including cytotoxic T lymphocyte-associated protein 4 (CTLA4), programmed cell death 1 (PDCD1) (PD1), CD274 (PD-L1), hepatitis A virus cellular receptor 2 (HAVCR2), and lymphocyteactivating 3 (LAG3), was performed to estimate the immunotherapy responses through the 10-gene signature. As expected, the gene expression levels of most ICPs were significantly upregulated in the high-risk group (Figure 7(d)).        16 Oxidative Medicine and Cellular Longevity levels of SP1, MYB, KEAP1, AIFM2, ITGB4, TGFBR1, MAP1LC3B, and NOX4 were significantly upregulated, while the expression of ALDH3A2 and ZFP36 was downregulated in GC tumor tissues in the training group (Figure 8(a)). The same results were also verified by GEPIA (Figure 8(b)). Similarly, compared to the normal gastric cell GES1, most of these 10 genes were also differentially expressed in GC cell lines (HGC-27 and MGC-803) using real-time PCR (Figure 8(c)). Besides, we further validated the mRNA or protein expression of these 10 genes in GSE29272 and Human Protein Atlas (HPA) datasets in Figure S4. Furthermore, the expression levels of hub genes SP1, KEAP1, AIFM2, and NOX4 were further verified in our GC samples. The results showed that the expression of SP1, KEAP1, AIFM2, and NOX4 all increased in tumor tissues (Figure 8(d)). 19 Oxidative Medicine and Cellular Longevity real-time PCR. The results indicated that the mRNA expression levels of eight genes (AIFM2, ALDH3A2, KEAP1, MAP1LC3B, MYB, NOX4, SP1, and TGFBR1) were decreased and two genes (ITGB4 and ZFP36) were increased after being treated with Erastin. However, in the HGC-27 cell line, there was no statistical difference in the mRNA expression level of SP1. ITGB4 and MAP1LC3B genes were also not statistically different in MGC803-Erastin cell line (Figures 9(g) and 9(h)). In addition, similar to the mRNA expression results, the different protein expression levels of the hub genes (AIFM2, KEAP1, NOX4, and SP1) were further confirmed by western blot except AIFM2 and SP1 in the HGC-27 cell line (Figure 9(i)). In summary, the potential roles of these 10 ferroptosis-related gene markers could also be verified in cell line experiment.

Discussion
In this study, the expression level of the ferroptosis-related genes in GC tumor and normal tissues and their associations      23 Oxidative Medicine and Cellular Longevity with OS were systematically investigated. A novel prognostic gene signature was established and validated in an external cohort. The independent prognostic factor, functional enrichment, somatic mutation, tumor microenvironment, and immunotherapy response analysis were performed and indicated that the ferroptosis-related gene signature can effectively predict the prognosis and clinical status for GC patients.
Ferroptosis is involved in various diseases, especially in malignant tumors [19]. Recently, several studies [12,[20][21][22] have proven that some ferroptosis-related genes play key roles in the process of tumorigenesis and progression of GC, but whether ferroptosis could predict the prognosis and clinical status of GC patients remains largely unknown. Usually, TNM stage system or some serum biomarkers, including CEA, CA19-9, and CA125, are used to monitor the progress and predict the prognosis of GC patients. However, these approaches are not satisfactory with low accuracy and high nonspecificity; especially, there is higher heterogeneity in GC patients. Meanwhile, with the impressive progress of bioinformatics and RNA-seq, many scholars around the

24
Oxidative Medicine and Cellular Longevity world have constructed some ferroptosis-related gene signature by public databases to further explore key molecular markers and better methods to accurately predict the prognosis and drug sensitivity in several malignant tumors, including uveal melanoma, lung cancer, hepatocellular carcinoma, pancreatic cancer, and glioma [23][24][25][26]. However, few studies on gene signature had been constructed in gastric cancer.
In this study, we first screened the key ferroptosis-related DEGs in GC from the public databases. As expected, more than half of the ferroptosis-related genes were differentially expressed between adjacent nontumorous and tumor tissues in GC patients, suggesting ferroptosis plays a significant role in GC. Then, 29 of them were proven related to OS by the univariate Cox analysis, indicating that constructing a prognostic signature with these FDEGs is feasible and reasonable. Using the LASSO Cox analysis, the novel prognostic signa-ture integrating 10 ferroptosis-related genes was identified, including SP1, MYB, ALDH3A2, KEAP1, AIFM2, ITGB4, TGFBR1, MAP1LC3B, NOX4, and ZFP36. To further explore the role of these genes in GC, we summarized their main molecular functions based on the results of this study and previous studies. SP1, a key member of the transcription factor SP family, plays important roles in tissue development, cell differentiation, and tumor molecular biology [27]. SP1 can directly positively regulate glutathione peroxidase 4 (GPX4), which is able to significantly influence the level of lipid peroxidation and inhibit ferroptosis [28]. AIFM2 belongs to the antiferroptotic genes and is renamed as ferroptosis suppressor protein 1 (FSP1). Recent studies indicated that AIFM2 plays a significant role in ferroptosis and can act parallel to GPX4 to inhibit ferroptosis [29,30]. ALDH3A2 is involved in preventing cellular oxidative damage by oxidizing long-chain 25 Oxidative Medicine and Cellular Longevity aliphatic aldehydes. A recent study showed that ALDH3A2 can protect progenitor and leukemic stem cells from ferroptosis. Inhibiting GPX4 expression can further enhance the ferroptosis-inducing influence of ALDH3A2 depletion [31]. Different from other types of tumors, the expression of GPX4 is positively correlated with the prognosis of GC patients [12]. Thus, the potential role of SP1 and FSP1 during the process of ferroptosis in GC patients remains to be further explored. MYB has also been reported as an important transcription factor in solid tumors, which can regulate ferroportin expression, iron-related cellular activities, and tumor cell growth by modulating myeloid zinc-finger 1 [32]. Notably, a specific study indicated that MYB could inhibit Erastin-induced ferroptosis which was restrained through interacting with CDO1 in GC cells [13]. KEAP1 interacts with nuclear factor erythroid 2-related factor 2 (NRF2) in a redox-sensitive manner, and the interaction can promote the expression of gamma-glutamylcysteine synthetase [33]. In recent study, the NRF2-KEAP1 pathway is activated and upregulates SLC7A11 to inhibit ferroptosis when the expression of KEAP1 is downregulated [34]. ITGB4, a member of the integrin family, mediates cell-cell adhesion or cell growth and plays a significant role in the biology of invasive carcinoma by associating with integrin alpha 6 (ITGA6) subunit [35]. Besides, it has been reported that the induction of ferroptosis depends on cell clustering in matrix-detached cells that lack ITGB4 and ITGA6 expression [36]. TGFBR1, also known as the activin receptor-like kinase (ALK4/5), is involved in oxidative stress responses [37]. In renal proximal tubular epithelial cells, the ALK4/5 signaling pathway has been proven to be correlated with ferroptosis and blockade of the ALK4/5 signaling pathway can suppress ferroptosis [38]. Recent evidences demonstrated that autophagy facilitates ferroptosis by degrading antiferroptosis factors [39]. MAP1LC3B [40] and ZFP36 [41], key proteins of autophagy, have been considered to be correlated with ferroptosis. NOX4 is the core enzyme in mediating lipid peroxidation and promoting ferroptosis, and inhibition of NOX4 can significantly block ferroptosis [42,43]. In summary, although these 10 genes were all correlated with ferroptosis, few studies were performed to explore their molecular functions during the process of ferroptosis in GC. Thus, to make these results more scientific, we also choose the 4 hub genes (SP1, KEAP1, AIFM2, and NOX4) of these 10 genes to further verify their expression levels in the GC cell lines and our 30-paired GC tissues by real-time PCR and immunohistochemistry. In summary, similar to the gene expression level of public databases and related studies, the expression level of them was also upregulated in GC cell lines and tissues.
Based on risk score of the 10-gene signature, GC patients can be divided into the low-and high-risk groups. Low-risk GC patients were proven to have the better prognosis and significantly longer OS than high-risk patients in both the training and testing groups. Furthermore, a series of analysis was applied to further explore the prognostic value of the signature; the results showed that the risk score was the independent prognostic factor of the OS in GC patients. Accurate nomogram prediction models can also be con-structed based on the risk score. In summary, these results reveal a favorable predictive efficacy of the signature in both the training and testing groups. Meanwhile, we performed GO, KEGG, and GSEA analysis to identify the enriched biological process and pathway based on the DEGs in the highand low-risk patients. The results showed that the cell cycle, ECM-receptor interaction, PI3K-Akt signaling pathway, and tumorigenesis-related pathways were significantly enriched both in the training and testing groups. Consistent with our results, recently, Lin et al. [44] reported that dihydroartemisinin can cause cell cycle arrest in head and neck carcinoma (HNC) cells by inducing ferroptosis. Some studies [45,46] also demonstrated that epigenetic reprogramming of EMT promotes ferroptosis in HNC cells and gambogenic acid-induced ferroptosis inhibits the EMT in melanoma cells. Yi et al. [47] found that mutation of PI3K-Akt signaling could protect cancer cells from oxidative stress and ferroptosis.
Recent studies [48,49] have proven that somatic mutation and tumor immune microenvironment significantly correlate with tumorigenesis, tumor progress, and drug resistance in GC patients. Wang et al. demonstrated that IFNγ released by CD8 + T cells could inhibit expression of glutamate-cystine antiporter system xc -, then induce tumor cell lipid peroxidation and ferroptosis, and finally improve antitumor efficacy of immunotherapy [50]. Besides, Hung et al. reported that tyrosine-protein kinase receptor TYRO3 (TYRO3) overexpression elicited anti-PD-1/PD-L1 resistance through protecting tumor cells from immunotherapyinduced ferroptosis [51]. However, the specific mechanisms of ferroptosis in tumor immunotherapy are largely unknown. In order to explore the potential mechanism of this signature, we further performed somatic mutation, tumor microenvironment, and immunotherapy response analysis. To our surprise, the somatic mutation frequency of the low-risk GC patients was higher than that of the high-risk patients and the TMB level was significantly higher in the low-risk group indicating that low-risk GC patients may be more sensitive to immunotherapy and can benefit from the immunotherapy. Based on the TMB results, we also confirmed that there were significant differences in the signature and immune checkpoints between the high-and low-risk groups. Low-risk patients had a better response to immunotherapy, suggesting that the signature has the potential to predict the immunotherapy response in GC patients. Meanwhile, the infiltration abundance of immune cells was significantly different between the low-and high-risk groups in this study. Most of the immune cells were highly infiltrated in the high-risk patients, while the abundance of NK cells was higher in the low-risk group. Previous reports indicated that increased abundance of NK cells can bring better immunotherapy efficacy [52], further suggesting the low-risk GC patients have a better response to immunotherapy. In fact, the relationship between ferroptosis and cancer immunotherapy has been reported in 2019 [50,53], showing the sensitivity of tumor cells to ferroptosis is parallel to immune functions, which may be the reasonable explanation for the better response to immunotherapy in low-risk patients. But the exact mechanisms of how these 10 ferroptosis-related genes interact to 26 Oxidative Medicine and Cellular Longevity affect tumorigenesis and immune process are still unclear and further studies are demanded.
With the gradual understanding of the role of ferroptosis in various cancers, different ferroptosis inducers have been developed as anticancer therapies, such as sulfasalazine, sorafenib, and Erastin [54,55]. Erastin, first identified as killing tumor cells expressing oncogenic RAS, is a classical agent to induce ferroptosis by suppressing cystine/glutamate antiporter (xCT) and leading to decreased cysteine and then inhibiting the function of glutathione peroxidase 4 (GPX4) [4]. To evaluate whether these key genes play key roles in ferroptosis of GC, we applied Erastin to trigger ferroptosis in two GC cell lines, HGC-27 and MGC-803. Using CCK-8 assay, we detected significant lethal toxicity to both GC cell lines by Erastin at low concentration, showing IC 50 value was 7.46 and 10.79 μM, for HGC-27 and MGC-803, respectively. Increased ROS level is one of the features for ferroptosis. After treatment with Erastin, both cell lines showed obvious increased ROS signal by flow cytometry and fluorescence microscope, indicating Erastin-induced cell death could be attributed to ferroptosis. Besides, we also explored whether the expression of these key ferroptosis-related genes was regulated by Erastin during ferroptosis. Not surprisingly, most of the genes were dysregulated at mRNA or protein expression level, while the underlying mechanism was still not clear. Therefore, these ferroptosis-related gene could be targeted to induce cancer cell ferroptosis for future personalized therapy.

Conclusion
Ferroptosis has great potential clinical value in tumor treatments. However, the relationship between ferroptosis and tumors such as GC remains largely unclear. Thus, we systemically explored the key ferroptosis biomarkers in GC and constructed a novel ferroptosis-related gene signature which could effectively predict the prognosis of GC patients. Besides, we also further explored the signature-related downstream molecular biology functions and pathways, demonstrating the potential clinical value of this signature in somatic mutation and immunotherapy. In addition, all results had been verified by various external datasets and expression of 4 hub genes was verified in our own clinical samples. Finally, the novel prognostic signature constructed in this study needs further validation and the underlying mechanisms of ferroptosis in GC should be explored in future studies.

Data Availability
The RNA sequencing (RNA-seq) data and corresponding clinical characteristics and molecular information of gastric cancer samples in training cohort (TCGA-STAD) were downloaded from The Cancer Genome Atlas (TCGA) database by the "TCGAbiolinks" R package in February 2021. The RNA-seq data of normal human stomach samples in GTEx database was downloaded from the University of California Santa Cruz (UCSC, https://xenabrowser.net/ datapages/). Besides, the gene expression data and corre-sponding clinical information of the external validation cohorts (GSE84437, GSE29272) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/). Furthermore, the somatic mutation data of the TCGA-STAD was downloaded from the websites (https://portal.gdc.cancer .gov/). The 261 ferroptosis-related genes were downloaded from the FerrDb website (http://www.zhounan.org/ferrdb/), updating on 10 March 2021.

Ethical Approval
The study was conducted according to the guidelines of the Declaration of Helsinki, approved by the Ethical Committee of Ruijin Hospital (Shanghai, China). Figure S1: PCA plots of the TCGA-STAD training and GSE84437 testing datasets. Figure S2: KEGG circular and pathway annotation plots of the TCGA-STAD training and GSE84437 testing datasets. Figure S3: the diagrams of the correlation analysis between these 10 FDEGs and the immune infiltration level in TCGA dataset by TIMER. Figure S4: validation of the mRNA or protein expression of these 10 genes in GSE29272 (except AIFM2) and HPA (except NOX4) datasets.