The somatic mutation, transcriptome, CNV, and sample phenotype data of EC were downloaded from the TCGA xena database. Among them, there were 530 samples with somatic mutation data; there were two sets of CNV data, one of which consisted of the CNV results compiled by xena, showing CNV at the gene level, and the other was the sample DNA copy file downloaded from http://portal.gdc.cancer.gov/, which is convenient for grouping and performing later gistic2 analysis. One set of RNA-Seq data, labeled TCGA-EC, consists of a total of 583 samples. After removing four samples without survival information, the expression data of 579 samples were finally obtained, including 544 tumor tissues (clinical information shown in Table 1) and 35 normal tissues. The downloaded expression profile data format was log2(FPKM + 1). The expression value was restored to FPKM by the formula FPKM = 2^{original expression value} −1 and then converted to TPM by the following formula: and finally, log2(TPM + 1) conversion was performed. This data was used for subsequent analysis.

The gene sets of 23 infiltrating cells were obtained from a reference (Zhang B. et al., 2020). The gene sets of angiogenesis, CD8 T effector, EMT1, EMT2, EMT3, and panfibroblast TGFb were downloaded from another reference (Mariathasan et al., 2018). From an additional reference (Liang et al., 2020), 60 ferroptosis genes were obtained (see Supplementary Table 1). The c2.cp.kegg.v7.1.symbols with a total of 186 gene sets were downloaded from MsigDB.

Using TCGA-EC RNA-Seq data, the expression values of 60 ferroptosis genes were extracted, and the R package ggpubr was used to draw box plots to show the differences in the expression of these genes between tumor tissues and normal tissues. The CNV results of ferroptosis genes were extracted, the frequency of amplification and deletion was counted, and a dot plot was drawn. The R package maftools was used to import the maf files of somatic mutations in 530 EC samples and draw the somatic mutation spectrum. The R package RCircos was used to draw a circos map of 60 ferroptosis genes, showing the positions of the genes on the reference genome. The R package pca3d was used to perform principal component analysis (PCA) on the expression matrix of 60 ferroptosis genes and draw a three-dimensional PCA map.

Using the expression matrix of ferroptosis genes in 544 EC samples as input files, the R package ConsensusClusterPlus was used to perform unsupervised clustering with the following parameters: m (maximum number of categories) = 6, reps (repeated sampling) = 1,000, pItem (proportion of items selected each time) = 0.8, pFeature (the proportion of features selected each time) = 1, the clusterAlg (clustering algorithm) = “pam”, distance (calculation distance) = “spearman.” The output results were synthesized, the K values were filtered, and the classification of each sample was obtained.

For Fe.cluster, grouped two by two, the expression matrix of all genes of the two samples was extracted, the c2.cp.kegg.v7.1.symbols gene set was combined as the input file, gene set variation analysis [GSVA; a GSE method that estimates the variation of pathway activity over a sample population in an unsupervised manner (Hänzelmann et al., 2013)] was performed, and the enrichment score of each sample for each gene set was obtained. Then, the R package limma was used to analyze the differences in the gene sets, and the threshold was a BH-corrected p-value < 0.05. The top 20 differential gene sets were extracted, the R package ComplexHeatmap was used to draw the differential gene set heat map, and group labels were added.

The gene expression matrix of 544 tumor samples and the gene set of 23 infiltrating cells were used as the input files of the R package GSVA for single-sample GSEA (ssGSEA). Using the enrichment score as the content of each cell, a box plot was drawn using the R package ggpubr, and the Kruskal–Wallis rank sum test was performed to show the differences in the content of 23 infiltrating cells between the Fe.cluster groups. Differential cells were selected, the R package survival was used to perform single-factor Cox risk regression analysis to obtain the hazard ratio (HR) and p-values of differential cells, and then the R package forestplot was used to draw forest plots to visually display the prognostic effects of differential cells.

The display of angiogenesis, CD8 T effectors, EMT1, EMT2, EMT3, and panfibroblast TGFb enrichment scores between different Fe.cluster groups was also performed for ssGSEA following the method described above; the enrichment score of each sample was obtained for these six biological functions, the R package ggpubr was used to draw box plots, and the Kruskal–Wallis rank sum test was performed.

The Fe.clusters were grouped in pairs, the expression matrix was extracted, and the R package limma was used for differential gene analysis, with a differential gene screening threshold of abs[log2(fold change)] > log2(1.25) and BH correction p-value < 0.05. Multiple sets of differential genes were obtained, and the intersection was taken as the final differential gene. The R package clusterProfiler was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, and the screening threshold was BH correction p-value < 0.05 and Q-value < 0.05. The top 20 enrichment items of biological processes (BP), cell components (CC), molecular function (MF), and KEGG were selected, and a bubble chart was drawn to display the results.

