All these expression profiles and corresponding clinical data were obtained from TCGA¹. Then, complete clinical data of 548 UCEC samples and 23 normal samples including survival time were filtered for further analysis. We integrated the transcriptome and complete clinical data to screen out 511 overall survival–related UCEC samples. Half of them (n = 256) were randomly split into training cohort. The entire patients (n = 511) were defined as testing cohort to verify the signature. The baseline information is exhibited in Supplementary Table 1. Then, 60 ferroptosis-related genes were retrieved from the gene list provided by previous literature (Liang et al., 2020). A total of 15 UCEC specimens and 15 adjacent tissues were obtained from the Affiliated Tumor Hospital of Nantong University. We obtained all the written informed consent from patients.

Firstly, we performed univariate Cox regression analysis to screen targeted ferroptosis-related genes with prognostic values. To reduce the risk of over-fitting, Lasso regression analysis, and univariate and multivariate Cox regression analysis were used to construct the prognosis model (Tibshirani, 1997). Lasso algorithm was used to select variables, and “glmnet” R package was used to shrink (Simon et al., 2011). The risk score of the FRPS was calculated according to the normalized expression level of each gene and its corresponding regression coefficient. The formula was FRPS risk score = Σ(the expression amount of each gene multiplied by the corresponding coefficient). According to the median risk score of the FRPS, we divided the patients into two groups. Then, principal component analysis (PCA) was performed by using “scatterplot3d” R package on the base of expression. The “survminer” R package was used for survival analysis of each gene. The R package “survival ROC” was used to evaluate the predictive ability of the signature.

Total RNA from 15 UCEC samples and 15 adjacent tissues was extracted using TRIzol reagent (Invitrogen). The residual genomic DNA from total RNA was removed by 4 × gDNA wiper Mix (Vazyme R323-01). The complementary RNA was synthesized using PrimeScript RT reagent kit. The SYBR Premix Ex Taq Kit (TaKaRa DRR041) was used to perform real-time quantification. The relative expression levels of target genes were normalized by GAPDH and estimated using the 2^−ΔΔCt method. The primers used in this research are listed in Supplementary Table 2.

Nomogram involving prognostic clinicopathological factors (age, stage, histological type, grade, and FRPS) was carried out for prognosis prediction. In validation, we used the calibration plots for calibration of the nomogram. The “rms,” “foreign,” and “survival” package in R were used to establish and validate a nomogram (Park, 2018).

CIBERSORT algorithm was used to obtain the fraction of 22 immune cell types based on RNA-Seq data (Newman et al., 2015). The correlation between FRPS and these immune cells was analyzed by Spearman.

We used an online tool Immune Cell Abundance Identifier (ImmuCellAI) to estimate the abundance of 24 immune cells in UCEC (Miao et al., 2020). The datasets including RNA-Seq and microarray data were used to predict the patient’s response to an existing immune checkpoint blockade therapy.

The results of mRNAsi in TCGA-UCEC were obtained from a previous study (Malta et al., 2018). The mutation data of UCEC patients were downloaded from TCGA. The mutation annotation format (MAF) and MAF tool helped us to obtain somatic variation data (Mayakonda et al., 2018). The tumor mutation burden (TMB) score was obtained as follows: TMB = (total mutant bases/total covered bases) × 10⁶ (Robinson et al., 2017). The functional enrichment analysis was conducted by single-sample gene set enrichment analysis (ssGSEA) using the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways.

Microsatellite instability sensor algorithm is a program that can report the percentage of unstable microsatellites (Niu et al., 2014). We used this algorithm to obtain the MSI status for all cases based on somatic mutation data downloaded from TCGA. The relationship between FRPS and MSI was analyzed using Spearman’s rank correlation coefficient.

