Identification of Alternative Splicing-Related Genes CYB561 and FOLH1 in the Tumor-Immune Microenvironment for Endometrial Cancer Based on TCGA Data Analysis

Background: Advanced and recurrent endometrial cancer EC remains controversial. Immunotherapy will play a landmark role in cancer treatment, and alternative splicing (AS) of messenger RNA (mRNA) may offer the potential of a broadened target space. Methods: We downloaded the clinical information and mRNA expression profiles from The Cancer Genome Atlas (TCGA) database. Hub genes were extracted from 11 AS-related genes to analyze the correlation between clinical parameters and the tumor-immune microenvironment. We also analyzed the correlations between the copy numbers, gene expressions of hub genes, and immune cells. The correlation between the risk score and the six most important checkpoint genes was also investigated. The ESTIMATE algorithm was finally performed on each EC sample based on the high- and low-risk groups. Results: The risk score was a reliable and stable independent risk predictor in the Uterine Corpus Endometrial Carcinoma (UCEC) cohort. CYB561|42921|AP and FOLH1|15817|ES were extracted. The expression of CYB561 and FOLH1 decreased gradually with the increased grade and International Federation of Gynecology and Obstetrics (FIGO) stage (p < 0.05). Gene copy number changes in CYB561 and FOLH1 led to the deletion number of myeloid DC cells and T cell CD8+. Low expression of both CYB561 and FOLH1 was associated with poor prognosis (p < 0.001). The checkpoint genes, CTLA-4 and PDCD1, exhibited a negative correlation with the risk score of AS in UCEC. Conclusion: AS-related gene signatures were related to the immune-tumor microenvironment and prognosis. These outcomes were significant for studying EC’s immune-related mechanisms and exploring novel prognostic predictors and precise therapy methods.


INTRODUCTION
Uterine corpus endometrial carcinoma (UCEC) is one of the three major malignancies of the female reproductive system (Bray et al., 2018). Although most patients are diagnosed at an early stage and have a good prognosis, there are still some patients who are at an advanced stage at first diagnosis or have recurrence and metastasis after treatment, with a 5-year survival rate of only 20-26% (Morice et al., 2016). Therefore, further studies on therapeutic monitoring and prognostic assessment of UCEC are crucial for both clinicians and patients. In 2013, tumor immunotherapy was regarded as an important scientific breakthrough and suggested that immunotherapy will play a landmark role in the field of cancer treatment (Couzin-Frankel, 2013). In recent years, immune checkpoint inhibitors have made breakthrough progress and have been written into the treatment guidelines for endometrial cancer (EC).
Alternative splicing (AS) is a universal mechanism to produce mRNA isomers using a limited set of genes, resulting in structurally and functionally different protein isoforms, modifying more than 95% of human genes (Gilbert, 1978;Buratti et al., 2006). Studies have shown that aberrant AS is closely associated with the occurrence, development, metastasis, and drug resistance of various cancers (de Necochea-Campion et al., 2016;Climente-González et al., 2017;Wang and Lee, 2018;Bonnal et al., 2020). AS has also become a hot topic in tumor immunotherapy and attracted the attention of researchers. Several forms of mRNA processing are dysregulated in cancer and offer promise concerning immunotherapy target expansion (Venables, 2004;Baralle and Giudice, 2017;Kahles et al., 2018). Although significant progress has been achieved in expanding the immunotherapy target space using tumor-specific mRNA processing events, much work is still needed (Frankiw et al., 2019).
What's more, splice factors might play a vital carcinogenic role in EC (Dou et al., 2020a;Li et al., 2020;Popli et al., 2020). By analyzing the whole genome of AS events in EC, studies have found several candidate splicing factors that may become therapeutic targets and predict patients' prognosis by constructing gene signatures (Wang C. et al., 2019;, which further demonstrates the importance of AS events in EC. A study found that with the increase of the ESTIMATE score and the infiltration of immune cells in UCEC patients, the prognosis would be better (Liu et al., 2021); however, further discussion was lacking.
Given the importance of immunotherapy in UCEC, characterization of immune infiltrating features is essential for further understanding the oncogenesis of UCEC and the development of a novel prognostic signature and therapeutic response classifier. In the present study, whole-genome analysis and prognostic model construction were firstly used to determine prognosis-related genes involved in the AS prognostic model. The characteristics of two AS-related genes in the tumor-immune microenvironment were then analyzed. Finally, the correlation between the six most important immune checkpoint genes and the risk score was also investigated. Our research provides a more comprehensive insight into precise immunotherapy for UCEC.

Data Acquisition and Curation Processing
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.

Screening for Prognostic AS Events in UCEC
TCGA SpliceSeq is a database based on TCGA RNA sequencing (RNA-seq) data. Seven types of selective splicing events were

Construction of Prognostic Models and Survival Analysis
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 n i 0 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).

