Immunological and prognostic significance of novel ferroptosis-related genes in soft tissue sarcoma

Background Ferroptosis has exhibited great potential in the treatment of cancer and has gained widespread attention in soft tissue sarcoma (STS). The aim was to explore the immunological and prognostic significance of novel ferroptosis-related genes in STS. Methods We identified ferroptosis-related differentially expressed genes (DEGs) in STS to construct the networks of enrichment analysis and protein-protein interaction. Subsequently, hub genes with prognostic significance were localized and a series of prognostic and immune analyses were performed. Results 40 ferroptosis-related DEGs were identified, of which HELLS, STMN1 EPAS1, CXCL2, NQO1, and IL6 were classified as hub genes and were associated with the prognosis in STS patients. In the results of the immune analysis, PDCD1, CTLA4, TIGIT, IDO1 and CD27 exhibited consistent intense correlations as immune checkpoint genes, as well as macrophage, neutrophil, cytotoxic cell, dendritic cell, interdigitating dendritic cell and plasmacytoid dendritic cell as immune cells. EPAS1 and HELLS might be independent prognostic factors for STS patients, and separate prognostic models were constructed by using them. Conclusions We recognized novel ferroptosis-related genes with prognostic value in STS. Furthermore, we searched out potential immune checkpoints and critical immune cells.


Introduction
Soft tissue sarcoma (STS) is a group of malignant tumors originating from mesenchymal tissue and containing multiple histological subtypes [1]. The prognosis of partial STS is poor with no effective treatment and the precise prediction of the prognosis for STS patients is a challenging topic [2]. The previous view was that immunotherapy was unpromising in STS, but this has been reversed in recent years [3].
Ferroptosis is an emerging phenotype of regulated cell death (RCD) which relies on reactive oxygen species deposition mediated by iron catalysis and lipid peroxidation [4]. Ferroptosis performs an essential role in the initiation, progression and prognosis of multiple diseases [5]. Meanwhile, ferroptosis has exhibited great potential in the treatment of cancer and has gained widespread attention in STS as well [6]. Recent studies have revealed that ferroptosis and tumor immunity can be mutually regulated [7,8].
In the present study, differentially expressed genes (DEGs) were identified through the Gene Expression Omnibus (GEO) database, the FerrDb database, the Immunology Database and Analysis Portal (ImmPort) database, and the networks of enrichment analysis and protein-protein interaction (PPI) were constructed. Prognostic and immune analyses were performed through the Cancer Genome Atlas (TCGA) database. The aim was to explore the immunological and prognostic significance of novel ferroptosis-related genes in STS.

Functional enrichment analysis and PPI networks construction
Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and Gene Set Enrichment Analysis (GSEA) were implemented through the clusterProfiler package of R [14]. FDR < 0.05 for the enriched item was considered statistically significant. After predicting the interactions between DEGs in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/ ) [15] by setting the combined score > 0.4, the PPI networks were built using Cytoscape [16] and cytoHubba [17] respectively.

Hub genes identification and prognostic models construction
Through the TCGA database, high expression and low expression groups were divided by the median of DEGs expression and survival analysis was performed with the survival package of R [18]. DEGs with potential prognostic significance were identified as hub genes by log-rank analysis and visualization was achieved through the survminer package. The Wilcoxon rank sum test was chosen for correlation analysis of hub genes expression with clinical variables, and the ggplot2 package was used for visualization. All clinical variables of STS were integrated into univariate Cox regression, parameters were included in overall survival (OS) and progression free survival (PFS), and variables that were significant for univariate analysis were integrated into multivariate Cox regression. After evaluating significant variables in the multivariate analysis by the timeROC package, they were incorporated into a nomogram to construct the model [19]. The population of the model was 263 patients with well-defined STS, from the TCGA database and screened with corresponding clinical information, and the model was validated by a calibration curve, with visualization implemented through the rms package. The results were considered statistically significant at p < 0.05.