Gene set enrichment analysis is a method used to determine whether a set of marker genes can predict a statistically significant difference between two different cohorts. Here, we analyzed the significant difference in survival between the two cohorts in the entire TCGA cohort divided by the risk score. Normalized p-values less than 0.05 and false discovery rate (FDR) less than 0.05 are considered significantly enriched (Canzler and Hackermuller, 2020).

To evaluate the response to chemotherapy drugs, we used public pharmacogenomics database Genomics of Drug Sensitivity in Cancer (GDSC)². The half-maximal inhibitory concentration (IC₅₀) is estimated by R package “pRRophetic” (Yang et al., 2013).

To identify the molecular subtypes in endometrial cancer, the TCGA UCEC cohort was divided into different groups by R package “Consensus Cluster Plus, 1000 iterations and resampling rate of 80%” (Yu et al., 2015). We performed the log-rank test and Kaplan–Meier curve to assess the overall survival (OS) difference between different groups. Chi-square test was a good assistant helping us to compare the distribution of age, grade, stage, and histologic type between different clusters.

All statistics and figures were analyzed using R 3.6.2. Wilcoxon’s test allowed us to evaluate the differential expression of genes related to ferroptosis between UCEC patients and controls. We used the χ² test to assess the relationship between FRPS and clinicopathological factors. P-value < 0.05 was considered to be statistically significant.

First, the expression profiles and survival data of UCEC patients in training cohort were filtered based on 60 ferroptosis-related genes through the univariate Cox regression analysis and 17 genes were correlated with overall survival (Table 1). Then, LASSO analysis (Figures 1A,B), and univariate and multivariate Cox regression analysis narrowed the screened scope to six genes (HMOX1, KEAP1, HSBP1, SAT1, CISD1, and GPX4) (Supplementary Material). The risk score of FRPS for OS = (0.002907 × HMOX1) + (0.013486 × KEAP1) + (−0.089640 × HSBP1) + (−0.001665 × SAT1) + (0.148239 × CISD1) + (−0.003060 × GPX4). Meanwhile, patients in the training cohort were divided into high-risk and low-risk groups according to the median risk score of the FRPS (Figure 1C). PCA indicated the patients in different risk groups were distributed in two directions (Figure 1D). The distribution of FRPS and survival status of patients in OS signature are shown in Figure 1E. The OS of patients in the low-risk group was significantly longer than those in the high-risk group according to the Kaplan–Meier curve (Figure 1F). Besides, the 1-year, 3-year, and 5-year ROC curves based on training cohort are plotted as Figure 1G, suggesting satisfactory prognostic value of the signature.

Then, the mRNA expression of these genes was validated by qPCR using the samples from the Nantong Third People’s Hospital Affiliated to Nantong University (Figure 2). The mRNA expression of KEAP1, HSBP1, SAT1, CISD1, and GPX4 were significantly different between tumor and the adjacent tissues. KEAP1 and HOMOX1 were low expressed in tumor than the para-carcinoma tissues, and the high expression of these two genes increased the risk of FRPS, which suggested the poor prognosis of patients. These results suggested the potential feasibility of this signature for clinical usage. Expression and Kaplan–Meier survival analysis of each gene in the signature were also performed and four genes were output significantly (Supplementary Figure 1).

To test the robustness of the aforementioned signature, the entire TCGA patients were divided into high-risk and low-risk groups by the same risk score (Figure 3A). PCA confirmed the similar results obtained from the training cohort; the patients in the two subgroups were distributed in discrete directions (Figure 3B). Similarly, patients in the high-risk group showed worse prognosis (Figure 3C). The expression of six genes in the signature are exhibited in Figure 3D. Kaplan–Meier survival analysis claimed the reduced survival time of patients in the high-risk group compared with those in the low-risk group (Figure 3E). Furthermore, subgroup analyses in age, grade, stage, and histological type in the signature were performed, demonstrating that patients with high-risk scores shared shorter OS in all the subgroups (p < 0.05) (Supplementary Figure 2). Besides, the 1-year, 3-year, and 5-year AUC of signature was 0.705, 0.676, and 0.713, respectively (Figure 3F). According to the Cox regression analysis, the histological type, tumor stage, and this signature were independent prognostic factors (Table 2). We further compared the prediction value of this signature with other existing signatures. As shown in Figure 4, the AUC for this signature was 0.705, which was higher than the existing ferroptosis-related signatures. This result revealed that our established signature was superior to other signatures in predicting patient’s survival information (Yu et al., 2015; Liu et al., 2021; Yang et al., 2021; Zhu et al., 2021).

