From the TCGA database (https://portal.gdc.cancer.gov/), we downloaded the raw mRNA expression data of CESC, including the normal group (n = 3) and the CESC group (n = 306). To analyze differentially expressed genes, FPKM data of mRNA level 3 were integrated and normalized using the Limma package. The differential gene screening conditions were $|$LogFC$|$ $>$ 1 and p value 0.05. Next, the Series Matrix File data of GSE44001 was downloaded from GEO public database, the annotation platform is GPL14951. The Series Matrix File data of GSE52903 were annotated with GPL6244. GeneCards (https://www.genecards.org/) database was used to extract 18,762 metabolism-related genes, and genes with Relevance scores $>$3 were analyzed as metabolic gene sets.

ClusterProfiler (R3.6, R Core Team, Vienna, Austria) was used to annotate the functions of differential genes and explore their functional relevance. Relevant functional categories were identified using GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes). We considered a statistical difference if p value and q value all less than 0.05.

The prognostic correlation models were constructed using Lasso regression once differential metabolism genes were selected. For each patient, a risk score formula was constructed by incorporating the expression values of each specific gene and weighting them with their regression coefficients in Lasso regression analysis. The risk score formula is shown below.

$\text{Risk score}={\displaystyle \sum _i}\text{regression coefficient of metabolic}$

$\text{gene}i\times \text{expression of metabolic gene}i\text{.}$

Where Metabolic gene i is the identifier of the i-th selected core metabolic gene. According to this formula, median risk score was used to distinguish low-risk and high-risk patients. A Kaplan-Meier survival curve was computed and then compared using log-rank statistics. Using Lasso regression and hierarchical multiple regression, we examined the role of risk score in predicting patient prognosis. We evaluated the model prediction performance using ROC curves.

We predicted each tumor sample’s chemotherapy susceptibility using the “pRRophetic” function in R software (R3.6, R Core Team, Vienna, Austria) based on the GDSC Cancer Drug Susceptibility Genomics Database. Based on the GDSC dataset, 10 cross-validations were conducted to determine the accuracy of the IC50 estimates for each specific chemotherapy drug treatment. We set all the parameter to default settings, including the “combat”, which removes batch effects, and the average of duplicate gene expression.

RNA-seq data from different subgroups of CESC patients was analyzed using the CIBERSORT algorithm to determine the relative proportions of immune infiltrating cells. Pearson correlation analysis is performed for gene expression and immune cell content, and statistical significance is defined as p 0.05.

By thoroughly scoring the gene sets of interest, Gene set variation analysis (GSVA) turns changes at the gene level into changes at the pathway level and establishes the biological function of the samples. In this study, molecular signatures database (v7.0) gene sets were downloaded, and GSVA scores were used to assess possible biological functional changes between samples.

R software (R3.6, R Core Team, Vienna, Austria) and SPSS 19.0 (SPSS Inc., Chicago, IL, USA) were used for all statistical analyses. Kaplan-Meier survival curves were generated and compared using log-rank. Multivariate analysis is performed using Cox proportional risk models. All statistical tests were two-sided, and we consider is as statistically significant if p 0.05.

We extracted a collection of 2235 metabolism-related regulatory factors from the downloaded mRNA expression data (FPKM). The differential expression between cervical cancer patients and control patients are investigated using the Wilcox test. According to the results, 378 metabolism-related genes, including 177 up-regulated and 201 down-regulated genes, were differentially expressed (logFC $>$1, logFC –1, and p 0.05). The changes in differential gene expression values (log2(expression value+1)) are shown by heat map in Fig. 1A, while the differential ploidy and p value of these differential genes are presented using a volcano plot in Fig. 1B.

Fig. 1.

The identification of metabolism-related genes with differential expression. (A) Heat map of the differentially expressed genes (DEGs). (B) Volcano plot of the DEGs.

As a result of GO and KEGG pathway enrichment analysis, a large number of pathways were significantly enriched for differentially expressed genes, such as coenzyme metabolic process, cytoplasmic vesicle lumen, coenzyme binding, etc. (Fig. 2A). In the KEGG enrichment process, there are central carbon metabolism in cancer, biosynthesis of amino acids and other metabolism-related pathways as illustrated in Fig. 2B. At the same time, through the Metascape database, we further analyzed candidate gene pathways. According to the results, these candidates were mainly enriched in hormone, organic hydroxy compound metabolic, and small molecule biosynthetic pathways (Fig. 2C). Fig. 2D shows the co-expression network analysis result of genes in the differential gene set.

Fig. 2.

Functional analysis of metabolism-related genes. (A) DEGs enriched by GO. (B) KEGG enrichment analysis of DEGs. (C) Related pathways of DEGs. (D) Protein network analysis of DEGs.