Immune analysis
Through the TCGA database, high expression and low expression groups were classified according to the upper and lower quartiles of DEGs expression and the GSVA package of R accompanied by Spearman correlation test was applied for immune analysis [20]. 7 popular immune checkpoint genes (ICGs) [21,22] and 24 immune cells composing the main tumor immune microenvironment [23] were included by applying the CIBERSORT deconvolution algorithm, and the ggplot2 package was used to construct co-expression plots. The results were considered statistically significant at p < 0.05.

Statistical analysis
Statistical analysis relied on R software (version 3.6.3) and Cytoscape software (version 3.8.2).

Ferroptosis-related DEGs identification in STS
A total of 927 DEGs for STS were identified in 183 samples from the GSE21122, GSE6481 and GSE2719 datasets, including 345 for up-regulation and 582 for down-regulation. The volcano plot covered all genes in differential analysis ( Fig 1A) and the heatmap displayed the top 20 DEGs for each of up-regulation and down-regulation ( Fig 1B). Among them, a total of 40 genes were associated with ferroptosis (Table 1), including 6 for up-regulation and 34 for down-regulation ( Fig 1C).

PLOS ONE
Novel ferroptosis-related genes

PPI networks construction
Interactions of ferroptosis-related DEGs in STS were predicted by STRING and a PPI network covering 38 nodes and 97 edges was structured by Cytoscape ( Fig 3A). Subsequently we used cytoHubba to further identify the top 25 genes and build a 25-node, 69-edge PPI network ( Fig 3B).

Hub genes identification
Survival analysis revealed the potential prognostic value of HELLS, STMN1 in up-regulation DEGs and EPAS1, CXCL2, NQO1, IL6 in down-regulation DEGs, with high expression of HELLS, STMN1 and low expression of EPAS1, CXCL2, NQO1, IL6 suggesting a short OS in STS patients ( Fig 4A). Accordingly, these 6 genes were identified as hub genes for further  study and the association between hub genes and STS clinical variables was analyzed (Fig 4B). Compared to male patients with STS, female patients exhibited high expression of HELLS, STMN1, NQO1 and low expression of EPAS1. Compared to other histological types, leiomyosarcoma showed high expression of HELLS, EPAS1, NQO1 and low expression of IL6. Besides, HELLS and STMN1 was highly expressed in STS metastatic patients compared to non-metastatic patients.

Hub genes GSEA analysis
The 263 STS samples from TCGA database were divided into low expression and high expression groups based on the median of hub gene expression respectively for GSEA analysis. GSEA manifested significant differences in enrichment of MSigDB Collection (FDR < 0.05) and significant-enriched gene sets were ranked based on normalized enrichment score (NES) values. The top-two most significant-enriched gene sets for HELLS were G alpha signaling events and olfactory transduction (Fig 5A). The top-two most significant-enriched gene sets for STMN1 were signaling by Rho GTPases and processing of capped intron-containing pre-mRNA ( Fig 5B). The top-two most significant-enriched gene sets for EPAS1 were M-phase and metabolism of amino acids and derivatives ( Fig 5C). The top-two most significantenriched gene sets for CXCL2 were neuronal system and neuroactive ligand receptor interaction ( Fig 5D). The top-two most significant-enriched gene sets for NQO1 were signaling by interleukins and Leishmania infection (Fig 5E). The top-two most significant-enriched gene sets for IL6 were signaling by interleukins and GPCR-ligand binding (Fig 5F).

EPAS1 and HELLS might be independent prognostic factors for STS patients
Clinical variables for STS and the expression of 6 hub genes were included in univariate Cox regression analysis and those factors of significance were further subsumed into multivariate Cox regression analysis. The results revealed that when the prognostic indicator was OS, high grade residual tumor, metastasis, positive margin status, high expression of HELLS and STMN1, low expression of EPAS1, CXCL2, NQO1, IL6 were associated with poor prognosis. Furthermore, residual tumor, metastasis status, margin status, EPAS1 expression might be independent prognostic factors for OS in STS patients (Table 2). When the prognostic indicator was PFS, high grade residual tumor, metastasis, positive margin status, HELLS high expression were associated with poor prognosis. Residual tumor, metastasis status, margin status, HELLS expression might be independent prognostic factors for PFS in STS patients ( Table 3).