Establishment and Validation of a Predictive Nomogram
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.
Immune Score Estimate and Immune Cell Infiltrating Proportion Inference 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.

Extraction of AS-Related DEGs in UCEC Samples
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.

Integration of AS-Related DEGs With Clinical Characteristics and Prognosis
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.

Analysis of the Relationship Between Stromal/Immune Scores and AS-Related DEGs in the EC Immune Microenvironment
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 highand 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.

Construction of a Potential SF-AS Regulatory Network
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 OSassociated 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.

Overview of AS Events in TCGA UCEC Cohort
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%).

Prognostic Index Models Featured by AS Events for UCEC
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 lowrisk 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 survivalassociated 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 crossvalidation 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  Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 770569 5 markers ( Figure 2C) are shown. The horizontal axis displays the patients' order of risk score from low to high (Figure 2).

Construction and Evaluation of the Nomogram
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 Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 770569 6 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

The Risk Score and AS Events Are Associated With the Infiltration of Immune Cells in the UCEC Microenvironment
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) Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 770569 9 cells, monocytes, M0 macrophages, M1 macrophages, resting dendritic cells, and activated dendritic cells] were significantly different between the two groups (p < 0.05).

The Risk Score is Associated With the Key Immune Checkpoint Genes in the UCEC Tumor-Immune Microenvironment
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).

Extraction of IRGs Depending on AS Events and Their Correlation With Clinical Parameters
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).

Associations Between CYB561 Expression and Immune Cell Infiltration
The landscape of 22 types of infiltrating immune cells in the lowand 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).

Correlations Between FOLH1 Expression and Immune Cell Infiltration
The landscape of 22 types of infiltrating immune cells in the lowand 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 lowand high-FOLH1 expression groups.
The Tumor Purity, ESTIMATE Score, Immune Score, Stromal Score and Between the Low-and High-Expression Groups of CYB561 and FOLH1 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 Quantity of Six Immune Cells in Gene Copy Number of CYB561 and FOLH1
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.