The differential gene expression matrix of 544 tumor samples was extracted as the input file, and the R package ConsensusClusterPlus was used to perform unsupervised clustering. The parameters were set as follows: m = 6, reps = 1,000, pItem = 0.8, pFeature = 1, clusterAlg = “pam,” and distance = “spearman.” The output results were synthesized, the K values were filtered, and the classification of each sample was obtained. Sample clinical data were integrated, survival analysis was performed on Fe.gene.cluster, and Kaplan–Meier curves were drawn. The expression matrix of 60 ferroptosis genes was extracted, and the R package ggpubr was used to draw a box plot to show the expression in different Fe.gene.clusters and perform the Kruskal–Wallis rank sum test.

For differentially expressed genes, the R package survival was used to perform single-factor Cox risk regression analysis, and genes were screened according to a p-value < 0.05. For the retained genes, the R package randomForestSRC was used to construct a random forest model, and then important feature variables were screened as ferroptosis gene signatures.

The expression matrix of ferroptosis gene signatures was selected, PCA was performed, the two principal components were selected as PC1 and PC2, and the ferroptosis score (Fe.score) was calculated according to the following formula:

Sample clinical data were integrated and divided into high Fe.score and low Fe.score according to the median Fe.score. The R package survival was used to analyze the survival of the Fe.score, and the R package survivalROC was used to draw a 5-year ROC curve to evaluate the Fe.score survival model.

The ferroptosis gene and the top 10 differential gene sets of the three groups of GSVA were integrated to perform ssGSEA, and then the enrichment score was combined with the Fe.score to calculate the Pearson correlation coefficient matrix and calculate the p-value. The R package corrplot was used to draw the correlation diagram, and a p-value < 0.01 was considered to be extremely significant. Dots of the corresponding colors were drawn in the figure according to the correlation coefficient.

The R package pRRophetic was used to perform GDSC drug IC50 prediction, and a box plot was drawn to show the IC50 difference between Fe.score groups. The bladder cancer data set (IMvigor210) (Geeleher et al., 2014) was used to calculate the Fe.score, and the differences in immunotherapy results between groups were analyzed. The model in this paper was compared with another model (Deng et al., 2020), and the R package survivalROC was used to draw a 5-year ROC curve to evaluate the Fe.score survival model.

The comparison of the two groups in the box plot in this paper uses the nonparametric Wilcox rank sum test; the comparison of multiple groups uses the Kruskal–Wallis rank sum test. The R package ggalluvial was used to draw the Fe.cluster, Fe.gene.cluster, and Fe.score grouping of 544 samples and survival status mulberry charts. The R package maftools was used to draw the somatic mutation map and CNV peak map, and the Gistic2 analysis result file of the sample was required when drawing CNV. However, there was no complete Gistic2 result in the TCGA database, so DNA copy was used to perform Gistic2 analysis first. For analysis methods and parameter settings, refer to http://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/. The R package ComplexHeatmap was used to draw heat maps. The Spearman correlation coefficient between ferroptosis genes was calculated, and screening was performed according to a p-value < 0.001 and an absolute value of the correlation coefficient greater than 0.2. Using consensus clustering to cluster genes, single-factor Cox risk regression was used to determine the prognostic effects of genes, and the results were sorted into tables and imported into Cytoscape (3.7.2) to draw gene interaction network diagrams.

We first screened 60 ferroptosis genes, and the positions of 60 ferroptosis genes in the reference genome hg38 are shown in Figure 1A. Then, we analyzed the expression of the genes and found that, using the expression values of 60 ferroptosis genes for PCA, normal tissues and EC tissues could be clearly distinguished (Figure 1B). Figure 1C shows the detailed gene expression differences between EC tissues and normal tissues. Except for the differences in the expression of the ACSL4, AKR1C3, ALOX5, CBS, EMC2, GCLM, GLS2, HSPB1, KEAP1, NOX1, and RPL8 genes, which were not significantly different between normal tissues and cancer tissues, the expression of the other genes was significantly different. We further studied the mutations of ferroptosis genes and found many missense mutations in the TP53 gene in the EC samples (Figure 1D). For the CNV mutation frequency of the ferroptosis gene, we found that GPX4, PGD, and CHAC1 had a higher frequency (16.8, 8.2, and 6.6%) of fragment deletions, and TFRC, KEAP1, PHKG2, and SQLE had a higher frequency (15.0, 8.8, 8.2, and 8.0%) of fragment amplifications (Figure 1E).