Validation of EPAS1 and HELLS prognostic value
Predictive efficacy of EPAS1 and HELLS for prognosis was internally verified using timedependent receiver operating characteristic (ROC) curves in TCGA database (Fig 6A and 6B). Subsequently, predictive efficacy of EPAS1 and HELLS for prognosis was externally validated using time-dependent ROC curves in GEO database, which exhibited similar prognostic value (Fig 6C and 6D).

Construction and evaluation of prognostic models for STS patients
The statistically significant results of the multivariate Cox regression analysis were used to construct the separate nomogram for prediction models of OS ( Fig 7A) and PFS (Fig 7B) in STS patients. For both patients with primary STS and metastatic STS, the indicators for each nomogram were derived from the primary tumor foci. The C-indexes for OS and PFS model were 0.756 (0.719-0.794) and 0.782 (0.756-0.808) respectively. Calibration curves for the models of OS ( Fig 7C) and PFS (Fig 7D) confirmed the consistency of the predicted prognosis with the actual outcome.

Association of hub genes expression and ICGs
CXCL2 and IL6 were shown to be immunologically relevant in 6 hub genes (Fig 8A). We correlated hub genes with ICGs and presented the results in co-expression heatmaps. CXCL2 and IL6 showed consistent results, with both CXCL2 and IL6 positively linked to the expression of PDCD1, CTLA4, TIGIT, IDO1 and CD27 (Fig 8B and 8C). Consistency of results and significant association with ICGs were not demonstrated in HELLS, STMN1, EPAS1 and NQO1 (Fig 8D-8G).

Association of hub genes expression and immune cells infiltration
6 hub genes were subsequently correlated with 24 immune cells in the tumor microenvironment ( Fig 9A). In addition to CXCL2 and IL6, we observed that HELLS was also strongly associated with immune cells and exhibited the consistent result with CXCL2 and IL6. CXCL2 and IL6 with down-regulated in STS were significantly positively related to macrophage, neutrophil, cytotoxic cell, dendritic cell (DC), interdigitating dendritic cell (iDC), plasmacytoid dendritic cell (pDC) (all r > 0.3) (Fig 9B and 9C), and HELLS with up-regulated in STS was comparatively negatively correlated with macrophage, neutrophil, cytotoxic cell, DC, iDC, pDC (all r < -0.3) (Fig 9D) Consistency of results and significant association with immune cells were not demonstrated in STMN1, EPAS1 and NQO1.

