We downloaded the mRNA expression profiles and corresponding clinical data of the UCEC cohort from the TCGA database (June 2021, https://portal.gdc.cancer.gov/). The AS event data for UCEC were obtained from the https://bioinformatics.mdanderson.org/TCGASpliceSeq/(Ryan et al., 2016). Since these data are publically available, there was no requirement for approval by an ethics committee. We fully assessed the availability of clinical information. In our research, a few patients were excluded because they met the following criteria: 1) Epithelial neoplasm disease type, nos in TCGA; and 2) incomplete clinical data (e.g., age, grade, FIGO stage, and survival data). The percent spliced in (PSI) value can be used to quantify each AS event, which is the ratio of normalized reads indicating the presence of a transcript element versus the total normalized reads for that event, with a rating from 0 to 1. PSI = splice in/splice in+splice out. We screened the AS data for PSI value >0.75, representing the association between gene expression and AS events. We then merged the gene expression and clinical profiles using Perl (v5.30.0, https://www.perl.org/), establishing genomics and clinical databases for further research. A total of 524 patients with complete AS events and clinical data were included in our analysis. The clinical features of the patients are summarized in Table 1.

TCGA SpliceSeq is a database based on TCGA RNA sequencing (RNA-seq) data. Seven types of selective splicing events were analyzed, including Alternate Acceptor site (AA), Alternate Donor site (AD), Alternate Promoter (AP), Alternate Terminator (AT), Exon Skip (ES), Mutually Exclusive Exons (ME), and Retained Intron (RI). We analyzed the distributions of all encoded genes using the UpSetR package in each of the seven different types of AS events and survival-related AS events in UCEC.

Different AS events in genes led to diversity in outcomes, and changes in gene expression affected survival time. To further understand the prognostic value of AS events in UCEC patients, univariate Cox regression analysis with R package “survival” was performed to determine the survival-related different expressed alternative splicing (DEAS) events, including overall survival (OS)-related DEAS events. Next, the least absolute shrinkage and selection operator (LASSO) regression was applied to identify the final elimination of potential predictors with non-zero coefficients using the R package “glmnet”, which can avoid model overfitting to obtain a better fitting model. Furthermore, based on the results of LASSO Cox regression, predictive models were constructed using multivariate Cox regression analysis. Based on the PSI values and multivariate Cox analysis, we calculated the risk scores of each patient and obtained the corresponding coefficients, respectively. The following formula obtained the risk score: Risk   score = ∑ i = 0 n PSI × β i where β is the regression coefficient of the AS events. A total of 524 EC patients were divided into high- and low-risk groups bound by the median of risk score, and Kaplan-Meier survival analysis was performed to determine whether they had completely different prognoses. Furthermore, receiver operating characteristic (ROC) curves of 1, 3, and 5 years were generated using the survival ROC package in R to show the discrimination of the predictive signatures (Heagerty et al., 2000).

All clinical factors, including the risk score, age, FIGO stage, and grade, were incorporated to construct a nomogram to evaluate the probability of 1-, 3-, and 5-year OS of UCEC patients in the entire set. Validation of the nomogram was evaluated by calibration plot using the “rms” package. The calibration curve of the nomogram was plotted to assess the nomogram-predicted probabilities against the actual rates.

Normalized RNA expression data were used to infer the Immune Score using the estimate package (Yoshihara et al., 2013) and quantify the infiltrating proportions of 22 types of immune cells using the “CIBERSORT” package (Newman et al., 2015). The infiltrating percentage of 22 types of immune cells was equal to 100%. Single sample gene set enrichment analysis (ssGSEA) was used to quantify and classify the immunity stage based on immune-related gene (IRGs) sets (He et al., 2018). Next, 47 immune checkpoint genes were analyzed, and 16 of them that differed from the tumor and normal samples were screened. The differences between the 16 hub immune checkpoints among the high- and low-risk groups were analyzed, and the correlations between the six most important immune checkpoint genes (CD274, PDCD1, PDCD1LG2, CTLA4, HAVCR2, and IDO1) and the risk score were determined.

Using R package “limma” with the threshold of |log2FC|>1 and p < 0.05, the 11 genes (Table 2) involved in the model construction were analyzed to observe whether their expression differed between the UCEC and normal samples.

High- and low-expression groups of gub genes were obtained according to the gene expression. These were then used to analyze the difference in clinical indicators, including age, grade, and FIGO stage. Finally, the prognosis of the hub genes in the two groups was judged using the “survival” and “survminer” packages.

The ESTIMATE algorithm was applied to analyze the Stromal Score, Immune Score, ESTIMATE Score, and Tumor Purity based on the transcriptome profiles of UCEC to determine the effect ssGSEA grouping. We further compared the Stromal Score, Immune Score, ESTIMATE Score, and Tumor Purity in the high- and low-expression groups of hub genes using the Limma.R and ggpubr.R packages. The relationships between the copy number of hub genes and the quantity of six immune cells (B cell, myeloid DC cell, macrophage, neutrophil, T cell CD4⁺, and T cell CD8⁺) were evaluated using the Tumor Immune Estimation Resource (TIMER) database.

Splicing factors (SF) are protein factors involved in the splicing process of RNA precursors, which are closely related to the development and treatment of cancer (Dvinge et al., 2016; Obeng et al., 2019). A total of 404 SFs data downloaded from the SpliceAid2 database were used to analyze the correlation between the expression level of SFs and the PSI values of OS-associated AS events by R packages (BiocManager, limma). An absolute value of the correlation coefficient >0.6 and p < 0.001 were considered statistically significant. Finally, Cytoscape software (v3.7.2, https://cytoscape.org/) was used to visualize the potential SF-AS regulatory network.

All statistical analyses were performed using R version 4.1.0 (R packages: survival, survminer, UpSetR, glmnet, estimate, ggpubr, e1071, rms, preprocessCore, vioplot, ggExtra, GSVA, GSEABase, reshape2, pheatmap, corrplot, ggplot, ggplot2, and BiocManager). For all analyses, a two-tailed p < 0.05 was regarded as statistically significant if not noted.

A total of 524 UCEC patients were identified, and the baseline characteristics of these patients are summarized in Table 1. The mRNA splicing data included in this study contains 28,281 AS events in 8,141 genes. Given the possibility of multiple splicing modes for a single gene, we created UpSet plots to analyze interactive sets of seven types of AS events quantitatively. As shown in Figure 1A, a single gene could have up to five different splicing modes, and most genes had more than one AS event. Exon skip (ES) was the most frequent splice type among the seven AS types (34.4%), followed by an alternate terminator (AT) (27.5%) and alternate promoter (AP) (15.7%).

To explore the prognostic utility of an AS signature in EC, AS events associated with OS were identified by fitting univariate Cox proportional hazard regression models after merging the clinical data in the training cohort using Perl. In total, 1,108 AS events were determined with p < 0.05, including 633 high-risk survival-associated AS events (hazard ratio HR > 1) and 475 low-risk survival-associated AS events (HR < 1). The AS events can be counted through the UpSet plot. An UpSet plot was generated to visualize the intersecting sets between different genes and AS events. The bar charts on the left showed the number of genes with some kind of AS events. The upper bar charts showed the number of intersecting genes, indicating the number of genes with a certain type or types of AS events (Figure 1A). Figure 1A indicates that one gene might have more than one survival-associated AS event. It is noteworthy that the three highest frequency survival-associated AS events were still ES, AT, and AP in the UCEC cohort.

After conducting univariate regression analysis, LASSO regression was performed to select the optimal survival-related AS events to construct the prediction models to avoid model overfitting based on OS. First, 15 AS events were screened out by LASSO regression, and then the AS events with the same contribution were optimized (Figure 1B). Finally, an 11-AS event signature was identified as a predictor of survival in EC through the Cox proportional hazards regression model (Table 2). Besides, the minimum adjusted estimate of cross-validation prediction error was 0.020, and the p value of bootstrap (K = 1000) was 4.82E-305.

Kaplan-Meier curves and log-rank tests were plotted to explore the relationship between risk score and survival status. The survival probability of low-risk patients was higher than that of high-risk patients; in other words, high-risk patients had a higher mortality rate, as illustrated in Figure 1C (p < 0.001). We then applied ROC analysis to compare the predictive power of these prognostic models. The larger the area under the curve, the higher the accuracy of the model to predict the prognosis of patients. Figure 1D showed a robust and significantly improved performance; the areas under the ROC curve (AUC) in 1, 2, and 3 years were all greater than 0.800. The result illustrated that the accuracy of using the model to predict the 1-, 2- and 3-year survival rate of patients was relatively high. Moreover, the AUC of the risk score model predicting the 1-year survival rate was larger than that of the age, grade, and FIGO stage (Figure 1D). It means that the accuracy of predicting the 1-year survival rate of patients by the model is better than that of using other clinical parameters (age, grade, stage) to predict the prognosis.

Meanwhile, the risk scores of each UCEC patient were calculated, and all patients were divided into low- and high-risk groups bound by the median risk score. The distribution diagram of survival risk score (Figure 2A), survival status of EC patients (Figure 2B), and clustering heatmap of the PSI levels of eleven-AS markers (Figure 2C) are shown. The horizontal axis displays the patients’ order of risk score from low to high (Figure 2).

The calibration curve demonstrated that the predicted values are satisfactorily consistent in the prediction of the 1-, 3-, and 5-year OS because the red lines in three pictures are almost overlap with the 45° dashed lines (Figure 3B). The box charts in Figure 3E show whether there are differences in patients’ risk score among different clinical index (age, grade, stage). The risk score of patients >65 years old was higher than that of patients ≤65 years old. With the increase of grade and stage of UCEC, the risk score increased gradually (Figure 3B).

Univariate and multivariate Cox regression methods were used and combined with patient clinical characteristics (age, grade, and FIGO stage) to analyze whether the 11-AS event signature could be an independent predictor of survival in patients with UCEC. When the p values of univariate analysis and multivariate analysis were both less than 0.05, it was considered that the model could be an independent prognostic factor. As depicted in Figures 3C,D, the results showed that the risk score could still be used as a reliable and stable independent risk predictor in the UCEC cohort (p < 0.001; Figures 3C,D). We then constructed a predictive nomogram based on the multivariate analysis (Figure 3A) that included risk scores and clinical characteristics. The results demonstrated that the risk score had satisfactory diagnostic ability and clinical characteristics (p < 0.05).

First, the immune score in 29 types of infiltrating immune cells and immune function was assessed by the ssGSEA method (He et al., 2018). Figures 4B,C show the immune score differences of each immune cell in the low and high-risk score groups (Figures 4B,C). We further explored the impact of the risk score on the infiltration of 22 types of immune cells in the tumor microenvironment using the CIBERSORT algorithm. The landscape of 22 types of infiltrating immune cells in the low and high-risk score groups is shown in Figure 4A. Differential analysis results showed that eight types of immune cells [CD8 T cells, regulatory T cells (Tregs), activated natural killer (NK) cells, monocytes, M0 macrophages, M1 macrophages, resting dendritic cells, and activated dendritic cells] were significantly different between the two groups (p < 0.05).

The difference in the expression level of 47 immune checkpoint genes in the low- and high-risk score groups was assessed, and 26 genes were found to have significant differences (Figure 5A). Next, R packages (limma, corrplot, ggpubr, and ggExtra) were used to screen the risk scores related to the six most important checkpoint genes (CD274, PDCD1, PDCD1LG2, CTLA4, HAVCR2, and IDO1). Two immune checkpoint genes, PDCD1 and CTLA4, with negative correlation with risk score were identified (p < 0.001; Figure 5B). The scatter plot displaying the correlation of these two genes and the risk score were plotted separately. Although two of the correlation coefficients did not reach 0.3, the scatter plot showed a negative correlation (Supplementary Figure S1). At the same time, we can find that the expression of PDCD1 and CTLA4 in the high-risk score group was lower than that in the low-risk score group (Figure 5A).

The expressions of 11 genes (Table 2) were identified to analyze the difference between UCEC and normal cohorts by the “limma” package (with the threshold of |log2FC| > 1 and p < 0.05), and two genes, CYB561 (|logFC| = 1.1892, p < 0.001) and FOLH1 (|logFC| = 1.0862, p < 0.001), were extracted for further analysis. Next, we divided the tumor patients into high- and low-expression groups according to the optimal cut-offs in CYB561 and FOLH1 (4.67 in CYB561, 2.64 in FOLH1) for clinical prognostic analysis.

The correlations between the expression of the two hub genes and the clinicopathological parameters were evaluated using R packages (limma, survival, and survminer). CYB561 and FOLH1 expression levels were significantly associated with grade and FIGO stage (p < 0.05). The expression of CYB561 and FOLH1 decreased gradually with increases in the grade and stage. However, no notable association between the two genes and age was observed (p ≥ 0.05) (Figures 6A,B). Survival analysis revealed that the low expressions of both CYB561 and FOLH1 were associated with poor prognosis (p < 0.001) (Figure 6C).

The landscape of 22 types of infiltrating immune cells in the low- and high- CYB561 expression groups are shown in Figure 7A. The two groups differed between resting dendritic cells, neutrophils, activated memory CD4 T cells, resting memory CD4 T cells, and Tregs. Figure 7B shows the immune score difference of each immune cell in the two groups (Figure 7B). We investigated the association between CYB561 expression and the tumor-infiltrating immune cells in UCEC using the TIMER database. The results demonstrated that CYB561 expression was positively correlated with B cell, CD4⁺ T cells, and CD8⁺ T cells, and was negatively correlated with myeloid dendritic cells, macrophages, and neutrophils (p < 0.05) (Supplementary Figure S2). However, we found no strong correlation between immune cell infiltration and CYB561 expression. Given that the risk score was related to tumor immunity, we finally appraised the correlation between the gene signature and the expression of immune checkpoints. Figure 7C shows the 23 immune checkpoints with differential expression in the low- and high-CYB561 expression groups. The gene expression of CD40LG, TNFSF14, TNFRSF14, CD276, VTCN1, HHLA2, TNFSF15, LGALS9, and CD44 was lower in the low-CYB561 expression groups (Figure 7C).

The landscape of 22 types of infiltrating immune cells in the low- and high- FOLH1 expression groups are shown in Figure 8A. Resting memory CD4 T cells, gamma delta T cells, resting NK cells, resting dendritic cells, activated dendritic cells, and neutrophils were different in the two groups. Figure 8B shows the immune score difference of each immune cell in the two groups. We further investigated the association between FOLH1 expression and the tumor-infiltrating immune cells in UCEC. The results showed that FOLH1 expression was positively correlated with macrophages, CD4⁺ T cells, and CD8⁺ T cells, and was negatively correlated with B cells, myeloid dendritic cells, and neutrophils (p < 0.05) (Supplementary Figure S3). However, we also found no strong correlation between immune cell infiltration and FOLH1 expression. Figure 8C shows the immune checkpoints with differential expression in the low- and high-FOLH1 expression groups.

First, the violin plot assessed the differences in Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score between the two groups, calculated using the ESTIMATE algorithm (Figure 9A). ESTIMATE Score and Immune Score were higher in the low-risk score group, while Tumor Purity in the low-risk score group was lower than that in the high-risk score group (p < 0.05). The Stromal Score showed no difference (p ≥ 0.05).

In order to determine the effectiveness of the grouping strategy between the low- and high-expression groups of CYB561 and FOLH1, the ESTIMATE method was applied to evaluate Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score. Compared with the high-CYB561 expression group, the low expression group had a higher Stromal Score (p < 0.05) (Figure 9B). The other parameters had no differences between the two groups in CYB561 and FOLH1 (p ≥ 0.05) (Figures 9B,C).

The correlations between the CYB561 copy number and six immune cells were also analyzed using the TIMER database. The number of cells was found to decrease with the increase in the gene copy number in myeloid DC cells and CD8⁺ T cells, and was found to decrease with the decrease in the gene copy number in myeloid DC cells and CD8⁺ T cells (p < 0.05) (Supplementary Figure S4). Next, we analyzed the correlations between the FOLH1 copy number and six immune cells. The number of cells was found to decrease with the increase in the gene copy number in myeloid DC cells, macrophage, CD8⁺ T cells, and the cells number was found to decrease with the reduce in the gene copy number in myeloid DC cells and CD8⁺ T cells (p < 0.05) (Supplementary Figure S5). Notably, we observed that the change in gene copy numbers in the two hub genes led to the deletion number of both myeloid DC cells and CD8⁺ T cells.

Thirty-six SFs (blue) were found to be significantly related to 19 survival-associated AS events, consisting of nine low-risk AS events (RPS24|12297|AA, ZNF169|86927|AT, MAT2A|54280|AT, C22orf39|61055|AT, VDAC1|73335|AP, FAM72A|9577|AT, PPP2R5C|29315|AT, TMEM33|69133|RI, SYTL1|1321|AP; purple) and 10 high-risk AS events (LEPROTL1|83274|AT, ANKHD1|73652|AT, MAT2A|54281|AT, BUB3|13390|AA, ATP5D|46401|RI, FAM72A|9578|AT, VDAC1|73334|AP, MARVELD3|37467|AT, NDUFA3|51776|AT,SYTL1|1320|AP; red). The majority of low-risk AS events were negatively correlated with SF expression (blue lines), and all of the high-risk AS events were positively correlated with SF expression (red lines) (Figure 10).