We first performed cluster analysis of ferroptosis genes and divided them into four clusters: regulator cluster A, regulator cluster B, regulator cluster C, and regulator cluster D (see Supplementary Table 1 for the detailed results of the ferroptosis clusters). GPX4, SAT1, and TP53 in the ferroptosis gene regulatory network were significant prognostic protective factors. The CBS, CHAC1, and CISD1 genes were significant prognostic risk factors (see Supplementary Table 1 for details). Then, the correlation coefficient between ferroptosis genes was calculated, and statistical tests were performed. When the p-value < 0.001 and the absolute value of the correlation coefficient was greater than 0.2, the genes were considered to have an interaction relationship. As shown in Figure 2A, the larger the correlation coefficient, the thicker the connected line.

Then, we performed a cluster analysis of 544 EC samples based on the expression values of 60 ferroptosis genes (see Supplementary Table 2 for the classification results) and finally divided the samples into four clusters. In Figure 2B, 1, 2, 3, and 4 correspond to Fe.cluster A, Fe.cluster B, Fe.cluster C, and Fe.cluster D, and the number of samples in each cluster are 183, 175, 153, and 33, respectively. As shown in Figure 2C, most of the ferroptosis genes were highly expressed in Fe.cluster B. The expression pattern of the ferroptosis gene in Fe.cluster D was quite different from that in the other three groups, and the expression of part of the ferroptosis gene (PGD, G6PD, AKR1C1, AKR1C2, AKR1C3, NQO1) was particularly high in Fe.cluster D.

Then, we explored how the Fe.clusters reflect the immune status of EC. As shown in Figure 2D, 21 kinds of immune cells demonstrated significant differences in different Fe.clusters except for activated CD4 T cells and type 2 T helper cells. We also observed that activated CD8 T cells, eosinophils, CD56dim natural killer cells, and activated B cells were good prognostic factors for EC (Figure 2E).

We next analyzed the specific molecular characteristics of Fe.clusters. The GSVA functional enrichment analysis results of Fe.cluster A, Fe.cluster B, and Fe.cluster C are shown in Figures 3A–C. It can be seen that most of the pathway enrichment scores of Fe.cluster B are higher than Fe.cluster A, and Fe.cluster A is higher than Fe.cluster C. We also focused on analyzing several key BP, and we found that the enrichment scores of angiogenesis, CD8 T effector, EMT2, and panfibroblast TGFb (Pan-F-TBRS) were significantly different in Fe.clusters (Figure 3D).

Then, we explored the relationship between clusters and clinical indicators. Figure 3E survival curve shows the survival difference between Fe.cluster groups: the log-rank test p-value is 0.021, the difference is significant, and Fe.cluster A has a higher survival time than the other three groups. In addition, as shown in Figures 3F–Q, for tumor stage, the proportion of Fe.cluster A and Fe.cluster B stage I patients is higher than the other two groups; for tumor grade, the proportions of G1 and of G2 patients in Fe.cluster A and Fe.cluster B are higher than in the other two groups. The proportion of patients with diabetes, hypertension, and more pregnancies in Fe.cluster A is lower than in the other three groups; menopausal status, age, radiotherapy, and BMI were not significantly different between the four groups. In Fe.cluster C, patients with PTEN mutation, TP53 wild-type, KRAS wild-type, and APC wild-type accounted for lower proportions than in the other three groups.

Considering that the ferroptosis gene expression pattern and clinical features of Fe.cluster D are quite different from those of the other three groups and the number of samples in this group is relatively small, a difference analysis only on the samples of Fe.cluster A, Fe.cluster B, and Fe.cluster C was performed, and 438 overlapping differential genes were obtained (see Figure 4A and Supplementary Table 3). Then, using the expression matrix of differential genes to perform consensus clustering, three groups were obtained, namely, Fe.gene.cluster A, Fe.gene.cluster B, and Fe.gene.cluster C. The sample numbers were 78, 147, and 319, respectively (see Supplementary Table 2 for the results of the Fe.gene.cluster). The expression heat map of 438 differentially expressed genes is shown in Figure 4A, and we can also see the relationship between Fe.gene.clusters and clinical indicators (stage, age, and status) from Figure 4B. Furthermore, we found that there was a very significant difference in survival between the Fe.gene.cluster groups, and Fe.gene.cluster A and Fe.gene.cluster B had obvious survival advantages (Figure 4C).

