Identification of Ferroptosis-Related Genes as Biomarkers for Sarcoma

Sarcomas are seen as mixed-up nature with genetic and transcriptional heterogeneity and poor prognosis. Although the genes involved in ferroptosis are still unclear, iron loss is considered to be the core of glioblastoma, tumor progression, and tumor microenvironment. Here, we developed and tested the prognosis of SARC, which is a genetic marker associated with iron residues. The ferroptosis-related gene expression, one-way Cox analysis, and least-selection absolute regression algorithm (LASSO) are used to track prognostic-related genes and create risk assessment models. Finally, immune system infiltration and immune control point analysis are used to study the characteristics of the tumor microenvironment related to risk assessment. Moreover, LncRNA–miRNA–mRNA network was contributed in our studies. We determined the biomarker characteristics associated with iron degradation in gene 32 and developed a risk assessment model. ROC analysis showed that its model was accurately predicted, with 1, 2, 3, 4, and 5 years of overall survival in TCGA cohort of SARC patients. A comparative analysis of settings found that overall survival (OS) was lower in the high-risk than that in the low-risk group. The nomogram survival prediction model also helped to predict the OS of SARC patients. The nomogram survival prediction model has strong predictive power for the overall survival of SARC patients in TCGA dataset. GSEA analysis shows that high-risk groups are rich in inflammation, cancer-related symptoms, and pathological processes. High risk is related to immune cell infiltration and immune checkpoint. Our prediction model is based on SARC ferritin-related genes, which may support SARC prediction and provide potential attack points.

INTRODUCTION SARC (sarcoma) is a type of stromal cancer with more than 100 subtypes. Sarcomas occur in all age groups and are more general in teens and young adults than in the elderly (Reed et al., 2019). SARC also includes many malignant tumors with different symptoms, behaviors, and results. For example, chemotherapy has only a known effect on a few sarcomas. Although confirmed sarcoma can be removed by surgery or radiotherapy, the incidence of metastatic sarcoma within 5 years is still as high as 50% (Brennan et al., 2014). Only 5% of metastatic patients have a 5-year survival period .
Although some of the diagnostic markers of SARC can effectively predict the progress of SARC, they are still at the molecular stage and not yet used in clinical practice. Therefore, diagnostic markers related to undiscovered genes are of great significance for the diagnosis and prediction-related analysis of SARC.
Ferroptosis, an iron-dependent form of programmed cell death, which is linked to pathophysiological conditions. Sensitivity to lipid peroxidation is attributed to many metabolic pathways, such as fatty acid metabolism, iron processing, methionine metabolism, and mitochondrial respiration (Zheng and Conrad, 2020;Jiang et al., 2021). Notably, ferroptosis is related to stroke, ischemic heart disease, liver and kidney damage (Benjamin et al., 2017). In particular, there is strong evidence that high levels of ferroptosis play an important role in the suppression of tumor growth, and the activation of iron levels will enhance the therapeutic effect of anticancer drugs . Cancer cells store large amounts of iron and active oxygen to stimulate metabolism and growth (Battaglia et al., 2020). Gene dysregulation involved in iron nucleocytosis may promote cancer cell proliferation, invasion, and metastasis (Jung et al., 2019). In addition, CD8 + T cells cause iron loss by controlling the link between cancer cell death mechanism and immune system activation, thereby inhibiting tumor growth (Tang et al., 2020).
In the past, there was an association between many subtypes of sarcoma and apoptosis-related genes such as GPX4, which induces iron loss by increasing MDA, ROS, and intracellular iron levels in osteosarcoma cells (Lin H et al., 2021). LOX promotes tumor deposition and catalyzes the production of lipid hydrogen peroxide in the plasma membrane. In the process of ptosis, ALOX15 protein is always present in the cell membrane (Shintoku et al., 2017). Pharmacological lipid peroxidation inhibitors (such as ferrostatin-1 and lipid peroxidin-1), reactive oxygen species scavengers (such as alpha-tocopherol and glutathione), and iron chelator deferoxamine can inhibit the accumulation of reactive oxygen species; lipid peroxidation oxidizes and heals swamps (Dächert et al., 2020). However, the side effects of ferroptosis limit their applications in treating sarcomas (Blatt, 1994). Therefore, sarcoma iron gout needs more understanding.
These studies focused on genes related to iron-related diseases. We perform extensive bioinformatics analyses based on The Cancer Genome Atlas (TCGA) clinical data to analyze the gene expression levels, DNA methylation, and copy number transfer models. The SARC risk assessment system in TCGA dataset was developed and validated by detecting regulated iron-related genes. In addition, the function and gene clusters were performed to identify possible pathways and mechanisms of ferroptosis metabolism. Our results indicate that four genetic markers can be used as independent predictors of OS in sarcoma patients.