In order to identify the key metabolic genes, we collected clinical information from CESC patients. According to Fig. 3A–C, we used Lasso regression and Cox univariate regression algorithms to identify cervical cancer signature genes. By using Cox univariate regression, 59 prognosis-related genes were identified: ISCU, TCN2, FOXP3, PGK1, UCP2, SPINT2, FASN, MOCS1, TFRC, ADH1B, TNF, MIR200A, GAPDH, MMP1, PFKFB4, GCH1, CDIPT, ACOT4, LDHA, LEPR, BCL2, SPP1, TNFRSF1B, TGFA, CPE, CA9, NPL, SPTBN1, GALNT3, DHCR7, FGFR2, PFKM, NME4, SDS, SQLE, SLC25A42, JUN, NT5E, TRPV4, HSPG2, PLA2G7, SLC25A10, CKB, APOC1, GART, SLC2A3, SCD, VDAC1, ENO1, APOD, SELP, HMGCS1, NR1D2, TREM2, PLIN2, FGFR3, PIP4K2B, COX5A and HK2. After randomly dividing TCGA patients into training and validation sets in an equal ratio, Lasso regression analysis was used to determine the best risk score. The metabolic gene prognostic risk score model was constructed as the formula below: Risk Score = FGFR2$\times$ (–0.4067) + FOXP3$\times$ (–0.2731) + NPL$\times$ (–0.2704) + CKB$\times$ (–0.1992) + APOD$\times$ (–0.1027) + TCN2$\times$ (–0.0673) + TNFRSF1B$\times$ (–0.0378) + FGFR3$\times$ (–0.0154) + HSPG2$\times$ 0.0067 + TFRC$\times$ 0.0089 + NT5E$\times$ 0.0220 + CA9$\times$ 0.0313 + CPE$\times$ 0.0438 + FASN$\times$ 0.0442 + NME4$\times$ 0.0828 + PIP4K2B$\times$ 0.0977 + SPP1$\times$ 0.0989 + MMP1$\times$ 0.1013 + JUN$\times$ 0.1311 + SPINT2$\times$ 0.1474 + GAPDH $\mathrm{\times }$ 0.1894 + MIR200A$\times$ 0.1902 + TNF$\times$ 0.1972 + SQLE$\times$ 0.2018 + ADH1B$\times$ 0.3305. Based on the median risk score, patients were divided into high and low risk groups. In both the training and test sets, the OS of the high-risk group was significantly lower (p 0.001) than that of the low-risk group, as shown in Fig. 3D,E. Furthermore, the ROC curve results showed a C-index of 0.85 and 0.79 in the training set and test set, both indicating good validation efficacy, as shown in Fig. 4A,B. Baseline characteristics of the patients is presented in Supplementary Table 1.

Fig. 3.

Prognostic risk model construction. (A,B) Lasso regression. (C) The regression coefficient of differentially expressed genes. (D,E) Training and test set survival analysis.

Fig. 4.

Validation of the risk model. (A) ROC curves of the training cohort (1-year AUC = 0.88, 3-year AUC = 0.88, 5-year AUC = 0.89). (B) ROC curves of the testing cohort (1-year AUC = 0.83, 3-year AUC = 0.86, 5-year AUC = 0.74).

The relationship between the tumor immune infiltration and the risk score was further investigated. The results showed that risk score was significantly and positively correlated with Macrophages M0 (p 0.01), Dendritic cells activated (p 0.05), Mast cells activated (p 0.01) content, and significantly negatively correlated with Mast cells resting (p 0.01), T cells CD8 (p 0.01), B cells naive (p 0.01), as shown in Fig. 5A. The interaction between immune cells is shown in Fig. 5D. A combination of surgery and chemotherapy is effective in treating cervical cancer at an early stage. We use the R package’s “pRRophetic” function to determine chemotherapy sensitivity of tumor samples from the GDSC database. A significant relationship was found between the risk score and the sensitivity to Bortezomib, Docetaxel, Erlotinib, Metformin, Mitomycin.C, and Paclitaxel (Fig. 5B). Additionally, we examined the mutation profiles of high-risk and low-risk patients. High-risk individuals had a significantly higher proportion of mutations of KMT2C and TTN than low-risk individuals, shown in Fig. 5C. At the same time, we also found that there were significant differences in microsatellite instability (MSI) and tumor mutation burden (TMB) between the two groups, and the expression of TMB and MSI was increased in the high-risk group (Fig. 5E,F).

Fig. 5.

Multiple omics maps of the risk model. (A,D) Immune infiltration and risk score. (B) Relationship between risk score and drug sensitivity. (C) Relationship between model and SNP mutation. (E,F) TMB and MSI were significantly different between the two groups.

We then studied two groups’ specific signaling pathways. Results of the GSVA revealed that ADIPOGENESIS, TNFA_SIGNALING_VIA_NFKB, ANDROGEN_RESPONSE, MYOGENESIS, and KRAS_SIGNALING_DN were the most enriched pathways for the two groups of patients, as shown in Fig. 6. It has been shown that altering these signaling pathways can affect the prognosis of cervical cancer patients.

Fig. 6.

Signaling mechanisms for prognostic model.

Data on processed CESC patients with survival information was downloaded from GEO (GSE44001, GSE52903). And Kaplan-Meier survival analysis was used to assess survival differences. Compared to the low-risk group, OS was significantly lower in the high-risk group (Fig. 7A,B). ROC curve analysis was performed on external data sets to verify the model’s accuracy. The results showed that the model has a strong predictive effect on the prognosis of patients (GSE44001-C-index = 0.68, and GSE52903-C-index = 0.73), as shown in Fig. 7C,D.

Fig. 7.

Robustness of the prognostic model. (A,B) Overall survival (OS) in the external data sets (GSE44001, GSE52903). (C,D) ROC curve analysis in the external data sets.

The samples were divided into high and low risk groups based on the risk core. Univariate and bivariate analyses of CESC patients showed that the risk score was an independent prognostic factor, as shown in Fig. 8A, B. Nomograms were created to present the results of the regression analysis, the different staging patterns of cervical cancer were significantly correlated with the distribution of risk score values obtained by our model, shown in Fig. 9A. At the same time, the three-year and five-year OS of CESC patients were also predicted and analyzed, and the results showed that there was little difference between the predicted OS and the observed OS, suggesting that the nomogram model had a good prediction effect (Fig. 9B).

Fig. 8.

The model has an independent prognostic value. (A) Univariate analysis. (B) Multivariate Cox analyses.

Fig. 9.

Nomogram. (A) Indicators of clinical risk are correlated with the risk score. (B) The model accurately predicted 3- and 5-year OS of CESC.