DISCUSSION
Dysregulation of AS can affect essential biological processes and thus drive disease-associated pathophysiology (Gamazon and Stranger, 2014). Emerging data have demonstrated that aberrant AS events are closely associated with cancer progression, metastasis, therapeutic resistance, and other oncogenic processes (Climente-González et al., 2017). Cancer cells have general and cancer type-specific and subtype-specific alterations in the splicing process, which can have prognostic value and contribute to every hallmark of cancer progression, including the cancer-immune responses (Bonnal et al., 2020). Moreover, substantial preclinical work has identified a variety of FIGURE 9 | Construction and verification of the immune-related UCEC subgroups. (A) The violin plot shows the differences in Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score between the high-and low-risk score groups, calculated using the ESTIMATE algorithm. Comparison of Stromal Scores, Immune Scores, ESTIMATE Scores, and Tumor Purity between two groups of CYB561 (B) and FOLH1(C).
small molecule compounds and genetic and other approaches to target the spliceosome or its products with potential therapeutic effects (Bonnal et al., 2020). Therefore, it is of great importance to further study the characteristics of AS in the immune microenvironment for UCEC immunotherapy. In recent years, the relationship between AS and UCEC has been studied. In endometrial cancer, AS of vascular endothelial growth factor A (VEGF-A) is regulated by RBM10 (Dou et al., 2020a). Popli et al. found that SF3B1 plays a crucial oncogenic role in the tumorigenesis of EC and hence may support the development of SF3B1 inhibitors to treat this disease (Dou et al., 2020a). XQ et al. showed that miR-335 modulates Numb AS via targeting RBM10 in EC (Dou et al., 2020b).
We extracted IRGs depending on AS events and examined their correlation with clinical parameters. Finally, two AS-related genes, CYB561|42921|AP and FOLH1|15817|ES, were extracted from the 11 genes involved in the AS prognostic model. CYB561 encodes the protein CYB561, named as such because of its optical absorbance at 561 nm CYB561 is a heme-containing enzyme that is necessary for the continuous regeneration of semidehydroascorbate to ascorbate inside chromaffin granules and neuropeptide secretory vesicles (van den Berg et al., 2018). It is widely expressed in the adrenal glands, prostate, and 23 other tissues, including the endometrium. However, data on the role of the CYB561 gene in human cancers are very limited. A metaanalysis showed that low mRNA expression of CYB561 was prognostic of a poor outcome in ovarian cancer (Willis et al., 2016). CYB561 serves as a potential prognostic biomarker and target for breast cancer (Yang et al., 2021). We found the expression of CYB561 decreased gradually with the increased grade and FIGO stage which indicated a lower survival probability (p < 0.001) in UCEC. Our results were consistent with the previous ones. Besides, alternate promoter of CYB561 was associated with the OS of UCEC patients (p = 0.0003) in our research. The changes in transcription is regarded as a defining feature of cancer. Most human protein-coding genes are regulated by multiple, distinct promoters, suggesting that the selection of promoter is closely related to the expression of target gene (Demircioğlu et al., 2019). How the AP contributes the low expression of CYB561 in endometrial carcinoma remains to be further explored.
FOLH1 is also known as prostate-specific membrane antigen (PSMA), which encodes a transmembrane glycoprotein that acts as a glutamate carboxypeptidase on different alternative substrates. In the prostate, this protein is up-regulated in cancerous cells and is used as an effective diagnostic and prognostic indicator of prostate cancer (Date et al., 2017). PSMA is highly and specifically expressed in the neovasculature of ovarian, endometrial, and cervical squamous carcinomas (Wernicke et al., 2017 which is consistent with our findings. Their research indicated that the loss of PSMA expression can be considered a prognostic marker in patients with endometrial adenocarcinoma and could be due to epigenetic silencing (Mhawech-Fauceglia et al., 2008). FOLH1 likely arose from a duplication event of a nearby chromosomal region. Alternative splicing gives rise to multiple transcript variants encoding several different isoforms (Watt et al., 2001;Zink et al., 2020). Our research found ES in FOLH1 was associated with the OS of UCEC patients (p = 0.0000). What's more, ES was also the most frequent splice type among the seven AS types (34.4%) in UCEC. If the normal exon can be restored into the exon of ES occurred, it will bring hope to the treatment of many diseases (Verhaart and Aartsma-Rus, 2019). However, the mechanism of ES in FOLH1 leading to a high stage and poor prognosis of UCEC is unknown. Next, two immune checkpoint genes, Cytotoxic Lymphocyte Antigen 4 (CTLA-4) and Programmed Cell Death 1 (PDCD1), showed negative correlations with the risk score of AS in UCEC. CTLA-4 is expressed on the surface of naive effector T cells and Tregs (Billeskov et al., 2017;Menéndez-Menéndez et al., 2019). Based on its role as a negative regulator of T cell activation, CTLA-4 has become an attractive target for therapies aiming to enhance the effector activity of T lymphocytes. The first targeted drug for CTLA-4, ipilimumab, was approved by the Food and Drug Administration (FDA) in 2011 to treat melanoma (Lipson and Drake, 2011). At present, both nivolumab and ipilimumab are undergoing phase II clinical trials in UCEC (Grywalska et al., 2019). In our study, the CTLA-4 gene expression, the number and immune score of Tregs all decreased in the high-risk score group, which predicted a worse prognosis. Therefore, we supposed that a high-risk score of AS might be related to the decreased immune activity of Treg cells and the low expression of CTLA-4. It is possible that the targeted regulation of AS can improve the immune activity of Treg cells and increase the expression of CTLA-4, which may be valuable in improving the survival rate of UCEC patients, although further confirmation is needed. PDCD1, also known as PD-1, functions primarily in peripheral tissues. It is expressed on the surface of activated T cells, Tregs, activated B cells, and NK cells (Page et al., 2014). In 2014, the first FDA-approved immune checkpoint inhibitor targeting PD-1 was nivolumab (Grywalska et al., 2019). During the 2015 annual meeting of the Society of Gynecologic Oncology, Herzog et al. reported that the highest PD-1 expression rates among studied cancer types were in EC (75.2%) (Page et al., 2014). We found that PDCD1 expression was suppressed in the high-risk score group, and the 5-year survival rate was lower than that in the low-risk score group. It has been confirmed that there are variable splicing events in the PD1 gene (Nielsen et al., 2005;Wang et al., 2021). Another research considered the AS events in PD-1 may be a novel source for diagnostic and therapeutic target on celiac disease (Ponce de León et al., 2021). Why the high-risk AS events in PDCD1 lead to a worse prognosis of endometrial cancer needs further study.
The current study also has several limitations that should be noted. Firstly, this study is based on bioinformatics analysis, and there are no recruited cohorts for prognostic verification. Secondly, the values of the two-gene signatures for immunotherapeutic drugs prediction have not been verified in patient cohorts.

CONCLUSION
This study assessed the heterogeneity of tumor-infiltrating immune cells in UCEC and identified two AS-related genes, CYB561 and FOLH1, from the 11 genes involved in the AS prognostic model. Two immune checkpoint genes, CTLA4 and PDCD1, were negatively correlated with the risk score. The outcomes of this study are significant for investigating the immune-related mechanisms of tumor progression and exploring novel prognostic predictors and precise therapy methods.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.