We obtained the transcriptomic data and corresponding clinical characteristics of UCECs from the TCGA database, consisting of a total of 583 UCEC samples. We excluded normal samples, duplicate samples, and samples with missing survival times, leaving us with 539 samples for analysis. The transcriptomic profiles of TCGA-UCEC in the fragments per kilobase per million (FPKM) were transformed into transcripts per kilobase million (TPMs).

We applied the CIBERSORT package (http://cibersort.stanford.edu/) to obtain a fraction matrix of ICI using the TPM data, which estimated the abundances of 22 tumour infiltration immune cells (TIICs) [15]. The ESTIMATE algorithm was used to infer the tumour purity and calculate the immune score and stromal score. Additionally, single-sample gene set enrichment analysis (ssGSEA) was performed on UCEC samples using the GSVA package (http://www.bioconductor.org/) based on the expression levels of immune-related markers. We then used the “ConsensusClusterPlus” package (http://www.bioconductor.org/) to divide the UCEC samples into two subgroups. To ensure the stability of the classification, the analysis process was repeated 1000 times using the unsupervised clustering “Pam” method. We used the “pheatmap” R package (https://cran.rstudio.com) to draw the heatmap.

We divided UCEC samples into two ICI clusters based on the consensus clustering algorithm. Differentially expressed genes (DEGs) between the different ICI subgroups were identified by using the “limma” package. The significance cut-off was set as an adjusted p value 0.05 and absolute fold-change $>$1. Then, gene ontology (GO) annotation was performed by using the “ClusterProfiler” package (http://www.bioconductor.org/).

The “survival” package was used to perform univariate Cox regression analysis to identify DEGs related to the OS of UCEC patients. Then, we conducted multivariate Cox regression to construct an ICI-related prognostic risk model, and the risk score was calculated by the following formula for each UCEC patient. Risk score = ${\sum }_i^n\text{coef}\left(\text{gene}i\right)\times \text{expr}(\text{gene}i)$. Coef (gene$i$) was defined as the coefficient of gene i and expr (gene$i$) was defined as the expression of gene i.

According to the molecular classification method based on TCGA [16], we divided the 539 UCEC patients into four molecular subtypes: POLE, MSI-H, CNL and CNH. To investigate the correlation between the risk score and molecular subtypes, we generated box plots using GraphPad Prism v8.0.2 software (GraphPad Software, Inc., San Diego, CA, USA) and performed unpaired t tests to determine statistical significance. We used the same method to examine the relationship between risk score and clinical stage, survival status, and age groups.

We divided the 539 UCEC patients into two groups based on the median risk score. To explore the correlation between the risk score and TIME, we used the Wilcox test and drew box plots to demonstrate the difference in TIME between the low- and high-risk groups.

Research has demonstrated that the expression levels of certain key genes related to immune checkpoint blockade (ICB) may be correlated with the clinical outcomes of immune checkpoint inhibitor blockade therapy. To investigate the potential role of the risk score in immunotherapy, we analysed the correlation between the risk score and the expression levels of ICB-related genes. Here, we focus on the six key ICB genes, namely, PD-L1 (CD274), PD-1 (PDCD1) [17], PD-L2 (PDCD1LG2), TIM-3 (HAVCR2), CTLA-4 [18], and IDO1 [19], that have been reported in previous studies. Additionally, we further compared the correlation between nine inflammation-related genes (TNF, TBX2, PRF1, IFNG, CD8A, CLCL9, CXCL10, GZMA, and GZMB) [20] and the risk score. To evaluate the potential impact of immunotherapy predicted by the risk score, we also analysed the correlation between the immunophenoscore (IPS) and the risk score of UCEC patients. The IPS data for UCEC patients are available for download from The Cancer Immunome Atlas (TCIA).

To enhance the interpretability of the prognostic risk model and facilitate patient evaluation, we developed a nomogram using the “rms” R package, which integrates the risk score and age to estimate OS at 1, 3, and 5 years. The nomogram assigns a score to each factor based on their respective contribution to the outcome variable and sums them to obtain the total score. To validate the prognostic value of the nomogram, we used calibration and receiver operating characteristic (ROC) curves.

To reveal the intrinsic link between the TIME and somatic mutations, we conducted an analysis to assess the correlation between TMB and risk score. Somatic mutation data were obtained from the TCGA database, and the “maftools” package (http://www.bioconductor.org/) was utilized to calculate the number of somatic nonsynonymous point mutations in each sample, which corresponded to the TMB. Furthermore, we generated separate waterfall plots for the low- and high-risk groups to visually display the distribution of gene mutations among patients in each risk score group.

All statistical analyses were conducted using R version 4.1.0 (University of Auckland, Auckland, New Zealand). The fraction of 22 TIICs was calculated using the “CIBERSORT” package in R. The Wilcoxon test was used to compare differences between two groups. DEGs were identified using the “limma” package in R. Survival analysis was performed by the Kaplan-Meier (K-M) method, and comparisons between groups were made using the log-rank test. The “survivalROC” and “timeROC” packages were used to generate the ROC curve. GSEA was performed using GSEA software version 4.1.0 (Broad Institute, Inc., Cambridge, MA, USA).

To further analyse the TIME in UCEC, we calculated the level of 22 immune cells in a total of 539 samples. Based on the ICI profile, the 539 samples were classified into different subtypes using the “ConsensusClusterPlus” package. To ensure sample equilibrium and clustering stability, we chose k = 2 (Supplementary Fig. 1), resulting in the division of samples into two ICI patterns, ICI Cluster A (329 samples) and ICI Cluster B (210 samples). We drew a comprehensive heatmap to further describe and explore the relationship between ICI patterns and clinical phenotypes, including age and survival status (Fig. 2A). Furthermore, we conducted Kaplan‒Meier survival analysis to compare the prognosis of the two ICI subgroups and found that ICI Cluster B had a better prognosis than ICI Cluster A, although the difference was not statistically significant (p = 0.31, Fig. 2B).

Fig. 2.

The landscape of ICI in UCEC. (A) Heatmap of immune cells in UCEC patients. (B) K-M curves for the OS of UCEC patients in two ICI clusters. (C) The levels of immune cells, immune scores, and stromal scores in two ICI clusters. (D) Correlation of immune cells and immune scores. (E–G) OS analysis between patients with different infiltration levels of immune cells. (E) CD8+ T cells. (F) Regulatory T cells. (G) Monocytes. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns, not significant. OS, overall survival; NK, natural killer; K-M, Kaplan-Meier.

To further investigate the association between immune scores and infiltrating immune cells, we compared the TIME between the two clusters. Our analysis revealed that ICI Cluster B displayed higher levels of plasma cells, CD8+ T cells, activated state CD4+ T cells, follicular helper T cells, regulatory T cells, and M1 macrophages, as well as high stromal and immune scores. These findings are consistent with previous studies showing that the abundance of these immune cells significantly correlates with antitumour immune status [21]. Furthermore, we compared the correlations between immune cells and the immune score, stromal score, and ESTIMATE score (Fig. 2C,D). Among these immune cells, plasma cells, CD8+ T cells, activated-state CD4+ T cells, regulatory T cells, and M1 macrophages were significantly and positively correlated with immune scores. Therefore, we defined ICI Cluster B as immunoinflammatory. In contrast, most immune cell infiltrates were significantly lower in ICI Cluster A, except for resting-state CD4+ T cells, M0 macrophages, activated-state dendritic cells, higher levels of neutrophils, and lower stromal and immune scores. Resting-state CD4+ T cells, M0 macrophages, and activated-state dendritic cells were significantly and negatively correlated with immune scores. Therefore, ICI Cluster A can be considered the immune desert type.

To further assess the value of ICI subgroups in predicting prognosis, we conducted Kaplan‒Meier survival analysis based on the different ICI subgroups. The results showed that regulatory T cells and CD8+ T cells enriched in ICI Cluster B were associated with a longer survival time and a good prognosis (Fig. 2E,F). In contrast, a high level of monocytes, which were enriched in ICI Cluster A, showed a poor prognosis (Fig. 2G).

Furthermore, we compared the expression levels of six key immune checkpoint blocking genes (PD-1, CTLA-4, TIM-3, IDO1, PD-L1, and PD-L2) in the two ICI clusters (Fig. 3). Our analysis revealed that the expression levels of ICB-related genes were significantly higher in ICI Cluster B than in ICI Cluster A, suggesting that ICI cluster B may be more suitable for immunotherapy.

Fig. 3.

The expression levels of six ICB genes in two ICI clusters. (A) PD-1. (B) PD-L1. (C) PD-L2. (D) CTLA-4. (E) TIM-3. (F) IDO1. **p 0.01; ****p 0.0001. ICB, immune checkpoint blockade.

Using the limma package, we identified 29 DEGs between the two ICI clusters, which were defined as ICI-related genes (Table 1). We performed GO enrichment analysis with the ICI-related genes and found that most biological functions were related to the immune response, such as immunoglobulin production, humoral immune response, immunoglobulin complex, and immunoglobulin receptor binding (Supplementary Fig. 2). We then performed univariate Cox and multivariate Cox regression analyses to identify the ICI-related genes associated with OS and to establish a prognostic risk model. Finally, an optimal prognostic model was obtained: risk score = (0.9842 $\times$ expr (LINC01871)) + (1.0121 $\times$ expr (LINC01281)) + (0.9983 $\times$ expr (CXCL13)) + (0.9990 $\times$ expr (IGKJ5)). Kaplan‒Meier analysis revealed that a high level of the four genes was significantly associated with longer survival time (Fig. 4).

Fig. 4.

Kaplan‒Meier OS curves of four ICI-related genes. (A) LINC01871. (B) LINC01281. (C) CXCL13. (D) IGKJ5. Patients with high expression of ICI-related genes showed a good prognosis (red curves), while patients with low expression of ICI-related genes showed a poor prognosis (blue curves).

To further validate the prognostic role of risk scores, risk curves and scatter plots were constructed according to the survival status of each sample (Fig. 5A,B). The risk curves and scatter plots visualized the risk scores and survival status of TCGA-UCEC patients, illustrating that patient death was closely associated with high-risk scores. Patients were divided into high- and low-risk groups with the median risk score as the cut-off. The Kaplan‒Meier survival curve showed that the low-risk group had a longer survival time (p 0.001, Fig. 5C). Additionally, we analysed the correlation between the risk score and molecular subtypes and clinical stages. We found that the POLE and MSI-H subtypes had a significantly lower risk score, while the CNH subtype had a significantly higher risk score (Fig. 5D). These findings suggest that the POLE and MSI-H subtypes had a better prognosis, and the CNH subtype had a worse prognosis, which is consistent with a previous report [16]. The stage III and stage IV groups had a higher risk score than the stage I and stage II groups, implying that early diagnosis is beneficial to prognosis (Fig. 5E). Moreover, we analysed the correlation between the risk score and clinical phenotypes (survival status and age). We found that the deceased group had a significantly higher risk score (Fig. 5F), but there were no significant differences between age groups (Fig. 5G).

Fig. 5.

Evaluation of the ICI-related prognostic model. (A) Distribution of the risk score for each patient. (B) Survival status for each patient. Red points correspond to surviving patients, and blue points correspond to deceased patients. (C) Survival analysis between patients with different risk scores. (D) Risk score in different molecular subgroups. (E) Risk score in different clinical stages. (F) Risk score in the surviving and deceased groups. (G) Risk score in different age groups. *p 0.05; **p 0.01; ****p 0.0001; ns, not significant.

As the risk score was calculated by the ICI-related prognostic model, we further explored the differences in the tumour immune environment between the high- and low-risk groups. Our results showed that the risk score was negatively correlated with the level of various immune cells closely related to the antitumour immune response, including plasma cells, CD8+ T cells, activated state CD4+ T cells, follicular helper T cells, regulatory T cells, and M1 macrophages (Fig. 6A). Patients with low risk scores had higher stromal scores, immune scores, and ESTIMATE scores than patients with high risk scores (Fig. 6B). The ssGSEA results showed that both ICI and immune signal enrichment were significantly higher in the low-risk group than in the high-risk group (Fig. 6C). The heatmap exhibited the enrichment of the immune features of the patients in the high-risk and low-risk groups correspondingly (Fig. 6D). Based on these results, the high-risk and low-risk groups exhibited significantly different ICI phenotypes. Patients with low risk scores showed an immune-inflamed phenotype with a large number of immune cell infiltrates, while patients with high risk scores showed an immune desert phenotype with less immune infiltration. These findings indicated that the risk score was closely related to the TIME and may play an important role in immunotherapy.

Fig. 6.

Association between the risk score and TIME characteristics. (A) Immune cell infiltrating abundances in the low- and high-risk score groups. (B) Differences in stromal score, immune score, and ESTIMATE score between the low- and high-risk score groups. (C) Difference in immune functions between the low- and high-risk score groups. (D) Heatmap of immune cell infiltrating abundances between the low- and high-risk score groups. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns, not significant; TIME, tumour immune microenvironment.

In this study, we analysed ICB-related genes and inflammation-related genes to assess tolerance and immunoreactivity in patients with different risk scores (Fig. 7A). Almost all genes were significantly upregulated in patients with a low risk score, except for TBX2 and TNF. These results imply that patients with a low risk score have better tolerability and immunoreactivity. To further reveal the potential role of the risk score in the ICB treatment of UCEC patients, we examined the correlation between the six key ICB genes and the risk score and drew pie charts (Fig. 7B). We found that the expression of all six genes was significantly negatively correlated with the risk score (Fig. 7C–H). These results suggest that the risk score may play an important role in predicting the response to ICB treatment in UCEC patients.

Fig. 7.

Correlation of risk score with immunotherapy. (A) Differences in the expression levels of ICB-relevant genes (CD274, PDCD1LG2, PDCD1, CTLA4, IDO1, and HAVCR2) and inflammation-related genes (CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF) between the low- and high-risk score groups. (B–H) Correlation between the risk score and immune check inhibitors. (I–L) Difference in the IPS score between the low- and high-risk score groups. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns, not significant.

Two subtypes of IPS, IPS-PD1-positive and IPS-CTLA4-positive, were used as proxies for immunotherapy response in UCEC patients to assess the effect of immunotherapy predicted by risk scores (Fig. 7I–L). Our results showed that the IPS-CTLA4 and PD1 blockade score, IPS-CTLA4 blockade score, IPS-PD1 blockade score, and IPS score were significantly higher in the low-risk group. This indicates that patients with low risk scores are more suitable for immunotherapy.

We constructed a nomogram to predict the 1-, 3-, and 5-year OS of UCEC samples using the risk score and age (Fig. 8A). To validate the performance of the prognostic nomogram, we plotted calibration curves (Fig. 8B–D). Furthermore, we drew an ROC curve to further evaluate the constructed prognostic model, and the area under curve (AUC) values of 1-, 3-, and 5-year OS were 0.642, 0.658, and 0.716, respectively (Fig. 8E). These results indicate that the model had high sensitivity and specificity.

Fig. 8.

Evaluation of the risk score in prognosis. (A) A nomogram was constructed using age and risk score to predict the OS of UCEC patients. (B) 1-year nomogram calibration curves. (C) 3-year nomogram calibration curves. (D) 5-year nomogram calibration curves. (E) ROC curves for 1-, 3-, and 5-year survival. OS, overall survival; UCEC, uterine corpus endometrial carcinoma; ROC, receiver operating characteristic; AUC, area under curve.

It has been reported that high levels of TMB are associated with high infiltration of CD8+ T cells [22]. CD8+ T cells recognize tumour neoantigens and induce an immune response to destroy tumour cells. Therefore, TMB may serve as a prognostic factor for the effectiveness of the response to antitumour immunotherapy. In our study, we assessed the TMB levels in patients with both low and high risk scores and found that patients with a low risk score had a higher TMB (Fig. 9A). Additionally, our correlation analysis revealed a significant negative correlation between TMB and the risk score (Fig. 9B), indicating that a lower risk score may be associated with higher TMB and thus with a better response to antitumour immunotherapy.

Fig. 9.

Correlation between the risk score and TMB. (A) Difference in TMB between the low- and high-risk score groups. (B) Scatterplots depicting the negative correlation between risk score and TMB. (C) Survival analyses for patients with low and high TMB using K-M curves. (D) Survival analyses for patients stratified by both the risk score and TMB using K-M curves. The waterfall plot of tumour somatic mutations was constructed using the high-risk score (E) and low-risk score (F). ****p 0.0001. TMB, tumour mutation burden.

Then, we divided the patients into two groups based on their TMB levels. Kaplan‒Meier survival analysis showed that patients with high TMB had a greater probability of survival (Fig. 9C). To further assess the prognostic value of the risk score and TMB, we evaluated the synergistic effect of these two factors in predicting the survival of patients with UCEC (Fig. 9D). Patients with both a high TMB and low risk score exhibited the best survival probability, while patients with a low TMB and high risk score had the worst survival probability. These findings suggest that the risk score may serve as an independent prognostic factor and may be utilized to predict the clinical outcome of antitumour immunotherapy.

We also analysed the distribution of mutations between the low- and high-risk groups. The integrated somatic mutation landscape was used to visualize the mutation patterns of the top 20 most frequently mutated driver genes (Fig. 9E,F). The mutation profile of significantly mutated genes showed a higher mutation rate of TP53 (41% vs. 33%) in the high-risk group and a higher mutation rate of TTN (49% vs. 27%) and KMT2D (35% vs. 18%) in the low-risk group. These findings provide deeper insight into the intrinsic link between the risk score and somatic mutations in UCEC immunotherapy.