Analysis of the Ferroptosis-Related Genes in SARC
Ferroptosis-related genes are obtained from the genset (the Molecular Feature Database (MSigDB) version 7.1) (Subramanian et al., 2005;Liberzon et al., 2015), including the 15 gene sets (Mou et al., 2020;Zhang S et al., 2020). After removing overlapping genes, we got a set of genes related to ferroptosis-related gene (FRG), including 32 genes.

Datasets and Data Processing
As of 30 November 2021, we have obtained data on RNA sequences and summary of 260 SARC patients' clinical characteristics from TCGA database. The RNA sequence data of 88 normal human ovarian samples were downloaded from the GTEx database. The RNA-seq data and clinical information from the external validation cohort were extracted from the GEO database (GSE21050).
A total of 32 FRGs were obtained from previous studies (Lin W et al., 2021;Ye et al., 2021). The "limma" packages demonstrate the difference between FRG expression in SARC and normal tissues. Then, we contribute 32 protein-protein interaction networks (PPI) to search for interaction finder (STRING).

Gene Mutation, Methylation, and Copy Number Variation
The "maftool" package is used to create single-nucleotide variation, copy number variation, and 32 FRG cascade plots for SARC patients. GSCA Lite is used to keep methylation alive and analyze the changes in copy number.

Gene Functional Enrichment Analysis
Genetic ontology (GO), which includes biological process (BP), cell composition (CC), and molecular function classification (MF), uses the "ggplot2" software package in R. In the same way, Kyoto Encyclopedia and Genome uses this software package to analyze (KEGG) genes. A hypergeometric distribution test was applied to detect enrichment terms, and p values were adjusted by the false discovery rate (FDR) method with a cutoff FDR <0.05.

Development of Ferroptosis-Related Prognostic (FRG) Model
Cox regression analysis (COX analysis) is used to investigate the effects of FRGs on prognosis of SARC. The Kaplan-Mayer curve uses the logarithmic test and one-way regression of Cox proportional hazards to generate p (HR) values with 95% confidence intervals (CIs). PRGs with important prognostic value were selected in further analysis. Based on these prognostic-related PRGs, LASSO-Cox regression analysis was used to contribute a prognostic model. According to the average risk score, TCGA-SARC patients were divided into low-risk and high-risk subgroups, and Kaplan-Meier analysis was used to compare the overall survival rates of the two subgroups. Time receiver performance (ROC) analysis is used to evaluate the prediction accuracy and risk assessment of each gene. Taking into account the clinical characteristics, we developed a nomosis prediction schedule to predict the overall survival rate at 1, 3, and 5 years. A forest was used to show the p-value, HR, and 95% CI of each variable through the "forest plot" R package.
To evaluate whether the risk score system can serve as an independent predictive index, univariate and multivariate Cox regression analyses were performed with clinicopathologic parameters. All independent prognostic parameters were used to construct a nomogram to predict the 1-, 3-and 5-year OS probabilities by the "rms" package. The concordance index (C-index), calibration, and ROC analyses were used to evaluate the discriminative ability of the nomogram (Alba et al., 2017).