Discussion
STS is a set of heterogeneous malignancies involving over 100 different histological types, with widely varying treatment outcomes [24]. In general, current therapies are only effective in a small proportion of STS, with limited efficacy in most STS and even recurrence in more than 50% of patients [25]. Although it was considered that STS was extremely insensitive to immune responses in the past, which precluded the application of immunotherapy to STS, recent studies have demonstrated a large degree of immune heterogeneity within the subclass of STS and some positive responses to immunotherapy have also been reported in successive clinical trials [26,27]. Partial STS subtypes, including dedifferentiated liposarcoma, leiomyosarcoma, embryonal rhabdomyosarcoma and undifferentiated pleomorphic sarcoma, have been identified as featuring high levels of immune cells infiltration and ICGs expression, and exhibit a potentially active reaction to immune checkpoint inhibitors (ICIs) therapy [3]. Consequently it is essential to locate critical ICGs and immune infiltration factors adapted to STS. Currently, the availability of immunotherapy alone is severely limited in patients with most tumor types. Since extensive crossover between immunotherapy and non-apoptotic RCD mechanisms has been detected, non-apoptotic cancer cell death accompanied by immunomodulation is considered an exceedingly promising strategy for cancer treatment [28]. Ferroptosis, a neoteric form of RCD with unique biological and morphological features, has been shown to interact with the tumor immune response and can influence immunotherapeutic efficacy on the one hand [8], and in turn is regulated by immune cells on the other [7]. In the present study, based on ferroptosis-related genes in STS, we identified potential ICGs including PDCD1, CTLA4, TIGIT, IDO1 and CD27, which might serve as important targets for immunotherapy. In addition, we explored a group of closely related immune cells including macrophage, neutrophil, cytotoxic cell, DC, iDC and pDC, which might act as pivotal regulators in the immune microenvironment of STS. Interestingly, we observed high concordance of immune analysis results for HELLS with CXCL2 and IL6, revealing for the first time a possible immunological effect of HELLS in tumor.
Among dedifferentiated liposarcoma, undifferentiated pleomorphic sarcoma and leiomyosarcoma, it has been confirmed that tumors with high immunogenic gene profiles are accompanied by high levels of PDCD1 expression [29]. PD-1, as the most researched immune checkpoint, is encoded by PDCD1 and also occupies an important position in STS study. More than half of the samples in a STS cohort had positive expression of PD-1 on immune cells [30], and PD-1 expression is also generally considered to be associated with the prognosis of STS patients [31,32]. Moreover, CTLA4, IDO1 and other ICGs have demonstrated varying degrees of value for STS management [3]. In terms of immune cells, macrophage has been established as a significant player in several sarcoma types [33], with the modification of the macrophage phenotype from tumor-promoting to tumor-suppressing regarded as a promising option for STS treatment [34]. And a range of immunotherapies targeting DC, iDC and pDC may be well tolerated in patients with refractory STS due to their excellent immunological response and safety profile, as well as offering the opportunity to prevent recurrence of sarcoma [35]. On balance, for most STS subtypes, immunotherapy may be required novel regimens and combinations [34].
In addition, we substantiated that EPAS1 and HELLS might act as independent prognostic predictors of STS, leading to the construction of two efficient prognostic models. For both patients with primary STS and metastatic STS, the indicators for each nomogram were derived from the primary tumor foci. However, the obtained model still needs to be further verified in an independent cohort. The expression of 6 hub genes was discovered to be associated with survival during the model construction, with EPAS1, STMN1, CXCL2, NQO1 being identified for the first time in STS. EPAS1 is a diver of ferroptosis [36], compared to normal tissue, which is expressed at lower levels in most human STS [37]. Zhu et al. found that up-regulation of EPAS1 significantly enhanced the growth inhibition of gastric adenocarcinoma and that targeting EPAS1 might be an alternative therapeutic approach for cancer [38]. Relatively, HELLS, NQO1 are suppressors of ferroptosis [39,40]. Law et al. suggested that HELLS mediated epigenetic silencing of various cancer suppressor genes and evidenced in hepatocellular carcinoma that its overexpression potentiated tumor cell migration and proliferation [41]. Huang et al. identified high expression and prognostic impact of HELLS in STS samples [42], which also underpinned our findings. In the TCGA database of STS samples, GSEA indicated that NQO1 was closely connected to interleukin-related signaling pathways. NQO1 has been confirmed to interact with interleukins in a variety of cancers, thereby affecting the inflammatory response and participating in the immune regulation associated with the tumor microenvironment [43,44]. As for STMN1, CXCL2 and IL6, they are currently treated as biomarkers of ferroptosis and their expression is monitored for down-regulation once ferroptosis occurs [45,46]. STMN1 is commonly recognized as an oncogene, and its up-regulation is tightly linked to the malignant behaviour and poor prognosis of various tumors [47]. In leiomyosarcoma, STMN1 has also been characterised by high expression and can be a sensitive biomarker with strong diagnostic efficacy [48,49]. Our study revealed the potential immunological relevance and clinical value of these novel ferroptosis-related genes, which might contribute to the precise treatment and prognostic prediction of patients with STS.

Conclusions
In conclusion, we identified novel ferroptosis-related genes with prognostic value in STS. Furthermore, we searched out potential immune checkpoints and critical immune cells.