Construction and validation of a metabolic gene-associated prognostic model for cervical carcinoma and the role on tumor microenvironment and immunity

Identified differentially expressed metabolic genes and conducted enrichment analysis

The cancer cell metabolic gene set comprising 2071 genes was downloaded from ccmGDB (https://bioinfo.uth.edu/ccmGDB/). A total of 495 CC-specific mRNAs (1320 upregulated and 1944 downregulated) were identified as differentially expressed metabolic genes (DEMGs), with Limma package using P value < 0.05 and |logFC| > 1 as screening criteria. The volcano map exhibits significant differences and distribution of mRNAs in CC based on the fold change (Figure 1A). GO analysis showed that the DEMGs were significantly associated with purine-containing compound metabolic, ribose phosphate metabolic, purine nucleotide metabolic, ribonucleotide metabolic and purine ribonucleotide metabolic processes (Figure 1B). And the pathways such as purine metabolism, biosynthesis of amino acids, pyrimidine metabolism, biosynthesis of cofactors and carbon metabolism were mainly enriched in the KEGG pathways (Figure 1C).

Figure 1. Filtered differentially expressed metabolic genes (DEMGs) and performed Enrichment analysis. (A) A volcano map presented DEMGs between the cervical cancer and the normal samples, Dots represent significantly down-regulated (Green) and up-regulated (Red) genes. (B) GO function annotations and (C) KEGG enrichment analysis.

Construction of a metabolic prognostic model based on fifteen DEMGs for cervical carcinoma

The survival difference of the DEMGs was analyzed by Kaplan-Meier method and log-rank analysis. There are 72 genes remarkably related to the prognosis of CC (Supplementary Figure 1). A total of 15 DEMGs were incorporated into the multivariate Cox regression analysis after univariate Cox and Lasso regression analysis (Figure 2A–2C). A multivariate Cox proportional hazards regression model was subsequently utilized to build the metabolic gene-associated prognostic model (MGPM) risk signature and calculate the risk score of each individual, patients were then divided into low- and high-risk groups using the median risk score as a cut-off. And nine out of fifteen variables ((HCCS, LDHC, PGK1, MSMO1, PLA2G7, LIPG, TUBB4B, AGPAT4 and GNG8) had significant statistical difference (P < 0.05) in MGPM (Supplementary Table 1). The survival analysis of the MGPM using the Kaplan–Meier method indicated that the patients in the high-risk group had a shorter overall survival compared with those in the low-risk group, and there was a significant difference between them (P < 0.001, Figure 2D). The ROC curve analysis represented a favorite prognostic value for this model in evaluating the overall survival of patients with CC (1-year AUC = 0.842, 3-year AUC = 0.861, 5-year AUC = 0.849; Figure 2E). The survival status of each patient was evaluated according to the risk score (Figure 2F). Compared with the low-risk group, patients have a shorter survival time and a higher incidence of death in the high-risk group (Figure 2F). The percentage of deaths in the low- and high-risk groups was 42% and 8%, respectively. There is a statistical difference between them (P < 2.22E-16, Figure 2G). Additionally, the nomogram was employed to depict the prospective survival of CC patients (Figure 2H), and the calibration plot was performed to assess the discernment of the MGPM-based risk signatures, and the result exhibited its predicted value was virtually coincided with the actual value (Figure 2I). The clinicopathological characteristics of CC patients in this model are detailed in Supplementary Table 2.

Figure 2. Construction of metabolic gene-associated prognostic model (MGPM) using a series of bioinformatics technology. (A, B) LASSO Cox regression analysis with 5-fold cross-validation and generation of coefficient outline based on the log (lambda) sequence were used to identify the potential independent prognostic risk signature genes. (C) Multivariate Cox regression analysis results show the p values and hazard ratios (HR) with confidence intervals (CI) of the fifteen DEMGs. (D) Kaplan-Meier survival curves show the overall survival (OS) rates of high-risk (n = 146) and low-risk (n = 147) CC patients. Patients in high-risk group had a shorter OS compared to those in low-risk group. (E) ROC curve analysis results show the accuracy and reliability of the MGPM in determining the 1-, 3- and 5-year survival outcomes (AUC values are shown in parentheses). (F) From the top to bottom, the survival status of each patient was sorted according to the low-risk (blue) and high-risk (red) scores. A scatter plot then exhibited the survival status and survival time of patients, the dots represent patients that have died (red) and alive (blue) at the time of analysis. (G) Analysis of the difference in survival status between the low- and high-risk groups. (H, I) The nomogram and calibration plot were employed to depict the expected survival of individual CC patients and to assess the discriminative ability of the MGPM-based risk signature, respectively.

Validation of the MGPM-based risk signature in the training and testing cohorts

Subsequently, we proceeded with internal verification in the training (N = 147) and testing (N = 146) groups to evaluate the accuracy and robustness of the MGPM comprising fifteen DEMGs (Figure 3). From top to bottom, the risk score of each sample was calculated and ranked based on the MGPM in the training and testing data groups, respectively. And then, the overall survival status of CC patients according to the risk score distribution was presented using a scatter plot. Both data sets showed patients with high-risk score were bound up with a higher mortality rate, and the increase of risk score was proportional to the death rate in CC patients. A heatmap presenting the expression profiles of those signature genes displayed those patients in the high-risk group were more geared to elevate NPR2, PGK1, MSMO1, LIPG and AGPAT4 levels, while those in the low-risk group were more inclined to elevate ALOX12B, HCCS, LDHC, PTGDS, SULT1E1, SULT1A3, TYMS, PLA2G7, TUBB4B and GNG8 levels, these finding in the training dataset showed consist with the testing dataset. Lastly, Kaplan-Meier survival curve analysis indicated that the overall survival of patients with high-risk score was remarkably shorter than those with low-risk score in both datasets (P < 0.001), and ROC curve analysis also indicated there was a favorite prognostic predictive value of the MGPM both in training set (1-year AUC = 0.842, 3-year AUC = 0.861, 5-year AUC = 0.849) and testing set (1-year AUC = 0.843, 3-year AUC = 0.85, 5-year AUC = 0.827).

Figure 3. The MGPM comprising fifteen metabolic genes was verified in the training and testing data set. From top to bottom, the ranked dot plot illustrating the predictor-score distribution, a scatter plot presenting the patients’ overall survival status, a heatmap showing the expression profile of the fifteen signature genes of CC patients, Kaplan-Meier survival curves indicating the OS rates of high- and low-risk groups and ROC curve analysis demonstrating the accuracy and reliability of the MGPM in determining survival outcomes.

Relationship between MGPM-based risk signatures and clinicopathological characteristics of CC patients

Next, the association between the MGPM-based risk signatures and the clinicopathological parameters was analyzed. A heatmap was employed to present the correlations between the MGPM-based risk signatures and clinicopathological variables. The results indicated that there were significant statistical differences between them regarding the tumor node metastasis classification (P < 0.01), survival status (P < 0.001) and estimate score (P < 0.0001) (Figure 4A). The differences in clinicopathological parameters were then evaluated between the low- and high-risk groups of MGPM by Wilcoxon rank-sum test. There are statistically significant differences in stage (P = 0.009) primary tumor (T, P = 0.00056) and lymph nodes (N, P = 0.031), while no statistical differences in remote metastasis (M, P = 0.97) (Wilcoxon signed-rank test, Figure 4B–4E).

Figure 4. Correlation analysis between the MGPM and clinicopathological features in low- and high-risk CC patients. (A) A heatmap showed the correlation between the MGPM-based risk score and tumor microenvironment (ESTIMATEScore) or clinicopathological parameters using the Chi-square test. There were statistically significant differences in ESTIMATEScore, survival status and tumor node metastasis classification (T, N and M stages) of low- and high-risk groups. (B–E) The distribution of risk scores based on MGPM stratified according to (B) stage (stages I–II vs. stages III–IV), (C) T stage (T1-2 vs. T3-4), (D) M stage (M0 vs. M1) and (E) N stage (N0 vs. N1-2).

To explore its clinical application and usability, we performed Cox regression analysis on common clinical characteristics. The results exhibited that the risk score of MGPM had statistically significant differences in both univariate- and multivariate-Cox regression groups (P < 0.001), which indicated the MGPM could as an independent index for prognosis to forecast the OS of CC patients (Figure 5A, Supplementary Table 3). After stratifying patients according to clinicopathological parameters, the overall survival rate of patients with the following characteristics will be significantly shortened in the high-risk group by Kaplan–Meier analysis, including age ≤ 65 (P < 0.001) and age > 65 (P = 0.022) (Figure 5B), stage I/II (P < 0.001) and stage III/IV (P < 0.001) (Figure 5C), T1-2 stage (P < 0.001) and T3-4 stage (P = 0.006) (Figure 5D), M0 stage (P < 0.001) (Figure 5E (a)), N0 stage (P = 0.006) and N1 stage (P = 0.001) (Figure 5F). However, there were no significance in M1 stage patients (P = 0.208, Figure 5E (b)) in both the high- and low-risk groups.

Figure 5. Cox regression and Kaplan-Meier survival analysis of risk score and clinicopathological parameters. (A) Cox regression analysis of risk score and clinicopathological parameters, which indicate the MPGM could as an independent prognostic factor. (B) Patients stratified by ages (age ≤ 65 and > 65), (C) Patients stratified by stage (stage I-II and stage III/IV), (D) Patients stratified by T stage (T1-2 and T3-4), (E) Patients stratified by M stage (M0 and M1) and (F) Patients stratified by N stage (N0 and N1).

The risk score based on MGPM was negatively associated with ESTIMATEScore of tumor microenvironment

ImmuneScore and StromalScore have been represented scores according to the ratio of immune and stromal components in the tumor, respectively. ESTIMATEScore was considered as the summation of them, and the tumor purity of each subject was calculated based on all three scores. The overall survival of patients was positively correlated with ESTIMATEScore (P = 0.0026, Figure 6A), while not significantly with StromalScore and ImmuneScore (Supplementary Figure 2). Then, the data combined with the clinical information and ESTIMATEScore of CC patients were analyzed to discover the difference between them, showing that ESTIMATEScore were notably declined in M1 group compared to M0 group with a statistical significantly difference (Wilcoxon signed-rank test, P = 0.041, Figure 6B). These findings indicated that the components of immune and stromal were involved in the progression of CC, especially in metastasis. In the constructed MGPM, the ESTIMATEScore was higher in the low-risk group than the high-risk group with a statistically significant difference, indicating a better prognosis in CC patients (Wilcoxon signed-rank test, P = 3E-05, Figure 6C). In addition, the Kaplan–Meier survival curves suggested the patients with low-ESTIMATEScore and high-risk score (based on the MGPM) had a shorter survival time and a higher incidence of death (P < 0.001, Figure 6D).

Figure 6. The relationship between the MGPM-based risk score and tumor microenvironment and immune cell infiltration. (A) Kaplan-Meier survival curve was utilized to evaluate the prognosis of the tumor microenvironment, the result showed low-ESTIMATEScore had a shorter OS. (B) The distribution of ESTIMATEScore according to clinicopathological features (M). (C) The difference of ESTIMATEScore in low- and high-risk groups. (D) Kaplan-Meier survival analysis of ESTIMATEScore in combined with a risk score. (E) Violin plot showed the ratio differentiation of 21 types of immune cells between CC tumor samples in low- or high-risk groups, and Wilcoxon rank-sum was applied for the significance test. (F) Prognosis-related immune-infiltrating cells. (G) The correlation analysis between 21 types of immune-infiltrating cells and risk score, ESTIMATEScore and MGPM-based risk signature genes. (H) Immunotherapy score of low- and high-risk groups based on TCIA database.

Relationship between the MGPM and the proportion of tumor-infiltrating immune cells

The proportion of tumor-infiltrating immune subtypes was analyzed by CIBERSORT method in CC cases (Supplementary Figure 3). Subsequently, we evaluated the important tumor-infiltrating immune cell differences in low- and high-risk groups through the Wilcoxon rank-sum test (Figure 6E). There are ten types of tumor-infiltrating immune cells with statistically significant differences, including T cells CD8 (P < 0.001), T cells memory activated (P = 0.010), T cells regulatory (Tregs) (P = 0.008), Macrophages M0 (P < 0.001), Macrophages M1 (P = 0.018), Dendritic cell activated (P = 0.005), Mast cells resting (P < 0.001), Mast cells activated (P < 0.001), Eosinophils (P = 0.007) and Neutrophils (P = 0.006). Among them, T cells CD8 (P = 0.011), T cells memory activated (P = 0.011), Mast cells resting (P = 0.029) and Mast cells activated (P < 0.001) were correlated prognosis in CC patients (Figure 6F). The co-expression analysis of twenty-one tumor-infiltrating immune cells and the MGPM-based risk signature was performed (Figure 6G). The results indicated T cells CD8 (R = −0.25, P = 2.37E-05), Macrophages M0 (R = 0.27, P = 2.24E-06), Mast cells resting (R = −0.29, P = 4.05E-07) and Mast cells activated (R = 0.35, P = 8.23E-10) were correlated with MGPM-based risk signature with statistically significant differences.

The co-expression analysis of twenty-one tumor-infiltrating immune cells and the MGPM-based risk signature and tumor microenvironment was then performed (Figure 6G). The results indicated T cell CD8 had a positive correlation with ESTIMATEScore (R = 0.34, P = 1.51E-09), PLA2G (R = 0.29, P = 2.8E-07) and GNG8 (R = 0.17, P = 0.003) expressions, a negative correlation with risk score of MGPM (R = −0.18, P = 0.002), LIPG (R = −0.21, P = 0.0003) and MSMO1 (R = −0.21, P = 0.00024) expressions. T cells regulatory (Tregs) had a positive correlation with the expressions of PGK1 (R = −0.17, P = 0.003), MSMO1 (R = −0.14, P = 0.014), LIPG (R = −0.16, P = 0.004), ALOX12B (R = −0.12, P = 0.041) and HCCS (R =−0.12, P = 0.045). Mast cells resting had a positive correlation with the LDHC expression (R = 0.15, P = 0.009), and a negative correlation with risk score (R = −0.19, P = 0.0012), LIPG (R = −0.19, P = 0.001) and PGK1 (R = −0.23, P = 4.25E-05) expressions. Mast cells activated had a positive correlation with risk score (R = 0.36, P = 1.2E−10), LIPG (R = 0.34, P = 2.32E-09) and PGK1 (R = 0.17, P = 0.003) expressions, a negative correlation with ESTIMATEScore (R = −0.16, P = 0.0049).

Correlation of MGPM with the response of immune checkpoint inhibitors

To assess the relationship between MGPM and immunotherapy responses, we explored the correlation of MGPM-based risk score with common immune checkpoint inhibitors based on the TCIA database, including cytotoxic T lymphocyte antigen 4 (CTLA4) and programmed cell death 1 (PD1). The results indicated that patients in low-risk group had better responses to immunotherapy (Wilcoxon signed-rank test, P < 0.05, Figure 6H).

Data selection and extraction of differentially expressed genes

RNA profiles comprising HTseq-count and fragments per kilobase per million reads (FPKM) of 306 primary CCs and 3 normal subjects were downloaded from the cancer genome atlas database (TCGA, https://portal.gdc.cancer.gov/, Supplementary Table 2). The corresponding clinical information, including general information (age, race, ect), survival outcome (survival status and days to the last follow-up) and tumor feature (stage, lymph nodes, metastases, etc) were also downloaded. After filtering the non-CC specific expression genes, the edgeR package (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) was applied to obtain differentially expressed metabolic-related genes using Padj < 0.05 and |logFC| > 1 as screening criteria.

Function and pathway enrichment analysis of the DEMGs

A rounded enrichment analysis concerning the functions and pathways of the DEMGs was performed in CC by the clusterProfiler packages. The cellular components (CC), biological processes (BP), and molecular functions (MF), as well as the KEGG signaling pathways. GO terms and KEGG pathways were presented by the bubble chart using the P-value < 0.05 as a screening criterion.

Identification of prognosis-associated metabolic genes and construction of risk model

The mRNA expression data of normal and para-carcinoma tissue samples, and with incomplete clinical information were deleted by using the Limma package. Univariate Cox and lasso regression analyses were identified the potential prognostic metabolic-related genes, and then multivariate Cox regression analysis was performed to construct metabolic gene-associated prognostic models (MGPM) and to assure that the multi-factor models were not overfitted. And the risk score of each sample was calculated through the following formula: risk score = [Coefi1] × Exp1 + [Coefi2] × Exp2 + ... [Coefin] × Expn, Coefi is the risk coefficient of each subject deduced from the LASSO-Cox mode, Exp is the expression of consensus genes. the receiver operating characteristic (ROC) curve was utilized to measure the discriminative ability of MGPM. A nomogram was formulated to present the risk score of each individual and the calibration plot was applied to evaluate the performance of the MGPM.

CIBERSORT estimation

The CIBERSORT algorithm was utilized to evaluate the expression abundance of 22 tumor-infiltrating immune cells (TICs) in CC using P < 0.05 as a cut off. And then, Wilcoxon rank-sum test was executed to calculate significant differences of the TICs in the proportion between low- and high-risk groups in MGPM. Kaplan–Meier survival analysis was then applied to explore the prognosis of all TICs using P < 0.05 as a threshold.

Statistical analysis

Statistical analyses were carried out by R software, with version number v3.5.2 (Package: GDCRNATools, Limma, ggplot2, rms, preprocessCore, glmnet, survminer, timeROC). All tests were two-tailed and P < 0.05 was supposed to have a statistical significance.

Data availability statement

Publicly available dataset (TCGA, https://portal.gdc.cancer.gov/) was analyzed in this study.