We continued to analyze 438 overlapping genes, and the GO and KEGG enrichment results are shown in Figures 4C–F. Differentially expressed genes were mainly involved in BPs such as organelle fission and cilium movement (Figure 4D), CCs such as microtubule-associated complexes and spindles (Figure 4E), MFs such as microtubule motor activity and motor activity (Figure 4F), and KEGG enrichment pathways such as prostate cancer and the HIF-1 signaling pathway (see Figure 4G and Supplementary Table 4 for the detailed GO and KEGG enrichment pathways).

Finally, we analyzed the differences in the expression of 60 ferroptosis genes between different Fe.gene.clusters. The expression differences of 14 ferroptosis genes in the Fe.gene.cluster group were not significant, and the remaining genes were significantly different (see Figure 4H).

Univariate Cox hazard analysis was performed on 438 differentially expressed genes, and according to p < 0.05, 204 statistically significant prognosis-related genes were screened. These 204 genes were used to build a random forest model and screen out 13 characteristic genes as ferroptosis gene signatures. PCA and calculation of the Fe.score were performed (see Supplementary Table 2 for details).

Then, we explored the relationship between Fe.score, Fe.gene.clusters, and Fe.clusters. Figure 5A shows that most samples of Fe.cluster C and Fe.cluster D correspond to Fe.gene.cluster C, and most samples of Fe.gene.cluster C belong to the high-Fe.score group, which has a survival disadvantage. Most samples of Fe.gene.cluster A and Fe.gene.cluster B belong to the low-Fe.score group, which has a survival advantage.

For a more detailed analysis, among the Fe.gene.cluster groups, the Fe.score of Fe.gene.cluster A and Fe.gene.cluster B was significantly lower than that of Fe.gene.cluster C (Figure 5B). Among the Fe.cluster groups, the Fe.score of Fe.cluster A and Fe.cluster B is extremely significantly lower than that of Fe.cluster C, and the difference between Fe.cluster A and Fe.cluster D is not significant (Figure 5C).

We first explored the correlation between the Fe.score and cell function enrichment pathways. Figure 6A shows that the Fe.score is positively correlated with mismatch repair and DNA damage repair 1 and is negatively correlated with CD8 T effectors and immune checkpoints. Further statistical analysis shows that the enrichment scores of antigen processing machinery, CD8 T effector, EMT1, EMT2, and EMT3 all have extremely significant differences in the different Fe.score groups, and the Low Fe.score is higher than the High Fe.score (Figure 6B). In addition, for some new molecular indicators, Fe.score is significantly positively correlated with HRD, CNA_frac_altered, LOH_frac_altered, and LST, and the correlation coefficient is approximately 0.5; Fe.score is significantly negatively correlated with mutLoad_nonsilent, mutLoad_silent, SNV_Neoantigens, and Indel_Neoantigens, and the correlation coefficient is about −0.15 (Figures 6C–J).

The somatic mutation spectrum of the High Fe.score group is shown in Figure 7A, and the somatic mutation spectrum of the Low Fe.score group is shown in Figure 7B. The mutated genes and the types of mutations in the two groups were significantly different. We found that PTEN mutation, TP53 wild-type, KRAS mutation, and APC mutation patients had lower Fe.scores (Figure 7C). Figures 7D,E show the somatic mutation spectrum of the high- and low-Fe.score groups, and the G-score of the low-Fe.score group is significantly lower than that of the high-Fe.score group; the CNV frequency of the low-Fe.score group is also significantly lower than that of the high-Fe.score group (Figure 7F).

For clinical indicators, we found that patients with stage I, G1, premenopausal, younger age, and fewer children had lower Fe.scores (Figure 7G). The Kaplan–Meier survival curve of the Fe.score shows that there is a very significant difference in survival between the high-Fe.score and low-Fe.score groups, and the low-Fe.score group has a significant survival advantage (see Figure 7H).

We used GDSC to predict the difference in IC50 of the drugs cisplatin, erlotinib, rapamycin, docetaxel, and temsirolimus between the Fe.score groups. Through the nonparametric Wilcoxon rank sum test, the differences reached a very significant level (Figure 8A). We calculated and grouped the Fe.score of the bladder cancer dataset (IMvigor210) and found that the low-Fe.score group had a very significant survival advantage (Figure 8B). The Fe.score difference between the CR, PR, SD, and PD groups was not significant (Figure 8C), and CR/PR patients accounted for a relatively high proportion. Chi-square tests were performed on CR/PR and SD/PD between the high-Fe.score and low-Fe.score groups, and the difference was not significant (p-value = 0.075) (Figure 8D). Figure 8E shows that the Fe.score is better than the previous prognostic model of endometrial carcinoma (35), and the AUC value is higher.