To better predict overall survival time, we integrated clinicopathological factors related to prognosis (age, stage, histological type, and FRPS) to establish a nomogram prediction model (Figure 5A). We compared the relationship between FRPS and clinicopathological factors (Supplementary Table 3). Quantifying the aforementioned variables as numerical points, 1-year, 3-year, and 5-year survival rates of UCEC patients can be calculated based on the total points of all risk factors. A calibration chart was also constructed to show the consensus of the predicted and observed results (Figure 5B). Meanwhile, ROC curve demonstrated the better predictive ability of FRPS in 1-year, 3-year, and 5-year OS than other clinical factors (Figure 5C). Combining clinical factors and FRPS accessed optimal predicting effect in UCEC patients based on 1-year, 3-year, and 5-year OS (Figure 5D).

Using ESTIMATE algorithm, patients in the high-risk group were found to have lower immune scores, stromal scores, and ESTIMATE scores (Figures 6A–C). On the contrary, patients accessed higher tumor purity scores in the high-risk group (Figure 6D). The aforementioned findings suggested that the tumor immune microenvironment was closely associated with the FRPS in UCEC patients. To find the major immune cells between the high-risk groups and low-risk groups, CIBERSORT algorithm was employed. The results showed that macrophages M1, macrophages M2, T cell follicular helper, and B cells naive were positively correlated with FRPS while NK cells activated, T cells regulatory (Tregs), neutrophils, and dendritic cells resting were negatively correlated with FRPS (Figure 6E). The distribution of immune cells and scores for each patient is exhibited in Figure 6F. In the present research, we also focused on the tumor infiltrating cells between subgroups. We found that B cells naive, macrophages M1, macrophages M2, T cells CD4 memory activated, T cells follicular helper, T cells gamma delta, NK cells activated, T cells regulatory (Tregs), and neutrophils infiltrated differently in different groups (Figure 6G).

Immune checkpoint proteins, playing important roles in the immune response, aroused our great interest to explore the relationship between FRPS and immune checkpoint modulators. The results demonstrated that CD40, PD-L1, and PD-L2 showed a positive correlation with FRPS while CD270, CD27, and CTLA4 were negatively related to FRPS (Figure 7A). The distribution of immune checkpoint proteins and risk scores for each patient is exhibited in Figure 7B. Immune checkpoint proteins between high-risk groups and low-risk groups were evaluated, and results demonstrated that CD27, CTLA4, PD-L2, B7-H4, CD40, PD-L1, and CD270 expressed differently between high- and low-risk groups (Figure 7C). Potential response to immunotherapy in each patient was further assessed by the online tool ImmuCellAI and patients in the low-risk group showed better reactivity to immunotherapy than those in the high-risk group (p = 0.016; Figure 7D).

Tumor mutation burden (TMB) is an important cause of tumor occurrence and development, can be used to predict the efficacy of immune checkpoint blockade, and has been shown to be a biomarker for patients who benefit from immunotherapy. In this study, we declared that FRPS was negatively correlated with TMB (Figure 8A). Lower TMB was observed in the high-risk group (Figure 8B). In addition, the mutant genes that showed the most significant difference in their mutation frequency between the two groups are shown in Figure 8C. TP53 and PPP2R1A were found to have higher mutation frequency in the high-risk group, and the rest of the genes showed higher mutation frequency in the low-risk group. Therefore, somatic mutation data were used to assess the TMB of patients. The order of somatic mutation frequency in the high-risk group was TP53 > PTEN > PIK3CA > TTN > ARID1A > PIK3R1 > KMT2D > CHD4 > MUC16 > PPP2R1A (Figure 8D); in the low-risk group, PTEN > ARID1A > PIK3CA > TTN > CTCF > PIK3R1 > CTNNB1 > KMT2D > ZFHX3 > KRAS (Figure 8E).

