An 11-GRG Predicting Model for UCEC Prognosis

doi:10.21203/rs.3.rs-3235391/v1

Background

Endometrial cancer (EC) is indeed one of the most prevalent gynecological malignancies. Further research is required to gain insights into the molecular pathways involved in EC tumorigenesis and to develop more accurate prognostic prediction methods.

Method

Data from the Cancer Genome Atlas( TCGA) database was used and validated using two GEO datasets, specifically GSE6008 and GSE17025. Various bioinformatics analyses were performed, including the Least Absolute Shrinkage and Selection Operator regression (lasso) regression, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis Gene Set Variation Analysis (GSVA), gene set enrichment analysis (GSEA), protein-protein interaction (PPI) network analysis, COX regression, calibration curves, and decision curve analysis (DCA). These analyses aimed to explore the associations and functions of GRGs in EC development, as well as develop a predictive model for prognosis assessment.

Results

There were 11 GRGs identified as significantly associated with EC by the Lasso regression, GSEA, and PPI. The further analysis identifies 61 miRNA molecules, 24 RBPs, 139 TFs, and 21 potential drugs or molecular compounds that might have links to these 11 key genes. These 11 GRGs were made into a Cox regression prediction model, among which the PGK2 shows significantly higher utility than other variables. Calibration analysis and DCA indicate that the clinical predictive performance of this 11-GRGs multivariate Cox regression model is highest at 5 years, followed by 3 years and 1 year. There were 6 genes (GPI, HK1, NUP188, PDHA1, PDHA2, PGK2) that exhibited higher predictive accuracy in time-dependent ROC curves.

Conclusion

The highly enriched GRGs that have been identified might provide a new understanding of the development of EC and its treatment. Moreover, the 11-GRGs model that has been constructed holds significant clinical implications for evaluating prognosis and providing specific therapy guidance to individuals with EC. Among the variables in the 11-GRG model, PGK2 demonstrates notably higher usefulness, highlighting its potential clinical value in EC.

Endometrial cancer

glycolysis-related genes

prognosis

PGK2

Warburg effect

Uterine corpus cancer is a prevalent type of gynecological tumor in developed countries such as the United States and Europe, causing approximately 417,000 new cases and 97,000 deaths worldwide in 2020.¹ Among these cases, uterine corpus endometrial carcinoma (EC) is the most frequently diagnosed form of uterine cancer, originating from the uterine epithelium.² Early-stage diagnoses of EC are common and have led to improved survival rates due to effective treatment. However, advanced EC remains largely incurable ^{1, 3}. Therefore, it is crucial to identify new molecular players that could serve as biomarkers for diagnosis, prognosis, and potential therapeutic targets in EC.

Unlike normal cells that predominantly use mitochondrial oxidative phosphorylation for energy production, cancer cells exhibit a metabolic shift towards glycolysis as their primary energy source, known as the "Warburg effect.”^{4, 5} Building upon this understanding, in this study, it was hypothesized that glycolysis-related genes (GRGs) could serve as markers for the progression of EC and patient prognosis.

Therefore, the primary objective of the study was to develop a prediction model based on GRGs that could accurately predict clinical prognosis and identify potential factors associated with key GRGs involved in EC progression and prognosis. The study aimed to provide novel insights into the field of precision medicine for EC.

1 Data acquisition