Immune Related Analysis in Ferroptosis-Related Prognostic (FRG) Model
We used the tumor immune assessment function to analyze the relationship between prognostic FRGs and immune system infiltration, and comprehensively analyze the immune cells infiltrating the portal. The TIM gene module can monitor the relation between the expression of SARC gene and the degree of immune infiltration. Spearman correlation analysis was used to calculate the correlation between gene expression, and the results of tumor mutation burden (TMB) and microsatellite instability (MSI). p values less than 0.05 are considered statistically significant.

GSEA
The aforementioned software package R was used to calculate the DEG training package. Then GSEA (http://software.broadinstitute. org/gsea/index.jsp) was used to identify the characteristics of highand low-risk groups.

TIMER Database Analysis
The result of the infiltration estimation generated by the TIMER algorithm consists of 6 specific subgroups of immune cells, including B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and dendritic cells (Sturm et al., 2019;. We extracted the results of the penetration assessment and evaluated the different results of the penetration assessment of immune cell subsets between high-risk and low-risk populations (Newman et al., 2015).

Statistical Analysis
All statistical analyses in this study were performed using software R (version 3.6.3) and GraphPad Prism (version 8.0.2). Kaplan-Mayer survival analysis uses the logarithmic series test. Risk factors (HR) and 95% confidence intervals (CI) should be provided when relevant. Both groups were compared with Student's t-test and the Kruskal-Wallis test. A two-tailed p value less than 0.05 was considered statistically significant without specific annotation.

Expression of FRGs in SARC
First, we use TCGA data to study the expression of 32 FRGs in SARC and normal tissues. The total amount of SARC is 4 FRGs. For example, the expression of SLC7A11, FANCD2, and CISD1 Frontiers in Cell and Developmental Biology | www.frontiersin.org March 2022 | Volume 10 | Article 847513 4 was higher in tumor tissues than that in normal tissues, but ATP5MC3 was lower. A GeneMANIA network showed that FANCD2, SLC7A11, CISD1, and ATP5MC3 were performed to detect these FRG interactions. The minimum required interaction score is 0.9. The heat map of results showed that FANCD2, SLC7A11, CISD1, and ATP5MC3 are key genes ( Figure 1, SupplementaryFigure S1, Supplementary Table S1).

Landscape of Genetic Variation, Methylation, and Copy Number Variation
We first analyzed the relationship between FRG and methylation. In all FRGs, we found that the methylation rate of SLC1A5 and ATL1 is associated with a methylation level ( Figure 1B). Then we analyzed the percentage plot of individual nucleotide variation on the chart. The ACSL4 mutation frequency is high ( Figure 1C). We summarized the copy number and frequency of 32 FRG somatic mutations in SARC. As shown in Figures 1D,E, we found that CISD1 and NCOA4 have a higher level of CNV pie distribution, and CS, CARS, FANCD2, and LPCAT3 have a higher CNV expression level. ALOX15 has a high circulating distribution of CVN (copy number change), while CISD1, ACSL4, and ALOX15 show high levels of homozygous SARC deletion. Supplementary Figure  S5C shows that TP53 and RB1 have higher levels of genomic alteration. The mutation count, TMB, fraction genome altered, aneuploidy score, MSIsensor score, radiation therapy, and Onc Tree code, and cancer type details also show significant difference in genome changes of FRGs (Supplementary  Table S3).

Functional Enrichment Analysis
To illustrate the role of FRGs, we uses GO and KEGG databases to analyze the FRGs pathway. We found that these 32 FRGs are mainly involved in the p53 transcriptional gene network, response to toxic substance, ER-nucleus signaling pathway, endoplasmic reticulum organization, cellular transition metal ion homeostasis, and peptide transport are involved in the process of ferroptosis in GO assay ( Figure 2A). For COVID related GO enrichment analysis, RNA lamers intestinal organoid expansion, RNA Sun calu-3 24 h, proteome stukalov, translatome luarant also play an important role in the ferroptosis related pathway. The cell type signature for example response to stimulus, signaling and negative regulation of biological process and transcription factor target such as lake adult kidney c19 collecting duct intercalated cells types a medulla also play potential role in ferroptisis related pathway. In

Ferroptosis-Related Prognostic Gene Model
One-way Cox regression analysis was used to select the prognostic FRGs for the genetic prognostic model. As a result, a total of 4 genes with prognostic values were identified. Kaplan-Mayer survival curves are shown in Figures 3A,B. For overall survival, we found that the high expression of SLC7A11(p = 0.023) and FANCD2 (p = 0.011) was correlated with good prognosis, but that of NCOA4 was linked with poor prognosis. For disease-free survival, we found that high expression of CISD1 (p = 0.0042) links with better prognosis ( Figure 2C). Based on these 4 prognostic FRGs, LASSO-Cox regression analysis was performed to construct a prognostic gene model (Figure 4). According to risk assessment, SARC patients are divided into two groups. The risk results, survival status, and expression distribution of these genes are shown in Figure 4. The higher the risk score, the higher the risk of death and the shorter the survival time ( Figure 4C). The Kaplan-Meier curve shows that the overall survival rate of high-risk patients is not only higher than that of low-risk patients (average time = 4.2 and 6.6 years, p = 0.000228, Figure 4D) but also 1 year, 2, 3, 4 year annual, and 5-year ROC curves are 0.572, 0.628, 0.609, 0.56, and 0.596, respectively.

Building a Predictive Nomogram
Considering the clinicopathologic and prognostic characteristics of FRGs, we also compiled the prognostic nomogram of the survival prediction table. One-dimensional and multidimensional analyses showed that NOD1 and pT staging, pN staging, and pM staging are independent factors affecting the prognosis of SARC patients ( Figure 5). The chart with predictable names shows that the 1, 2, 3, and 5-year overall survival of the entire cohort can be better predicted than the ideal model (index C: 0.625 (0.536-1), p< 0.01). Immune Infiltration and Immune Survival in SARC Immunohistochemical results suggested SLC7A11 and FANCD2 staining in melanoma ( Figure 6A). Ferroptosis plays an important role in the creation of tumor immune microenvironment. In our study, we also used the TIMER database to clarify the relationship between the expression of prognostic FRG (FANCD2, SLC7A11, CISD1, and ATP5MC3) and SARC immune filtration and somatic copy number variation. For example, B cells in FANCD2 are significantly expanded, lacking CD4 + T cells and macrophage branching. However, FANCD2 have no significant difference in SLC7A11, CISD1, and ATP5MC3 ( Figure 6). SLC7A11 was negatively correlated with CD4+ T cells, and FAND2 was negatively correlated with CD4+ T cells and macrophages. The data show that FANCD2 and CISD1 are positively correlated with B cells. CISD1 was also positively correlated with CD8+ T cells but negatively correlated with CD4+ T cells. ATP5G3 was negatively correlated with macrophages ( Figure 7A). Immune survival ( Figure 7B) shows that lower CD4 + T cells and neutrophils level also have better prognosis. Moreover, a lower FANCD2 level also predicted poor prognosis.

Drug Sensitivity, TMB, and MSI in SARC
The instability of TMB and microsatellites can be used as biomarkers to predict the effect of tumor immunotherapy (Campbell et al., 2017;Latham et al., 2019;Samstein et al., 2019;Chan et al., 2020). The aforementioned analysis shows that FRGs are negatively correlated with tumor immune system infiltration. To determine whether this FRG can also be used as a biomarker for drug screening, we analyzed the correlation between FRG and TMB and MSI in SARC. The results show that MSI is positively correlated with SLC7A11 (p = 0.010) and FANCD2 (p = 7.06E-5). In addition, the TMB also negatively correlated with FANCD1, SLC7A11, CISD1, and ATP5MC3. SLC7A11 has negatively correlated with ESTIMATE score (p = 7.45E-7), immune score (p = 3.7E-5), and Stromal Score (p = 5.16E-8). FANCD2 was negatively correlated with Stromal Score (p = 0.018). CISD1 negatively correlated with ESTIMATE score (p = 0.044) and Stromal Score (p = 4.81E-4). Immune score shows that FANCD2 was correlated with Analyzing the correlation between gene expression and available drugs is key to setting treatment goals. In our study, we analyzed 14 genes that are significantly related to some or most anticancer drugs, such as FANCD2, HSPA5, NFE2L2, ACSL4, and DPP4 in the Cancer Therapeutic Response Portal Database (Supplementary Figure S2).

Contribute an miRNA-LncRNA-mRNA Network
To elucidate the molecular mechanism of FRGs in SARC, we constructed an mRNA-miRNA-lncRNA interaction network. We first investigate the correlation hub gene, and Ven plot found that 101 hub genes have been analyzed, as shown in Supplementary Figure S3A and Table 1. The ferroptosis-related gene pathway analysis found the energy produced by the oxidation of organic compounds and cellular respiration in the BP pathway and mitochondrial inner membrane, mitochondrial protein complex in CC pathway, electron transfer activity, NADH dehydrogenase activity in MF and Parkinson disease, oxidative phosphorylation, and non-alcoholic fatty liver disease in KEGG (Supplementary Figure S3A,B).

DISCUSSION
Ferroptosis is a recently found form of cell-regulated death related to the metabolism pathway in human health (Liang et al., 2019). Accumulation of iron (II) and lipid peroxidation play important roles in causing ferroptosis. Iron chelators and lipophilic antioxidants can inhibit ferroptosis (Su et al., 2020). However, there is increasing evidence that iron poisoning plays an important role in neuro-diseases, including neurodegenerative diseases and cardiovascular diseases . In addition, a thorough FRG landscape analysis of soft tissue sarcoma revealed a new type of FRG related to cancer and prognosis (Huang et al., 2021). Ferroptosis has recently become a treatment target and potential biomarkers for sarcoma .
We performed this study for the first time using gene information of SARC from open databases. First, we chose 32 genes, which are related to ferroptosis. Among them, four ferroptosis genes were analyzed using Cox and LASSO regression to select potential prognostic markers of one variant and then to construct a prognostic model. Among them, the expression level of two genes (SLC7A11 and FANCD2) was negatively correlated with OS, and the expression level of one gene (NCOA4) was positively correlated with OS. These data are consistent with previous results. Sun et al. found that the inhibition of GSH and SLC7A11 is the main cause of EMT and iron deficiency in A549 cells . SLC7A11 is becoming a potential new cancer treatment target. This article briefly introduces the structure and function of SLC7A11 (Lin et al., 2020). FANCD2 has played a new role in reducing ferroptosis. FANCD2 can be used as a target for the development of new cancer therapies aimed at reducing the side effects of ferroptosis (Song et al., 2016). ATP5MC3 is better predictor in a risk prognosis model and also be defined as potential drug for treatment of prostate cancer.
We also analyzed the functional enrichment of these FRGs and found that these 32 FRGs are mainly involved in the regulation of the p53 gene transcription network, response to toxic substance, ER-nucleus signaling pathway, endoplasmic reticulum organization, cellular transition metal ion homeostasis, and peptide transport (Kang et al., 2019). It has been confirmed that the level of glutathione in intestinal tissue is the most influential metabolic pathway related to glutathione/GPx 4 after exposure to microcystin . Colitis mice is evidently induced by ferroptosis, which is mediated by stress signals from the endoplasmic reticulum . The enigmatic lipid peroxidation product is believed to be the direct agent of ferritin metabolism, a special death program caused by glutathione peroxidase 4 (GPX4) insufficiency (Kagan et al., 2017).
In recent years, it has fundamentally changed the way cancer is treated (Dupont et al., 2021), which pointed out that by integrating multi-dimensional biological data and clinical characteristics, highly heterogeneous tumors can be classified into more specific subtypes for personalized treatment . In fact, the accumulation of evidence based on molecular profiling has established itself in the subgroup of cancer patients representing different phenotypes, prognosis, and treatment response (Tan et al., 2019). According to the characteristics of immunogen expression, patients with sarcoma can be divided into two subtypes: prognosis and clinical importance. Patients with high-risk immune subtypes are more sensitive to immune checkpoint blockers (Francescutti and Skitzki, 2012;Tsukahara et al., 2016).
In addition, we developed a ferritin-related gene-based nomogram prediction model for predicting OS in SARC patients. The calibration table and ROC analysis show that the nomogram has a reliable predictive ability for the queue TCGA operating system. The Nomap model can be used to determine the patient's prognosis and make follow-up plans.
GSEA showed that the drug is more sensitive to immune response and tumor progression in the high-risk group. In Kaposi's sarcoma, chemotherapy with vincristine, bleomycin, and etoposide appears to improve overall survival, survival, and quality of life (Macken et al., 2018). Methotrexate also plays an important role high-grade osteosarcoma in children and young adults (van Dalen et al., 2011;Abe et al., 2021). The overall biological response in Ewing's sarcoma exposure shows differences in response (Aryee et al., 2013). miRNA and lncRNA play a role in a variety of biological and pathological processes, such as apoptosis, cell cycle, migration, invasion, regulating proliferation, metastasis, EMT, and resistance to drugs. It plays an important role and has been extensively studied. mRNAs are target through transcription or post-transcription . Therefore, the LncRNA-miRNA-mRNA network was contributed in our studies. Perilous studies found that LncRNA TTN-AS1 regulates apoptosis and drug resistance in osteosarcoma cells through the miR-134-5p/MBTD1 axis (Fu et al., 2019). Zhang et al. also investigated that LncRNA axis SNHG3/ miRNA-151a-3p/RAB22A regulates osteosarcoma invasion and migration (Zheng et al., 2019).
Advances in cancer treatment are increasingly recognizing a more promising approach to ferroptosis in the development of effective combination therapies (Friedmann Angeli et al., 2019). As first-line treatment options for patients with SARC are developed, the biology of the tumor and the tumor microenvironment should be considered in order to obtain optimal benefit from treatment strategies. The expression of NCOA4 may be a potential new factor that promotes the stratification of SARC and/or immunotherapy in patients with iron hook, which may be an important factor in predicting the recurrence of SARC patients. In addition, the modulation of SLC7A11 overexpressed in many types of cancers and is associated with patients' poor prognosis (Lin et al., 2020). In addition, FANCD2 genes were correlated with the diagnostic and prognostic factors of low-grade glioma and breast cancer (Fagerholm et al., 2013;Liu et al., 2020). CISD1 was also regarded as a potentially effective tool for prognostic of pancreatic cancer (Yang et al., 2022). It means that ferroptosis-related genes show great potential during many cancer therapies (Xu et al., 2019).
The major limitation in this study is the lack of available analysis of data heterogeneity and platform differences based on a large number of normal and tumor samples. Our study shows that the disease-related level and extreme harshness of TNM staging are not independent prediction-related factors for OS in SARC patients. This difference may be due to heterogeneity of the data or different classification and classification models (Warren and Harrison, 2018). In addition, the group cancer type information in TCGA database is mainly limited to sarcoma groups, so it is difficult to extrapolate these results to different types and locations of sarcomas. The development of this field also requires further efforts to verify the results of bioinformatics predictions, including the detection of Western blotting proteins or immunohistochemical staining, to promote the analysis of iron dependence in vivo and in vitro and immunotherapy functions.

CONCLUSION
In general, ferroptosis induction and immunotherapy are now the main advances in the treatment of SARC. With a deeper understanding of bleomycin therapy biology and resistance mechanisms, ferroptosis-based combination therapy has received more and more attention. We first found that the high expression of SLC7A11, FANCD2, CISD1, and ATP3MC3 genes related to ferroptosis, which is related to the signal transduction of the infiltrating immune cell receptor in SARC. Therefore, it is expected that SLC7A11, FANCD2, CISD1, and ATP3MC3 will be used as new markers to identify patients who may undergo adequate ferroptosis induction therapy or combined immunotherapy.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
Conception and design: ZG and SL. Acquisition, analysis, and interpretation of the data: ZG, SL, ZW, SFL, and ZG. Drafting and writing: ZG. Final approval of the article: ZG, SL, ZW, SFL, ZG, and KT.

ACKNOWLEDGMENTS
We would like to thank all participants and our hospital.