Several researches had illustrated that MSI can affect the effect of immunotherapy in several cancers. In this research, we also investigated the MSI status between groups. The results revealed that MSI status was negatively correlated with FRPS (Figure 8F). Besides, according to ImmuCellAI, higher immunotherapy response rate was observed in the low-risk group compared with patients in the high-risk group (Figure 8G), which implied that patients in the low-risk group might benefit from immunotherapy.

In recent years, the role of m6A methylation in cancer has attracted widespread attention. More and more evidence showed that the genetic changes and expression disorders of m6A RNA are associated with the tumor occurrence, progression, and treatment resistance. The expression levels of HNRNPC, YTHDC1, ZC3H13, YTHDF2, FTO, YTHDF1, YTHDF3, METTL14, RBM15, WTAP, KIAA1429, FMR1, and HNRNPA2B1 were dramatically higher in UCEC high-risk group (p < 0.05) (Supplementary Figure 3A). In addition, the expression levels of mRNAsi (p = 3.52e-10) and EREG-mRNAsi (p = 0.032) in high-risk group were also higher (Supplementary Figure 3B). To further explore the correlation between FRPS and immune status, we used ssGSEA to quantify the enrichment scores of various immune cell subgroups, related functions, or pathways (Supplementary Figures 3C,D). Interestingly, the antigen presentation process (aDC and iDC) was significantly different, and the enriched cytokine–cytokine receptor interaction in the KEGG analysis scored higher in the high-risk group (p < 0.05, Supplementary Figure 3D).

In addition, GSEA analyzed the transcription information of patients in the high-risk and low-risk subgroups. Based on normalized p-value < 0.05, FDR < 0.05, and NES, we filtered the most significant enrichment biological approach. Representative KEGG pathways were related to some essential signaling pathways including cell cycle, DNA replication, mismatch repair, alpha linolenic acid metabolism, and ribosome and tyrosine metabolism (Supplementary Figure 4). The aforementioned results suggested the potential mechanism of the occurrence and development of UCEC.

Using the pRRophetic algorithm, IC₅₀ of 35 common chemotherapeutic agents was predicted in high- or low-risk group (PD.0332991, Nutlin.3a, X17.AAG, Bryostatin.1, PD.0325901, SB.216763, Bicalutamide, AZD6244, PF.02341066, LFM.A13, Temsirolimus, NVP.BEZ235, FTI.277, RDEA119, BMS.536924, MG.132, PF.562271, Roscovitine, AZ628, Vinblastine, EHT.1864, Tipifarnib, BMS.754807, Lapatinib, KIN001.135). All 25 drugs had higher IC₅₀ in high-risk patients, indicating that the low-risk patients were more sensitive to these 25 drugs (Wilcoxon test, all p < 0.05; Supplementary Figure 5).

Consensus clustering was analyzed based on the expression levels of six targeted genes. We chose K = 2 as the most optimal clustering because the clustering was suboptimal when divided into more than two clusters (Figures 9A–C). UCEC patients in clusters 1 and 2 showed significant differences in age, stage, and histological type, but did not show any significant differences in grade (Figures 9D,E). Moreover, the OS was significantly shorter in the UCEC patients of cluster 1 based on Kaplan–Meier curve compared with those of cluster 2 (Figure 9F). We also compared the significantly enriched KEGG pathways between two clusters; four pathways were identified, including fatty acid metabolism, graft-versus-host disease, protein export, and ribosome. These mechanisms may involve in the pathogenesis of UCEC (Figure 9G).