We downloaded the endometrial carcinoma (EC) dataset (TCGA-UCEC) from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) using the TCGAbiolinks package⁶. After excluding normal samples and samples lacking clinical information, we obtained a total of 552 EC samples with clinical information for analysis. This dataset includes 522 samples from endometrial carcinoma patients (group: EC). The count sequencing data of the TCGA-UCEC EC dataset was also normalized to Fragments Per Kilobase per Million (FPKM) format, while the corresponding clinical data for the samples were obtained from the UCSC Xena database (http://genome.ucsc.edu)⁷.

In addition, we also downloaded expression profile datasets of endometrial carcinoma (EC) patients from the GEO database (https://www.ncbi.nlm.nih.gov/geo/)⁸ using the R package GEOquery⁹. The datasets include GSE6008^10-14 and GSE17025^{15, 16}. The GSE6008 dataset consists of 37 samples of endometrioid carcinomas, 41 samples of serous carcinomas, 13 samples of mucinous carcinomas, 8 samples of clear cell carcinomas, and 4 normal samples. The data was obtained through chip-based expression profiling (Expression profiling by array). In total, there are 103 samples in this dataset, and the data platform used is GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. The GSE17025 dataset includes 79 samples of endometrioid carcinomas, 12 samples of serous carcinomas, and 12 samples of atrophic endometrium from postmenopausal women. The data was obtained through chip-based expression profiling (Expression profiling by array). There are a total of 103 samples in this dataset, and the data platform used is GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.

The probe annotation of the dataset uses the chip GPL platform files. The expression profile data of endometrioid carcinomas, serous carcinomas, mucinous carcinomas, and clear cell carcinomas samples from the dataset GSE6008 are included for subsequent analysis, with a total of 99 samples. The expression profile data of endometrioid carcinomas and serous carcinomas samples from the dataset GSE17025 are also included for subsequent analysis, with a total of 91 samples. The datasets GSE6008 and GSE17025 are used as validation sets for further verification. For detailed information about the datasets, please refer to Table 1.

Glycolysis is an important process in biological metabolism. We collected glycolysis-related genes (GRGs) from the GeneCards database (https://www.genecards.org/)¹⁷ using "Glycolysis" as the search keyword. Only "Protein Coding" genes with a relevance score greater than 0.500 were retained as glycolysis-related genes (NRGs). Additionally, we obtained a set of glycolysis-related genes (GRGs) from published literature^18-20 on the PubMed website using "Glycolysis" as the search keyword. The GeneCards database contains a total of 1,514 glycolysis-related genes (GRGs), and the published literature on PubMed contains 291 glycolysis-related genes (GRGs). Furthermore, by searching "Glycolysis" in the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/)²¹, we obtained 21 reference gene sets related to glycolysis, including BIOCARTA_GLYCOLYSIS_PATHWAY, which together comprise 733 glycolysis-related genes (GRGs). After merging and taking the intersection, we obtained a total of 113 glycolysis-related genes (GRGs). Please refer to Table 2 for specific gene names.

2 prognosis model for GRGs

To obtain a prognostic model for endometrial cancer Glycolysis-Related Genes (GRGs), we performed Least Absolute Shrinkage and Selection Operator (Lasso) regression on the TCGA-UCEC dataset using 10-fold cross-validation and a p-value of 0.05. We ran 1000 cycles to prevent overfitting. Lasso regression is commonly used to build prognostic models. It builds upon linear regression by adding a penalty term (lambda times the absolute value of the coefficients) to reduce overfitting and improve the model's generalization ability. We visualized the results of the Lasso regression and further explained the grouping of each sample based on risk scores, survival outcomes, and the molecular expression patterns of glycolysis-related prognostic genes within each group using a risk factor plot.

3 The risk score of the model is grouped high and low

To identify differentially expressed genes associated with the high and low-scoring groups in the Lasso regression prognostic model, we performed differential analysis on the gene expression profiles of the TCGA-UCEC dataset using the limma package [8]. We selected genes with |logFC| > 1 and adjusted P-value < 0.05 as the criteria for identifying differentially expressed genes for further study. Genes with logFC > 1 and adjusted P-value < 0.05 were considered up-regulated genes, indicating an increase in expression. Genes with logFC < -1 and adjusted P-value < 0.05 were considered down-regulated genes, indicating a decrease in expression. The results of the differential analysis will be used for subsequent analysis.

4 GO and KEGG analysis

Gene ontology (GO) ²² analysis is a commonly used method for large-scale functional enrichment studies. It includes biological processes (BP), molecular functions (MF), and cellular components (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG)²³ is a widely used database that stores information related to genomes, biological pathways, diseases, drugs, and more. We used the R package clusterProfiler²⁴ to perform GO annotation analysis on differentially expressed genes (DEGs). The criteria for selecting entries were P.adj < 0.05 and a false discovery rate (FDR) value (q-value) < 0.05, which were considered statistically significant. The Benjamini-Hochberg (BH) method was used for P-value correction.

5 GSEA and GSVA analysis between high and low-risk groups

Gene Set Enrichment Analysis (GSEA)²⁵ is used to evaluate the distribution trend of a pre-defined gene set within a ranked gene list associated with a phenotype, to assess its contribution to the phenotype. In this study, we first divided all genes in the TCGA-UCEC dataset GRGs prognostic model into two groups based on the ranking of gene logFC values: high-risk and low-risk groups. Then, we used the clusterProfiler package to perform enrichment analysis on all genes between the high-risk and low-risk groups in the TCGA-UCEC dataset GRGs prognostic model. The parameters used in this GSEA enrichment analysis are as follows: seed = 2023, permutation times = 1000, minimum number of genes per gene set = 10, maximum number of genes per gene set = 500, and the P-value correction method is Benjamini-Hochberg (BH). We obtained the gene set c2.all.v2022.1.Hs.symbols.gmt from the Molecular Signatures Database (MSigDB). The criteria for significant enrichment were P.adj < 0.05 and a false discovery rate (FDR) value (q-value) < 0.25.

Gene Set Variation Analysis (GSVA)²⁶, also known as gene set variation analysis, is a non-parametric and unsupervised analytical method. It primarily converts the gene expression matrix between different samples into a gene set expression matrix to assess the enrichment of gene sets in chip-based transcriptomes. This allows the evaluation of whether different pathways are enriched among different samples.

We conducted GSVA analysis on the gene expression matrix of high-risk and low-risk groups in the TCGA-UCEC dataset using the reference gene set "h.all.v7.4.symbols.gmt" obtained from the Molecular Signatures Database (MSigDB) [https://www.gsea-msigdb.org/]. Furthermore, we calculated the functional differences of the enriched pathways in the GSVA enrichment analysis results between the high-risk and low-risk expression groups in the TCGA-UCEC dataset. A P.adj < 0.05 was considered statistically significant.

6 Protein-protein interaction network， PPI network

Protein-protein interaction network (PPI network) is composed of individual proteins interacting with each other to participate in various biological processes such as signal transduction, gene expression regulation, energy and substance metabolism, and cell cycle control. The STRING database²⁷ is a database that searches for known and predicted protein interactions. In this study, we utilized the STRING database to construct a protein-protein interaction network (PPI network) using the GRGs and set the minimum required interaction score at low confidence (0.700). We visualized the PPI network model using Cytoscape. All nodes in the PPI network were considered key genes (hub genes) for subsequent research.

7 Construct mRNA-miRNA, mRNA-RBP, mRNA-TF, and mRNA-drugs interaction networks

The miRDB database²⁸ is used for miRNA target gene prediction and functional annotation. We used the miRDB database to predict the interaction between miRNAs and key genes (hub genes), and then we compiled the data with a Target Score > 80 from the miRDB database to construct an mRNA-miRNA interaction network.

ENCORI²⁹database (https://starbase.sysu.edu.cn/) is the 3.0 version of the starBase database. The interactions in the ENCORI database, including miRNA-ncRNA, miRNA-mRNA, ncRNA-RNA, RNA-RNA, RBP-ncRNA, and RBP-mRNA, are based on CLIP-seq and degradome sequencing data mining, specifically for plants. It provides various visualization interfaces for exploring microRNA targets. Additionally, we use the ENCORI database to predict RNA binding proteins (RBPs) that interact with key genes (hub genes), and we filter and plot mRNA-RBP interaction networks using the criteria of clusterNum > 2 and clipExpNum > 2.

The CHIPBase database (version 3.0)³⁰ (https://rna.sysu.edu.cn/chipbase/) identifies thousands of binding motif matrices and their binding sites from ChIP-seq data of DNA-binding proteins. It also predicts the transcriptional regulatory relationships between millions of transcription factors (TFs) and genes. The hTFtarget database³¹ (http://bioinfo.life.hust.edu.cn/hTFtarget.) is a comprehensive database that includes data on human transcription factors (TFs) and control target regulation. We used the CHIPBase database (version 3.0) and the hTFtarget database to search for transcription factors (TFs) that interact with the hub genes and visualize them using Cytoscape software.

In addition, we also use the Comparative Toxicogenomics Database (CTD)³² (http://ctdbase.org/) to predict potential drugs or small molecular compounds that interact with key genes (hub genes). Then, we visualize the mRNA-miRNA, mRNA-RBP, mRNA-TF, and mRNA-drugs interaction networks using Cytoscape software.

8 Prognostic clinical correlation analysis.

To obtain a multivariable Cox model for the high- and low-risk groups in the TCGA-UCEC dataset, we first conducted univariate Cox regression analysis on the key genes. Then, we included all variables for multivariable Cox regression analysis and constructed a multivariable Cox model. Based on the results of the multivariable Cox regression analysis, we used the R package "rms" to construct a column line graph (nomogram). The column line graph (nomogram)³³, also known as a nomogram graph, is a graphical representation of the functional relationship between multiple independent variables using a cluster of non-intersecting line segments on a Cartesian coordinate system. It represents the variables within the multivariable regression model by assigning scores based on a specific scale, and ultimately calculates a total score to predict the probability of an event occurring.

Finally, we evaluated the accuracy and discriminatory power of the multivariable Cox model using a calibration curve. A calibration curve (Calibration Plot)³⁴ is used to assess the predictive performance of a model by plotting the fit between the actual probabilities and the predicted probabilities under different scenarios. It is mainly used to analyze the goodness-of-fit of the Cox regression model to the observed outcomes. Additionally, we utilized Decision Curve Analysis (DCA)³⁵, which is a simple method for evaluating clinical prediction models, diagnostic tests, and molecular markers. In the end, we assessed the accuracy and discriminatory power of the model using a DCA plot. We utilized the R package "ggDCA"³⁶ to generate the DCA plot for evaluating the effectiveness of the multivariable Cox model.

9 Statistics

All data processing and analysis in this article were conducted using R software (Version 4.2.2). For the comparison of two groups of continuous variables, an independent Student's t-test was used to estimate the statistical significance of normally distributed variables, while the Mann-Whitney U test (also known as the Wilcoxon rank-sum test) was used to analyze differences between non-normally distributed variables. Chi-squared test or Fisher's exact test was employed to compare and analyze the statistical significance between two groups of categorical variables.

Survival analysis was performed using the "survival" package in R. Kaplan-Meier survival curves were used to display survival differences, and the log-rank test was applied to assess the significance of differences in survival time between the two groups of patients. Both univariable and multivariable Cox analyses were conducted based on the "survival" package in R, and lasso analysis was performed using the "glmnet"³⁷ package in R.

All reported statistical p-values are two-sided, and a significance level of P < 0.05 was considered statistically significant.

1 Standardized processing of endometrial cancer datasets

The whole study flowgram is presented in Figure 1. We first performed standardization on three endometrial cancer (EC) datasets: the TCGA-UCEC dataset, the GSE6008 dataset, and the GSE17025 dataset using the R package "limma" (Figure 2A-F). The TCGA-UCEC dataset consists of 552 EC samples (grouped as EC) (Figure 2A-B), the GSE6008 dataset includes 99 EC samples (grouped as EC) (Figure 2C-D), and the GSE17025 dataset contains 91 EC samples (grouped as EC) (Figure 2E-F). From Figure 2A-F, it can be observed that the expression profiles of the three EC datasets (TCGA-UCEC, GSE6008, and GSE17025) tend to exhibit consistent data expression patterns after undergoing standardization among the samples.

2 Construction of prognostic models of GRGs and clinical correlation analysis

To determine the prognostic value of 113 GRGs (Table 2) in the TCGA-UCEC dataset, we employed LASSO regression analysis to construct a prognostic model for GRGs (Figure 3A). Additionally, we visualized the results of the LASSO regression analysis to obtain the LASSO variable trajectory plot (Figure 3B). From the plot, it can be observed that our constructed prognostic model for GRGs includes 13 genes: AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B. Furthermore, we visualized the sample grouping in the constructed prognostic model using a risk factor plot (Figure 3C). The risk factor plot consists of three parts: 1. Risk group: The samples in the dataset were grouped based on the RiskScore (risk score) predicted by the GRGs prognostic model, with the median score as the cutoff; 2. Survival outcome: The survival time and survival outcome of clinical samples in the TCGA-UCEC dataset are displayed using a dot plot; 3. Heatmap: The expression patterns of the prognostic GRGs in the GRGs prognostic model samples are visualized through a heatmap.

To analyze the differences in gene expression between high and low RiskScore groups in the TCGA-UCEC dataset based on the GRGs prognostic model, we performed differential analysis using the limma package. The analysis identified differentially expressed genes (DEGs) between the high and low RiskScore groups in the TCGA-UCEC dataset. The results are as follows: a total of 59,427 DEGs were obtained in the TCGA-UCEC dataset, among which 1,499 genes met the threshold criteria of |logFC| > 1 and P.adj < 0.05. Under this threshold, there were 688 genes with higher expression in the high RiskScore group (lower expression in the low RiskScore group, positive logFC), and 811 genes with lower expression in the high RiskScore group (higher expression in the low RiskScore group, negative logFC). We generated a volcano plot to visualize the differential analysis results of the high and low RiskScore groups in the TCGA-UCEC dataset (Figure 3D). Furthermore, we analyzed the differential expression of the 13 GRGs (AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B) between the high and low RiskScore groups in the TCGA-UCEC dataset and created a heatmap using the R package pheatmap to illustrate the specific differential expression patterns (Figure 3E). From Figure 3E, it can be observed that the expression patterns of GRGs AKR1A1, GCK, NUP155, NUP188, PFKM, and PPP2R1B exhibit significant differences between the high and low RiskScore groups in the TCGA-UCEC dataset.

3 Analysis of the expression difference of GRGs prognostic model in patients with endometrial cancer

To analyze the expression differences of 13 GRGs (AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B) between the high-risk group (group: High) and low-risk group (group: Low) of endometrial cancer (EC) patients in three endometrial cancer datasets (TCGA-UCEC dataset, GSE6008 dataset, and GSE17025 dataset), we multiplied the coefficients of the LASSO model with the gene expression levels of each sample in the GSE6008 dataset and GSE17025 dataset, and then summed them up. This process yielded risk scores (RiskScores) for each sample in the GSE6008 dataset and GSE17025 dataset. Based on the risk scores, we divided the GSE6008 dataset and GSE17025 dataset into high-risk group (group: High) and low-risk group (group: Low) using the median grouping method.

We used the Mann-Whitney U test (Wilcoxon rank sum test) to analyze the expression differences of 13 GRGs in different groups of samples from three endometrial cancer (EC) datasets. The results were presented through group comparison plots (Figure 4A, Figure 4D, Figure 4G). In the TCGA-UCEC dataset in FPKM format, the expression level of GP1 among the 13 GRGs did not show statistically significant differences among different groups (P > 0.05). However, AKR1A1 (P < 0.001), NUP155 (P < 0.001), PC (P < 0.001), PDHA1 (P < 0.001), PDHA2 (P < 0.05), and PFKM (P < 0.001) showed significantly upregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low). On the other hand, AKR1A1 (P < 0.001), ALDOB (P < 0.001), HK1 (P < 0.01), NUP188 (P < 0.001), PGK2 (P < 0.05), and PPP2R1B (P < 0.001) showed significantly downregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low) (Figure 4A).

In the GSE6008 dataset, the expression levels of AKR1A1, ALDOB, GCK, GPI, PC, and PGK2 did not show statistically significant differences among different groups (P > 0.05). However, NUP155 (P < 0.001), PDHA1 (P < 0.01), PDHA2 (P < 0.05), and PFKM (P < 0.001) showed significantly upregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low). Additionally, KK1 (P < 0.05), NUP188 (P < 0.05), and PPP2R1B (P < 0.05) showed significantly downregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low) (Figure 4D).

In the GSE17025 dataset, AKR1A1, AKR1A1, HK1, NUP188, and PC showed no statistically significant difference in expression among different groups (P > 0.05), while GPI (P < 0.05), NUP155 (P < 0.001), PDHA2 (P < 0.001), and PFKM (P < 0.01) showed significantly upregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low). ALDOB (P < 0.05), PDHA1 (P < 0.01), PGK2 (P < 0.05), and PPP2R1B (P < 0.01) showed significantly downregulated expression in the high-risk group (group: High) compared to the low-risk group (group: Low) (Figure 4G).

Then, we further plotted the ROC curves of GRGs with significant statistical significance (P<0.05) in the high-risk group (group: High) and low-risk group (group: Low) within three endometrial cancer (EC) datasets (TCGA-UCEC dataset, GSE6008 dataset, and GSE17025 dataset) to present the results (Figure 4B-C, Figure 4E-F, Figure 4H-I). The selection criteria were set as having a certain diagnostic accuracy (AUC>0.7) in at least 2 datasets. According to Figure 4, it can be observed that in the TCGA-UCEC dataset, the expression of NUP155 (AUC = 0.796, Figure 4B) has a certain diagnostic accuracy for the high-risk group (group: High) and low-risk group (group: Low) in the TCGA-UCEC dataset, and the expression of PFKM (AUC = 0.823, Figure 4C) has a certain diagnostic accuracy for the high-risk group and low-risk group in the TCGA-UCEC dataset. In the GSE6008 dataset, the expression of NUP155 (AUC = 0.793, Figure 4E) has a certain diagnostic accuracy for the high-risk group and low-risk group in the GSE6008 dataset, and the expression of PFKM (AUC = 0.733, Figure 4F) has a higher diagnostic accuracy for the high-risk group and low-risk group in the GSE6008 dataset. In the GSE17025 dataset, the expression of NUP155 (AUC = 0.742, Figure 4H) and PFKM (AUC = 0.681, Figure 4I) has lower diagnostic accuracy for the high-risk group and low-risk group in the GSE17025 dataset.

To further explore the correlation of 13 GRGs (AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B) in the TCGA-UCEC dataset, we used the Spearman algorithm to analyze the expression levels of these 13 GRGs and presented the results using a correlation heatmap (Figure 4J). The results showed that the majority of correlations between these 13 GRGs in the TCGA-UCEC dataset were statistically significant (P<0.05). We then selected the gene pairs with the highest and lowest correlations and presented the results through correlation scatter plots (Figure 4K-L). Among them, we showed the correlation analysis results of the gene pair with the highest positive correlation, PDHA1 and GPI, in the TCGA-UCEC dataset (Figure 4K, r = 0.370, P < 0.001), and the gene pair with the highest negative correlation, AKR1A1 and NUP155, in the TCGA-UCEC dataset (Figure 4L, r = -0.319, P < 0.001), using the correlation scatter plot.

4 GRGs GO and KEGG analysis

To analyze the biological processes, molecular functions, cellular components, and biological pathways related to endometrial cancer (EC) for the 13 GRGs (AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B), we first performed Gene Ontology (GO) gene function enrichment analysis on these GRGs (Table 4). The selection criteria for enriched terms were P.adj < 0.05 and FDR value (q-value) < 0.05, indicating statistical significance.

The results showed that the 13 GRGs were mainly enriched in biological processes (BP) such as the pyruvate metabolic process, monosaccharide metabolic process, hexose metabolic process, glucose metabolic process, and glycolytic process in endometrial cancer (EC). They were also enriched in cellular components (CC) such as nuclear pore, mitochondrial matrix, sperm flagellum, 9+2 motile cilium, as well as molecular functions (MF) such as monosaccharide binding, carbohydrate kinase activity, carbohydrate binding, glucose binding, structural constituent of nuclear pore. We visualized the results of the GO function enrichment analysis using a bubble plot (Figure 5A).

Additionally, we presented the BP pathways, CC pathways, and MF pathways from the results of the GO gene function enrichment analysis in the form of network diagrams (Figure 5B-D).

Next, we performed KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis on 13 GRGs (Table 5). The results showed that the 13 GRGs were significantly enriched in KEGG pathways such as Glycolysis/Gluconeogenesis, Carbon metabolism, HIF-1 signaling pathway, Central carbon metabolism in cancer, and Pyruvate metabolism. We presented the gene expression profiles of the top 5 KEGG pathways in the enrichment results using bar graphs (Figure 5E) and network diagrams (Figure 5F).

5 GSEA enrichment analysis of differentially expressed genes between high and low-risk groups

To determine the impact of gene expression levels on the development of endometrial cancer (EC) in high-risk and low-risk groups based on the GRG prognostic model in the TCGA-UCEC dataset, we conducted Gene Set Enrichment Analysis (GSEA) to analyze the association between the expression of all genes, biological processes involved, cellular components affected, and molecular functions exerted. We used a significant enrichment screening criterion of P.adj.<0.05 and FDR value (q-value) < 0.25. The results showed that genes in the high-risk and low-risk groups of the GRG prognostic model in the TCGA-UCEC dataset were significantly enriched in pathways such as il6 deprivation (Figure 6B), senescence tp53 targets (Figure 6C), apoptosis by cdknia via tp53 (Figure 6D), hypoxia (Figure 6E), cell cycle targets of tp53 and tp73 (Figure 6F) (Figure 6B-F, Table 6).

6 GSVA enrichment analysis of differentially expressed genes between high and low-risk groups

To explore the differences of hallmark gene sets between high-risk and low-risk groups in the TCGA-UCEC dataset, we performed Gene Set Variation Analysis (GSVA) on the differentially expressed genes (DEGs) between the high-risk and low-risk groups in the TCGA-UCEC dataset (Figure 7, Table 7). The GSVA analysis results of the TCGA-UCEC dataset between the high-risk and low-risk groups showed significant differences (P.adj.<0.05) in 29 hallmark gene sets, including adipogenesis, allograft rejection, apoptosis, bile acid metabolism, coagulation, and others. Among them, 12 pathways such as DNA repair, E2F targets, and G2M checkpoint obtained significantly higher enrichment scores in the high-risk group compared to the low-risk group (P.adj.<0.05), while the other 17 pathways such as adipogenesis, allograft rejection, and apoptosis showed significantly higher enrichment scores in the low-risk group compared to the high-risk group (P.adj.<0.05). Based on the GSVA results, we analyzed the differential expression of the 29 HALLMARK pathways between different groups in the TCGA-UCEC dataset and used the R package pheatmap to generate a heatmap to visualize the specific differential analysis results (Figure 7A). We used the Mann-Whitney U test to analyze the degree of group differences in the HALLMARK pathways between different groups in the TCGA-LUAD dataset and presented the results through comparative graphics (Figure 7B). The results showed that the differential expression of the 29 HALLMARK pathways between different groups in the TCGA-UCEC dataset had significant statistical differences (P<0.05).

7 Protein-protein interaction networks (PPI networks), mRNA-miRNA, mRNA-RBP, mRNA-TF, and mRNA-drugs interaction networks of GRGs

We used the STRING database to perform protein-protein interaction (PPI) analysis on 13 gene-related genes (GRGs) (AKR1A1, ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2, PPP2R1B) in the high-risk and low-risk groups of the TCGA-UCEC dataset of the GRGs prognosis model. We only kept the GRGs that had connections with other nodes and constructed a PPI network consisting of 11 GRGs (ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2). The PPI network, which includes all the nodes, serves as key genes (hub genes). We visualized the PPI network using Cytoscape software (Figure 8A).

We used the mRNA-miRNA data from the miRDB database to predict the interaction between 11 key genes (hub genes: ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) and miRNAs. We then visualized the mRNA-miRNA interaction network by drawing it using Cytoscape software (Figure 8B). In the mRNA-miRNA interaction network, the light blue ellipses represent mRNA, and the dark green circles represent miRNA. From the mRNA-miRNA interaction network, we identified that our mRNA-miRNA network consists of 7 hub genes (GPI, HK1, NUP155, NUP188, PC, PDHA1, PFKM), and 61 miRNA molecules, resulting in a total of 67 interactions between mRNA and miRNA. The specific mRNA-miRNA interaction relationships can be found in Table S1.

We also used the ENCORI database to predict the interaction between 11 key genes (hub genes: ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) and RNA binding proteins (RBPs). We then filtered out the mRNA-RBP interaction pairs with clusterNum > 4 and clipExpNum > 4 and visualized the mRNA-RBP interaction network using Cytoscape software (Figure 8C). In the mRNA-RBP interaction network, the light blue circles represent mRNA, and the purple circles represent RBPs. From the mRNA-RBP interaction network (Figure 8C), we identified that our mRNA-RBP network consists of 7 hub genes (GPI, HK1, NUP155, NUP188, PC, PDHA1, PFKM), and 24 RBP molecules, resulting in a total of 123 interactions between mRNA and RBP. Specifically, the hub gene GPI interacts with 23 RBP molecules. The specific mRNA-RBP interaction relationships can be found in Table S2.

We searched for transcription factors (TFs) that bind to the key genes (hub genes) using the CHIPBase database (version 3.0) and the hTFtarget database. After downloading the interaction data from both databases, we intersected it with the 11 key genes, resulting in a final set of 8 hub genes (ALDOB, GCK, HK1, NUP155, NUP188, PC, PDHA1, PFKM) and 139 transcription factors (TFs). We visualized this interaction data using Cytoscape software, where the light blue circles represent mRNA and the light green circles represent transcription factors (TFs) (Figure 8D). In the mRNA-TF interaction network, the DEGs BAK1 has the highest number of interactions with TFs, while the PC gene has a total of 108 pairs of mRNA-TF interactions. The specific mRNA-TF interaction relationships can be found in Table S3.

The CTD database was used to identify potential drugs or molecular compounds for the 11 key genes (hub genes). We searched the CTD database for 21 potential drugs or molecular compounds that are associated with 9 key genes (ALDOB, GCK, GPI, HK1, NUP155, PC, PDHA1, PFKM, PGK2), as shown in the mRNA-drug interaction network. The light blue oval shapes represent mRNA, while the pink circular shapes represent drugs (Figure 8D). Among them, we found 9 drugs or molecular compounds targeting the GCK gene. The specific mRNA-drug interaction relationships can be found in Table S4.

8 Construction of a multivariate Cox prognostic model for the TCGA-UCEC dataset

To determine the prognostic value of the 11 GRGs (ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) in the TCGA-UCEC dataset, we first performed univariate Cox regression analysis on the expression levels of the 11 GRGs in the TCGA-UCEC dataset. Variables with a univariate analysis P<0.1 were included in the multivariate Cox regression analysis to construct a prognostic model for the high-risk and low-risk groups of GRGs (Table 8). Then, we presented the results of the univariate Cox regression analysis in the form of a forest plot (Figure 9A) after organizing them.

Next, we conducted nomogram analysis to assess the prognostic ability of the model and plotted a column line chart (Figure 9B). The column line chart, also known as a nomogram chart, is based on a multivariate regression analysis. It represents the variables within the multivariate regression model through a set scale to score and predict the probability of events occurring. From Figure 9B, it can be observed that PGK2 among the 11 genes in the TCGA-UCEC dataset's GRGs prognostic model has significantly higher utility compared to other variables. Additionally, we performed calibration analysis for the multivariate Cox regression model using the column line chart (nomogram chart) for 1-year (Figure 9C), 3-year (Figure 9D), and 5-year (Figure 9E) prognosis predictions. We also plotted calibration curves (Figure 9C-E), where the x-axis represents the predicted survival probabilities of the model, and the y-axis represents the actual survival probabilities. The lines and points of various colors represent the model's predictions at different time points. The closer the lines of different colors are to the gray ideal line, the better the predictive performance at that time point.

We then used Decision Curve Analysis (DCA) to evaluate the clinical utility of the constructed multivariate Cox regression model at 1-year (Figure 9F), 3-year (Figure 9G), and 5-year (Figure 9H) time points and presented the results (Figure 9F-H). In the DCA chart, the x-axis represents the probability threshold or threshold probability, while the y-axis represents the net benefit. By observing the stability and higher values of the model's lines compared to the "All positive" and "All negative" lines, we can determine the range of x-values that indicates better performance. The larger the range of x-values, the better the model's effectiveness. The results show that the clinical predictive performance of our constructed multivariate Cox regression model is ranked as follows: 5-year > 3-year > 1-year.

9 Prognostic analysis of GRGs

We performed a prognostic analysis on the 11 GRGs (ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) of the multivariate Cox regression model constructed from the TCGA-UCEC dataset (Table 9). We considered molecular factors with a significance level of P<0.05 to have statistical significance.

We plotted survival KM curves for each of the 11 GRGs (ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) using a screening threshold of P<0.05. We obtained five immune-cell death-related prognostic genes that met the criteria (Figure 7A-E). These five genes are: GCK (P<0.001, Figure 10A), NUP155 (P=0.041, Figure 10B), NUP188 (P=0.022, Figure 10C), PC (P=0.026, Figure 10D), and PDHA1 (P=0.049, Figure 10E).

In addition, we also performed a time-dependent ROC curve analysis to evaluate the correlation between the 11 GRGs (GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) and the survival information (OS event) of endometrial cancer (EC) patients in the TCGA-UCEC dataset. The genes with an AUC>0.6 were selected for further analysis. The results showed that GCK (AUC1=0.665, AUC3=0.559, AUC5=0.598, Figure 10F), NUP155 (AUC1=0.617, AUC3=0.624, AUC5=0.648, Figure 10G), PC (AUC1=0.577, AUC3=0.593, AUC5=0.614, Figure 10H), and PFKM (AUC1=0.588, AUC3=0.635, AUC5=0.674, Figure 10I) had relatively low accuracy in predicting the survival outcome (OS) of EC patients in the TCGA-UCEC dataset.

The research in bioinformatics has contributed to a better understanding of complex diseases, including EC, and has led to the identification of key genes involved in EC pathogenesis. In this particular analysis, 11 genes (ALDOB, GCK, GPI, HK1, NUP155, NUP188, PC, PDHA1, PDHA2, PFKM, PGK2) were uncovered using Lasso regression and PPI network nodes as methods. Moreover, the analysis identified 61 miRNA molecules, 24 RBPs (RNA-binding proteins), 139 TFs, and 21 potential drugs or molecular compounds that might be linked to these 11 key genes. This indicates their potential roles in the development and progression of EC. Calibration analysis and Decision Curve Analysis were conducted, which revealed that the multivariate Cox regression model demonstrates its highest clinical predictive performance at 5 years, followed by 3 years and 1 year. Additionally, according to time-dependent ROC curves, six genes (GPI, HK1, NUP188, PDHA1, PDHA2, PGK2) exhibited higher predictive accuracy. Furthermore, Kaplan-Meier curves were utilized to analyze the prognosis of the identified genes. Among them, five genes (GCK, NUP155, NUP188, PC, PDHA1) showed associations with immunogenic cell death-related prognostic information in EC patients.

Our study findings are consistent with several published studies. Among the 11 GRGs examined, GPI is a predictor of outcomes in UCEC patients^36–38 as well as other cancers^39–42. Similarly, other genes such as HK1^{43, 44}, PDHA1^45–47, PGK2⁴⁸, NUP155⁴⁹, PC^50–52, ALDOB^53–56, and PFKM^{57, 58} have been implicated in carcinogenesis and have demonstrated effects on patient prognoses in various cancer types. We observed that PDHA1 and GPI exhibit the highest positive correlation. Both genes play significant roles as specific metabolic enzymes in governing the Warburg effect and glutaminolysis processes⁵⁹. Additionally, our initial findings indicate that PGK2 holds greater utility compared to other variables in the 11-GRG model, thereby highlighting its potential value in predicting EC prognosis and guiding interventions. Furthermore, we identified TFs, miRNAs, and RBPs, as well as potential drugs and small molecule compounds that interact with the key GRGs. These findings provide a foundation for future research focused on uncovering novel mechanisms and developing therapies for EC.

Nonetheless, our research had some limitations. First, this study is retrospective and exclusively focuses on the TCGA and GEO datasets, necessitating additional validation through preliminary experiments. Second, the predictive ability of the prognostic signature, consisting of 11 GRGs, and its potential role in EC intervention needs to be validated through testing with other clinical samples. Lastly, the study is hindered by limited human samples, necessitating additional clinical samples and experiments for result validation shortly.

Our study has revealed potential biomarkers and targets for the development of EC, including 11 GRGs. Among these genes, PGK2 stands out as particularly valuable within the 11-GRG model, indicating its potential clinical significance in EC. However, further research is needed to explore their specific biological functions.

Ethics Approval and Consent to Participate

All the data analyzed in the current study are published data, which is not suitable for ethics approval in this study.

Consent for Publication

All authors agree for the publication.

Availability of Data and Materials

The datasets analyzed during the current study are available publicly. The endometrial carcinoma (EC) dataset (TCGA-UCEC) from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). The corresponding clinical data for the samples were obtained from the UCSC Xena database (http://genome.ucsc.edu). The expression profile datasets of endometrial carcinoma (EC) patients from the GEO database with GSE600810-14 and GSE17025.

Competing Interests

None.

Funding

This work was supported by my own personal sources.

Authors’ contributions

Hailan Su and Hailan Su had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

Concepet and design: Hailan Su

Drafting of the manuscript: Hailan Su

Critical revision of the manuscript: Hailan Su

Statistical analysis：Hailan Su and Zhenguo Qiao

Obtained funding：Hailan Su

Administrative, technical, or material support: Hailan Su

Supervision: Hailan Su

Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. May 2021;71(3):209-249. doi:10.3322/caac.21660
Urick ME, Bell DW. Clinical actionability of molecular targets in endometrial cancer. Nat Rev Cancer. Sep 2019;19(9):510-521. doi:10.1038/s41568-019-0177-x
Weiderpass E, Antoine J, Bray FI, Oh JK, Arbyn M. Trends in corpus uteri cancer mortality in member states of the European Union. Eur J Cancer. Jun 2014;50(9):1675-84. doi:10.1016/j.ejca.2014.02.020
Cairns RA, Harris IS, Mak TW. Regulation of cancer cell metabolism. Nat Rev Cancer. Feb 2011;11(2):85-95. doi:10.1038/nrc2981
Liberti MV, Locasale JW. The Warburg Effect: How Does it Benefit Cancer Cells? Trends Biochem Sci. Mar 2016;41(3):211-218. doi:10.1016/j.tibs.2015.12.001
Colaprico A ST, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. . TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. . Nucleic Acids Res 2016 May 5;44(8):e71
Goldman MJ CB, Hastie M, Repečka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN, Zhu J, Haussler D. . Visualizing and interpreting cancer genomics data via the Xena platform. . Nat Biotechnol 2020 Jun;38(6):675-678
Barrett T TD, Wilhite SE, Ledoux P, Rudnev D, Evangelista C, Kim IF, Soboleva A, Tomashevsky M, Edgar R. . NCBI GEO: mining tens of millions of expression profiles--database and tools update. . Nucleic Acids Res 2007 Jan;35(Database issue):D760-5.
Davis S MP. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. . Bioinformatics 2007 Jul 15;23(14):1846-7.
Bommer GT FY, Iura A, Giordano TJ, Kuick R, Kadikoy H, Sikorski D, Wu R, Cho KR, Fearon ER. . IRS1 regulation by Wnt/beta-catenin signaling and varied contribution of IRS1 to the neoplastic phenotype. . J Biol Chem 2010 Jan 15;285(3):1928-38
Hendrix ND WR, Kuick R, Schwartz DR, Fearon ER, Cho KR. . Fibroblast growth factor 9 has oncogenic activity and is a downstream target of Wnt signaling in ovarian endometrioid adenocarcinomas. . Cancer Res 2006 Feb 1;66(3):1354-62.
Liu Y KR, Hanash S, Richardson B. . DNA methylation inhibition increases T cell KIR expression through effects on both promoter methylation and transcription factors. . Clin Immunol 2009 Feb;130(2):213-24.
Wu R H-LN, Kuick R, Zhai Y, Schwartz DR, Akyol A, Hanash S, Misek DE, Katabuchi H, Williams BO, Fearon ER, Cho KR. . Mouse model of human ovarian endometrioid adenocarcinoma based on somatic defects in the Wnt/beta-catenin and PI3K/Pten signaling pathways. . Cancer Cell 2007 Apr;11(4):321-33
Wu R ZY, Kuick R, Karnezis AN, Garcia P, Naseem A, Hu TC, Fearon ER, Cho KR. . Impact of oviductal versus ovarian epithelial cell of origin on ovarian endometrioid carcinoma phenotype in the mouse. . J Pathol 2016 Nov;240(3):341-351
Day RS MK, Chandran UR, Lisovich A, Conrads TP, Hood BL, Kolli VS, Kirchner D, Litzi T, Maxwell GL. . Identifier mapping performance for integrating transcriptomics and proteomics experimental results. . BMC Bioinformatics 2011 May 27;12:213
Day RS MK. A decision theory paradigm for evaluating identifier mapping and filtering methods using data integration. . BMC Bioinformatics 2013 Jul 15;14:223
Stelzer G RN, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y, Kaplan S, Dahary D, Warshawsky D, Guan-Golan Y, Kohn A, Rappaport N, Safran M, Lancet D. . The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. . Curr Protoc Bioinformatics 2016 Jun 20;54:1301-13033.
Bi J BF, Pan X, Yang Q. . Establishment of a novel glycolysis-related prognostic gene signature for ovarian cancer and its relationships with immune infiltration of the tumor microenvironment. . J Transl Med 2021 Sep 8;19(1):382
Zhang D ZY, Yang S, Li Y, Wang M, Yao J, Deng Y, Li N, Wei B, Wu Y, Zhu Y, Li H, Dai Z. . Identification of a Novel Glycolysis-Related Gene Signature for Predicting Breast Cancer Survival. . Front Oncol 2021 Jan 8;10:596087.
Zheng J GJ, Zhu L, Zhou Y, Tong J. . Comprehensive analyses of glycolysis-related lncRNAs for ovarian cancer patients. . J Ovarian Res 2021 Sep 24;14(1):124
Liberzon A BC, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. . The Molecular Signatures Database (MSigDB) hallmark gene set collection. . Cell Syst 2015 Dec 23;1(6):417-425.
G. Y. Gene Ontology Semantic Similarity Analysis Using GOSemSim. . Methods Mol Biol 2020;2117:207-215
Kanehisa M GS. KEGG: kyoto encyclopedia of genes and genomes. . Nucleic Acids Res 2000 Jan 1;28(1):27-30.
Subramanian A TP, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. . Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. . Proc Natl Acad Sci U S A 2005 Oct 25;102(43):15545-50.
Wang HQ YX, Liu ZH, Cheng X, Samartzis D, Jia LT, Wu SX, Huang J, Chen J, Luo ZJ. . Deregulated miR-155 promotes Fas-mediated apoptosis in human intervertebral disc degeneration by targeting FADD and caspase-3. . J Pathol 2011 Oct;225(2):232-42.
Szklarczyk D GA, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. . STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. . Nucleic Acids Res 2019 Jan 8;47(D1):D607-D613.
Chen Y WX. miRDB: an online database for prediction of functional microRNA targets. . Nucleic Acids Res 2020 Jan 8;48(D1):D127-D131.
Li JH LS, Zhou H, Qu LH, Yang JH. . starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. . Nucleic Acids Res 2014 Jan;42(Database issue):D92-7.
Zhou KR LS, Sun WJ, Zheng LL, Zhou H, Yang JH, Qu LH. . ChIPBase v2.0: decoding transcriptional regulatory networks of non-coding RNAs and protein-coding genes from ChIP-seq data. . Nucleic Acids Res 2017 Jan 4;45(D1):D43-D50
Zhang Q LW, Zhang HM, Xie GY, Miao YR, Xia M, Guo AY. . hTFtarget: A Comprehensive Database for Regulations of Human Transcription Factors and Their Targets. . Genomics Proteomics Bioinformatics 2020 Apr;18(2):120-128
Davis AP GC, Johnson RJ, Sciaky D, Wiegers J, Wiegers TC, Mattingly CJ. . Comparative Toxicogenomics Database (CTD): update 2021. . Nucleic Acids Res 2021 Jan 8;49(D1):D1138-D1143
SY. P. Nomogram: An analogue tool to deliver digital knowledge. . J Thorac Cardiovasc Surg 2018 Apr;155(4):1793.
Perkins NJ WJ, Mumford SL, Sjaarda LA, Mitchell EM, Pollack AZ, Schisterman EF. . Combining Biomarker Calibration Data to Reduce Measurement Error. . Epidemiology 2019 Nov;30 Suppl 2(Suppl 2):S3-S9
Van Calster B WL, Verbeek JFM, Verbakel JY, Christodoulou E, Vickers AJ, Roobol MJ, Steyerberg EW. . Reporting and Interpreting Decision Curve Analysis: A Guide for Investigators. . Eur Urol 2018 Dec;74(6):796-804.
Engebretsen S BJ. Statistical predictions with glmnet. . Clin Epigenetics 2019 Aug 23;11(1):123.
Chen J, Wang G, Zhang D, Luo X, Zhang D, Zhang Y. Construction of novel hypoxia-related gene model for prognosis and tumor microenvironment in endometrial carcinoma. Front Endocrinol (Lausanne). 2022;13:1075431. doi:10.3389/fendo.2022.1075431
Cheng Y, Li X, Dai Y, Dong Y, Yang X, Wang J. Identification of an immune-related risk signature and nomogram predicting the overall survival in patients with endometrial cancer. J Gynecol Oncol. May 2021;32(3):e30. doi:10.3802/jgo.2021.32.e30
Fan Y, Li X, Tian L, Wang J. Identification of a Metabolism-Related Signature for the Prediction of Survival in Endometrial Cancer Patients. Front Oncol. 2021;11:630905. doi:10.3389/fonc.2021.630905
Djafarzadeh R, Noessner E, Engelmann H, et al. GPI-anchored TIMP-1 treatment renders renal cell carcinoma sensitive to FAS-meditated killing. Oncogene. Mar 9 2006;25(10):1496-508. doi:10.1038/sj.onc.1209188
Jiang WW, Zahurak M, Zhou ZT, et al. Alterations of GPI transamidase subunits in head and neck squamous carcinoma. Mol Cancer. Nov 21 2007;6:74. doi:10.1186/1476-4598-6-74
Naghibalhossaini F, Ebadi P. Evidence for CEA release from human colon cancer cells by an endogenous GPI-PLD enzyme. Cancer Lett. Mar 28 2006;234(2):158-67. doi:10.1016/j.canlet.2005.03.028
Wu G, Guo Z, Chatterjee A, et al. Overexpression of glycosylphosphatidylinositol (GPI) transamidase subunits phosphatidylinositol glycan class T and/or GPI anchor attachment 1 induces tumorigenesis and contributes to invasion in human breast cancer. Cancer Res. Oct 15 2006;66(20):9829-36. doi:10.1158/0008-5472.Can-06-0506
Chen QT, Zhang ZY, Huang QL, et al. HK1 from hepatic stellate cell-derived extracellular vesicles promotes progression of hepatocellular carcinoma. Nat Metab. Oct 2022;4(10):1306-1321. doi:10.1038/s42255-022-00642-5
Jin J, Guo D, Wang Y, Jiao W, Li D, He Y. Artesunate Inhibits the Development of Esophageal Cancer by Targeting HK1 to Reduce Glycolysis Levels in Areas With Zinc Deficiency. Front Oncol. 2022;12:871483. doi:10.3389/fonc.2022.871483
Deng L, Jiang A, Zeng H, Peng X, Song L. Comprehensive analyses of PDHA1 that serves as a predictive biomarker for immunotherapy response in cancer. Front Pharmacol. 2022;13:947372. doi:10.3389/fphar.2022.947372
Liu Z, Yu M, Fei B, Fang X, Ma T, Wang D. miR‑21‑5p targets PDHA1 to regulate glycolysis and cancer progression in gastric cancer. Oncol Rep. Nov 2018;40(5):2955-2963. doi:10.3892/or.2018.6695
Zhong C, Long Z, Yang T, et al. M6A-modified circRBM33 promotes prostate cancer progression via PDHA1-mediated mitochondrial respiration regulation and presents a potential target for ARSI therapy. Int J Biol Sci. 2023;19(5):1543-1563. doi:10.7150/ijbs.77133
Wu ST, Liu B, Ai ZZ, et al. Esculetin Inhibits Cancer Cell Glycolysis by Binding Tumor PGK2, GPD2, and GPI. Front Pharmacol. 2020;11:379. doi:10.3389/fphar.2020.00379
Holzer K, Ori A, Cooke A, et al. Nucleoporin Nup155 is part of the p53 network in liver cancer. Nat Commun. May 14 2019;10(1):2147. doi:10.1038/s41467-019-10133-z
Lin Q, He Y, Wang X, et al. Targeting Pyruvate Carboxylase by a Small Molecule Suppresses Breast Cancer Progression. Adv Sci (Weinh). May 2020;7(9):1903483. doi:10.1002/advs.201903483
Sellers K, Fox MP, Bousamra M, 2nd, et al. Pyruvate carboxylase is critical for non-small-cell lung cancer proliferation. J Clin Invest. Feb 2015;125(2):687-98. doi:10.1172/jci72873
Shu Y, Yang N, Cheng N, et al. Intervening pyruvate carboxylase stunts tumor growth by strengthening anti-tumor actions of tumor-associated macrophages. Signal Transduct Target Ther. Feb 2 2022;7(1):34. doi:10.1038/s41392-021-00807-w
He J, Jin Y, Chen Y, et al. Downregulation of ALDOB is associated with poor prognosis of patients with gastric cancer. Onco Targets Ther. 2016;9:6099-6109. doi:10.2147/ott.S110203
Jia W, Wu Q, Yu X, et al. Prognostic values of ALDOB expression and (18)F-FDG PET/CT in hepatocellular carcinoma. Front Oncol. 2022;12:1044902. doi:10.3389/fonc.2022.1044902
Shao Y, Wu B, Yang Z, et al. ALDOB represents a potential prognostic biomarker for patients with clear cell renal cell carcinoma. Transl Androl Urol. Apr 28 2023;12(4):549-571. doi:10.21037/tau-22-743
Xu JX, Qin SL, Wei HW, Chen YY, Peng YC, Qi LN. Down-regulation of ALDOB during metabolic reprogramming mediates malignant behavior in hepatocellular carcinoma and insensitivity to postoperative adjuvant transarterial chemoembolization. Clin Sci (Lond). Feb 27 2023;137(4):303-316. doi:10.1042/cs20220661
Ahsan H, Halpern J, Kibriya MG, et al. A genome-wide association study of early-onset breast cancer identifies PFKM as a novel breast cancer gene and supports a common genetic spectrum for breast cancer at any age. Cancer Epidemiol Biomarkers Prev. Apr 2014;23(4):658-69. doi:10.1158/1055-9965.Epi-13-0340
Gao W, Huang M, Chen X, et al. The role of S-nitrosylation of PFKM in regulation of glycolysis in ovarian cancer cells. Cell Death Dis. Apr 15 2021;12(4):408. doi:10.1038/s41419-021-03681-0
Huang X, Gan G, Wang X, Xu T, Xie W. The HGF-MET axis coordinates liver cancer metabolism and autophagy for chemotherapeutic resistance. Autophagy. Jul 2019;15(7):1258-1279. doi:10.1080/15548627.2019.1580105

Table 1. List of Uterine Corpus Endometrioid Carcinoma Information.

	TCGA-UCEC	GSE6008	GSE17025
Platform	TCGA	GPL96	GPL570
Species	Homo sapiens	Homo sapiens	Homo sapiens
Tissue	Endometrial tissues	Endometrial tissues	Endometrial tissues
Samples in EC group	552	99	91
Reference	/	Fibroblast growth factor 9 has oncogenic activity and is a downstream target of Wnt signaling in ovarian endometrioid adenocarcinomas IRS1 regulation by Wnt/beta-catenin signaling and varied contribution of IRS1 to the neoplastic phenotype Mouse model of human ovarian endometrioid adenocarcinoma based on somatic defects in the Wnt/beta-catenin and PI3K/Pten signaling pathways Impact of oviductal versus ovarian epithelial cell of origin on ovarian endometrioid carcinoma phenotype in the mouse DNA methylation inhibition increases T cell KIR expression through effects on both promoter methylation and transcription factors	Identifier mapping performance for integrating transcriptomics and proteomics experimental results A decision theory paradigm for evaluating identifier mapping and filtering methods using data integration

TCGA：The cancer genome atlas；UCEC：Uterine corpus endometrioid carcinoma。

Table 2. Glycolysis related genes related genes list.

Glycolysis related genes
AAAS	ADPGK	AKR1A1	ALDOA	ALDOB	ALDOC
BPGM	CD44	DLAT	DLD	ECD	EGFR
EGLN3	ENO1	ENO2	ENO3	FBP1	FBP2
G6PC1	G6PD	GAPDH	GAPDHS	GCK	GCKR
GNPDA1	GNPDA2	GOT1	GOT2	GPI	GYS1
HK1	HK2	HK3	HSPA5	IDH1	IRS2
KIF2A	LDHA	LDHAL6B	LDHB	LDHC	LHPP
MDH1	MDH2	MET	NDC1	NUP107	NUP133
NUP153	NUP155	NUP160	NUP188	NUP205	NUP210
NUP214	NUP35	NUP37	NUP42	NUP43	NUP50
NUP54	NUP58	NUP62	NUP85	NUP88	NUP93
NUP98	PC	PCK1	PCK2	PDHA1	PDHA2
PDHB	PDK3	PFKFB1	PFKFB2	PFKFB3	PFKFB4
PFKL	PFKM	PFKP	PGAM1	PGAM2	PGK1
PGK2	PGM2L1	PGP	PKLR	PKM	PLOD2
POM121	PPIA	PPP2CA	PPP2CB	PPP2R1A	PPP2R1B
PPP2R5D	PRKACA	PRKACB	PRKACG	RAE1	RANBP2
SEC13	SEH1L	SLC16A3	SLC25A13	SLC37A4	SOD1
STC2	TKTL1	TPI1	TPR	VEGFA

Table 3. Differntial expression genes of glycolysis between high and low risk group.

GeneSymbol	logFC	AveExpr	t	P.Value	adj.P.Val	B
AKR1A1	-0.44436	12.27991	-10.4157	2.61E-23	3.33E-21	42.16041
ALDOB	-0.49031	2.343936	-3.69638	0.000241	0.001183	-0.28922
GCK	1.410065	3.397963	8.780242	2.06E-17	9.83E-16	28.82982
GPI	0.062912	13.72358	1.125659	0.2608	0.402255	-6.36443
HK1	-0.20793	11.94134	-4.19413	3.19E-05	0.00019	1.607566
NUP155	0.714279	10.06487	14.43725	2.66E-40	1.22E-36	80.69049
NUP188	-0.32278	11.33827	-6.82173	2.38E-11	4.29E-10	15.1919
PC	0.216571	8.984535	2.883538	0.004086	0.014626	-2.89376
PDHA1	0.256345	11.1641	4.7602	2.48E-06	1.85E-05	4.036105
PDHA2	0.049934	0.143489	1.414204	0.157867	0.27618	-6.00145
PFKM	1.008421	9.910202	14.53331	9.67E-41	5.22E-37	81.68844
PGK2	-0.29014	0.719913	-3.46286	0.000576	0.002584	-1.10064
PPP2R1B	-0.46775	10.01081	-9.90734	2.10E-21	1.97E-19	37.84582

Table 4. GO enrichment analysis results of GRGs.

Ontology	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust
BP	GO:0006090	pyruvate metabolic process	9/13	106/18800	2.86e-18	9.88e-16
BP	GO:0005996	monosaccharide metabolic process	10/13	261/18800	6.17e-17	1.07e-14
BP	GO:0019318	hexose metabolic process	9/13	242/18800	5.72e-15	6.59e-13
BP	GO:0006006	glucose metabolic process	8/13	201/18800	1.83e-13	1.58e-11
BP	GO:0006096	glycolytic process	6/13	81/18800	8.87e-12	5.52e-10
CC	GO:0005643	nuclear pore	2/13	93/19594	0.0017	0.0318
CC	GO:0005759	mitochondrial matrix	3/13	473/19594	0.0033	0.0318
CC	GO:0036126	sperm flagellum	2/13	136/19594	0.0035	0.0318
CC	GO:0097729	9+2 motile cilium	2/13	149/19594	0.0042	0.0318
MF	GO:0048029	monosaccharide binding	5/13	71/18410	9.29e-10	3.53e-08
MF	GO:0019200	carbohydrate kinase activity	3/13	21/18410	3.63e-07	6.9e-06
MF	GO:0030246	carbohydrate binding	5/13	270/18410	7.64e-07	9.68e-06
MF	GO:0005536	glucose binding	2/13	11/18410	2.52e-05	0.0002
MF	GO:0017056	structural constituent of nuclear pore	2/13	28/18410	0.0002	0.0013

GO：Gene ontology；GRGs：Glycolysis related genes； BP：Biological process；CC：Cellular component；MF：Molecular function；。

Table 5. KEGG enrichment analysis results of GRGs.

Ontology	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust
KEGG	hsa00010	Glycolysis / Gluconeogenesis	9/13	67/8164	6.74e-17	3.03e-15
KEGG	hsa01200	Carbon metabolism	9/13	115/8164	1.09e-14	2.44e-13
KEGG	hsa04066	HIF-1 signaling pathway	6/13	109/8164	7.85e-09	1.18e-07
KEGG	hsa05230	Central carbon metabolism in cancer	5/13	70/8164	4.89e-08	5.5e-07
KEGG	hsa00620	Pyruvate metabolism	4/13	47/8164	6.64e-07	5.97e-06

KEGG：Kyoto Encyclopedia of Genes and Genomes；GRGs：Glycolysis related genes。

Table 6. GSEA enrichment analysis results of GRGs Prognosis Model high and low risk group genes.

Description	setSize	enrichmentScore	NES	p.adjust	qvalues
CROONQUIST_IL6_DEPRIVATION_DN	95	0.759988	2.367029	0.002049	0.023994
TANG_SENESCENCE_TP53_TARGETS_DN	56	0.7642	2.181746	0.002141	0.023994
WU_APOPTOSIS_BY_CDKN1A_VIA_TP53	53	0.76208	2.156821	0.00211	0.023994
MANALO_HYPOXIA_DN	291	0.554778	1.94721	0.002227	0.023994
SCIAN_CELL_CYCLE_TARGETS_OF_TP53_AND_TP73_DN	22	0.781271	1.84911	0.002033	0.023994
REACTOME_REGULATION_OF_GLUCOKINASE_BY_GLUCOKINASE_REGULATORY_PROTEIN	32	0.649762	1.674905	0.004065	0.036881
CHUANG_OXIDATIVE_STRESS_RESPONSE_DN	11	0.818894	1.659336	0.003992	0.036881
REACTOME_REGULATION_OF_TP53_ACTIVITY_THROUGH_PHOSPHORYLATION	92	0.515639	1.597948	0.002045	0.023994
DEBIASI_APOPTOSIS_BY_REOVIRUS_INFECTION_UP	300	0.447573	1.578698	0.002183	0.023994
CONCANNON_APOPTOSIS_BY_EPOXOMICIN_DN	174	0.419929	1.386314	0.004505	0.039296

GSEA：Gene set enrichment analysis；GRGs：Glycolysis related genes。

Table 7. GSVA enrichment analysis results of GRGs Prognosis Model high and low risk group genes.

Description	logFC	AveExpr	t	adj.P.Val
HALLMARK_G2M_CHECKPOINT	-0.34003	-0.01344	-10.2356	9.76E-23
HALLMARK_E2F_TARGETS	-0.3618	-0.00498	-10.1585	1.92E-22
HALLMARK_XENOBIOTIC_METABOLISM	0.210475	-0.06142	9.844826	2.87E-21
HALLMARK_MYC_TARGETS_V1	-0.2683	-0.01393	-8.21011	1.41E-15
HALLMARK_MYC_TARGETS_V2	-0.27032	-0.04134	-8.13241	2.51E-15
HALLMARK_SPERMATOGENESIS	-0.16527	-0.0413	-8.11092	2.94E-15
HALLMARK_UV_RESPONSE_UP	-0.17053	-0.06774	-8.06417	4.16E-15
HALLMARK_MITOTIC_SPINDLE	-0.20483	-0.0493	-7.48147	2.70E-13
HALLMARK_ESTROGEN_RESPONSE_EARLY	0.180062	-0.06063	7.413842	4.31E-13
HALLMARK_ESTROGEN_RESPONSE_LATE	0.148396	-0.07473	6.637184	7.28E-11
HALLMARK_P53_PATHWAY	0.137985	-0.07991	6.165324	1.31E-09
HALLMARK_PEROXISOME	0.140932	-0.06688	6.161537	1.34E-09
HALLMARK_BILE_ACID_METABOLISM	0.135999	-0.06445	6.105933	1.86E-09
HALLMARK_FATTY_ACID_METABOLISM	0.134003	-0.04991	5.548473	4.36E-08
HALLMARK_KRAS_SIGNALING_DN	-0.09892	-0.06617	-4.68067	3.55E-06
HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY	0.114317	-0.0692	4.432693	1.11E-05
HALLMARK_DNA_REPAIR	-0.09133	-0.07313	-3.77311	0.000178
HALLMARK_MTORC1_SIGNALING	-0.10357	-0.0326	-3.56447	0.000394
HALLMARK_IL2_STAT5_SIGNALING	0.085204	-0.05111	3.309744	0.000991
HALLMARK_UNFOLDED_PROTEIN_RESPONSE	-0.09026	-0.03906	-3.2785	0.001105
HALLMARK_INTERFERON_ALPHA_RESPONSE	-0.11776	-0.0223	-3.05395	0.002361
HALLMARK_WNT_BETA_CATENIN_SIGNALING	0.087447	-0.04918	2.985886	0.002945
HALLMARK_COAGULATION	0.076256	-0.04605	2.964792	0.003152
HALLMARK_GLYCOLYSIS	0.068585	-0.0566	2.85593	0.004443
HALLMARK_APOPTOSIS	0.066269	-0.05902	2.652662	0.008202
HALLMARK_ALLOGRAFT_REJECTION	0.085361	-0.02689	2.607811	0.009344
HALLMARK_INFLAMMATORY_RESPONSE	0.075519	-0.04237	2.508724	0.012385
HALLMARK_ADIPOGENESIS	0.05611	-0.07552	2.447098	0.014693
HALLMARK_TNFA_SIGNALING_VIA_NFKB	0.070009	-0.05535	2.294833	0.022095

GSVA：Gene set variation analysis；GRGs：Glycolysis related genes。

Table 8. Patient Characteristics of UCEC patients in the TCGA datasets.

Characteristic	levels	Overall
n		552
Clinical stage, n (%)	Stage I	342 (62%)
	Stage II	51 (9.2%)
	Stage III	130 (23.6%)
	Stage IV	29 (5.3%)
Race, n (%)	Asian	20 (3.9%)
	Black or African American	108 (21.3%)
	White	379 (74.8%)
Age, n (%)	<=60	206 (37.5%)
	>60	343 (62.5%)
Histological type, n (%)	Endometrioid	410 (74.3%)
	Mixed	24 (4.3%)
	Serous	118 (21.4%)
OS event, n (%)	Alive	458 (83%)
	Dead	94 (17%)
DSS event, n (%)	Alive	487 (88.5%)
	Dead	63 (11.5%)
PFI event, n (%)	Alive	423 (76.6%)
	Dead	129 (23.4%)

TCGA：The cancer genome atlas；UCEC：Uterine corpus endometrioid carcinoma；OS：Overall survival；DSS：Disease specific survival；PFI：Progression free interval。

Table 9. COX regression to hub genes associated with OS in TCGA-UCEC.

Characteristics	Total(N)	Univariate analysis		Multivariate analysis
Characteristics	Total(N)	Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
ALDOB	552	1.714 (1.034 - 2.841)	0.037	1.577 (0.932 - 2.667)	0.089
GCK	552	1.527 (1.128 - 2.068)	0.006	1.219 (0.842 - 1.765)	0.295
GPI	552	1.356 (0.991 - 1.856)	0.057	1.267 (0.907 - 1.769)	0.165
HK1	552	0.711 (0.504 - 1.002)	0.051	0.576 (0.394 - 0.842)	0.004
NUP155	552	1.726 (1.311 - 2.272)	< 0.001	1.549 (1.100 - 2.181)	0.012
NUP188	552	0.636 (0.442 - 0.917)	0.015	0.526 (0.359 - 0.770)	< 0.001
PC	552	1.420 (1.039 - 1.939)	0.028	1.464 (1.006 - 2.131)	0.047
PDHA1	552	1.478 (1.087 - 2.010)	0.013	1.060 (0.743 - 1.510)	0.749
PDHA2	552	15.840 (1.789 - 140.274)	0.013	2.232 (0.147 - 33.863)	0.563
PFKM	552	1.557 (1.218 - 1.989)	< 0.001	1.277 (0.968 - 1.684)	0.084
PGK2	552	0.014 (0.000 - 0.816)	0.040	0.010 (0.000 - 0.916)	0.046

OS：Overall survival；TCGA：The cancer genome atlas；UCEC：Uterine corpus endometrioid carcinoma。

No competing interests reported.

SupplementaryMaterial.docx

An 11-GRG Predicting Model for UCEC Prognosis

Status:

Version 1

Abstract

Figures

Introduction

Method