Our study for publicly available CESC gene-expression data including 291 patients was yielded on TCGA, which downloaded from the UCSC Xena browser (GDC hub: https://gdc.xenahubs.net), and analyzed using the R software (version 3.6.1) and R Bioconductor packages.

The unsupervised clustering method of assessed metabolic genes was employed to classify patients into multiple clusters for further assessment by using the ConsensusClusterPlus R package. Then the value for k, where the cophenetic correlation coefficient magnitude began to fall was selected as the optimal cluster number (Hartigan and Wong, 1979). This analysis has been confirmed the classification stability for repeating 1,000 times (Monti et al., 2003).

The consensus ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression) algorithm with ESTIMATE R package was employed to measure ESTIMATE, immune and stromal scores, which reflected the immune and stromal cell gene signatures enrichment (Yoshihara et al., 2013).

Single-sample GSEA (ssGSEA) with GSVA R package was used for estimating immune infiltration in different clusters, and then an enrichment score indicated the extent to which genes were coordinately up or down-regulated within a single sample (Barbie et al., 2009).

Next, DEGs among the CESC METAclusters were identified using the R limma package. Only the genes with | logFC| > 1 (adjusted p < 0.01) were regarded as DEGs. Based on the above differential genes, genes with significant prognostic value were utilized for gene clustering by using the ConsensusClusterPlus R package.

Principal component analysis (PCA) was done and PC1 was selected as the signature score. After acquiring the prognostic value of each gene biosignature score, the following formula was used to describe the METAscore of every CESC patient: METAscore  =   ∑ PC 1 i - ∑ PC 1 j (1) which i is the signature score of clusters with positive Cox coefficient, and j is the expression of genes with negative Cox coefficients.

Gene set variation analysis (GSVA) is a unsupervised and nonparametric gene set enrichment approach that estimates biosignature scores or pathways based on transcriptomic data (Hänzelmann et al., 2013). We downloaded the gene sets from MSigDB (Broad Institute) (Subramanian et al., 2005), and chose gene ontology (GO) gene sets to quantify pathway activity. Pathway analysis was processed based on METAclusters and METAscore by using the GSVA R package in this study.

The TIDE (tumor immune dysfunction and exclusion) algorithm was applied to predict the potential response to immune checkpoint blockade (ICB) therapy of METAscore. For the melanoma dataset (GSE78220, N = 28), expression patterns (FPKM normalized) and phenotypes were transformed into transcripts per kilobase million (TPM), converting the FPKM data to values more comparable between samples to calculate METAscore (Wagner et al., 2012).

The normality of the variables was verified using the Shapiro-Wilk normality test (Ghasemi and Zahediasl, 2012). For comparisons between two groups, statistical significance was estimated using the unpaired Student t-tests and Wilcoxon tests for normally distributed variables and non-normally distributed variables, respectively. For comparisons between more than two groups, Kruskal-Wallis tests and one-way analysis of variance (ANOVA) were employed as nonparametric and parametric techniques, respectively (Hazra and Gogtay, 2016). Pearson and distance correlation analyses were conducted for correlation coefficients. Two-sided Fisher exact assessments were conducted to examine contingency tables. Based on dichotomized METAscore, patients were grouped into high and low METAscore groups and R package sva was employed to diminish the computational batch effect. These data were all visualized using the ggplot2 package in R. Benjamini–Hochberg method was used in the differential gene analysis to identify significant genes by converting the p values into FDRs (Benjamini and Hochberg, 1995). OncoPrint was applied to depict the mutation landscape of the TCGA dataset using maftools package in R (Gu et al., 2016). Cluster survival curves were generated using the Kaplan-Meier evaluation, and the log-rank test was employed to determine the differences statistical significance. Hazard ratios were computed using the univariate and multivariate Cox proportional hazards regression models. Independent prognostic factors were determined using the R survival package. Survival curves were generated using the survminer package. Heatmaps were generated using pheatmap function (https://github.com/raivokolde/pheatmap). All statistical and computational analyses were conducted using R programming (https://www.r-project.org/). These tests were two-sided and p value less than 0.05 signified statistical significance.

A flow diagram for the steps of this study was presented in Supplementary Figure S1. The 2,752 metabolic genes, encoding all human small molecule transporters and metabolic enzymes obtained from literature screening, were subjected to metabolism-related cluster classification by unsupervised clustering in the 291 CESC samples from TCGA. After assessing, a total of 315 candidate genes were sorted out, and clustering of the TCGA CESC samples was performed based on these genes using the ConsensusClusterPlus package in R. The optimal k number was determined as, compared with k = 2 or 4, the consensus matrix heatmap presented distinct and sharp boundaries when k = 3, supporting the robust and stable sample clustering. Thus, three distinct METAclusters were determined that 90 cases in METAcluster 1, 168 cases in METAcluster 2 and 33 cases in METAcluster 3 (Figure 1A; Supplementary Figures S2A–D). Survival analysis revealed the significant difference in patient overall survival (OS) time among these METAclusters, hinting the prognostic value in CESC (Figure 1B).

ESTIMATE is a tool using gene expression data to predict tumor purity and the presence of tumor immune/stromal cell infiltration. The ESTIMATE algorithm mainly generates three score-patterns to quantify the overall infiltration: 1) an ESTIMATE score that signifies tumor purity, 2) an immune score that infers the invasion of immune cells, and 3) a stromal score that denotes the presence of stromal cells. Significant differences in ESTIMATE and immune score, but not stromal score, were presented among the three METAclusters (Figure 1C). We next evaluated immune infiltration to describe their immune landscape. An abundance of 28 immune-correlated cell populations was computed using the ssGSEA. In accordance, result showed an obvious differential expression in immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, immature dendritic cells, macrophage, mast cell, MDSC, neutrophils, monocyte, and T helper cell) among the METAclusters. These data distinctly indicated these three METAclusters maintained different immune-relevant signatures (Figure 1D).

With the remarkable difference in immune characteristic identified, to further typify the transcriptomic and metabolic behavior differences among these metabolic patterns, we applied the GSVA enrichment analysis. Pathway analysis revealed that key carcinogenic activation pathways in CESC including WNT, HIPPO, NOTCH, NF-κB, and TGFβ (Maliekal et al., 2008; Ramos-Solano et al., 2015; Zhu et al., 2016; Bahrami et al., 2017; Tilborghs et al., 2017; Rodrigues et al., 2019; Wang et al., 2020) (Figure 2A), and metabolic pathways (Figure 2B) were differentially activated among these METAclusters which emphasized the genetic significance of the metabolism-based classification.

Then subsequent analysis investigated the expression of key immune checkpoints which have been selected based on current clinical trials drug inhibitors in other specific cancer types. As shown, this analysis revealed discriminable expression in ligand, cell adhesion, co-inhibitor, co-stimulator, antigen present, receptor and other checkpoints (Figures 2C–F, Supplementary Figures S2E–G). Considering these immune checkpoints were in charge of regulating the immune activation through modulating the T-cell in the immune respond process, we inferred that in CESC, respective METAcluster possessed different immune checkpoint blockade efficacy presumably.

To affirm metabolism-phenotype distinction of each METAcluster, unsupervised cluster analysis of 207 most representative DEGs among three METAclusters obtained using the limma package (Smyth, 2004) was completed to separate CESC patients into genomic subtypes (Figure 3A). The optimal cluster number supported the existence of three distinct and robust geneclusters in CESC patients (Supplementary Figure S3). Among these three geneclusters, the prominent difference in OS was strikingly consistent with the result of METAclusters (Figure 3B). Also, the expression of ESTIMATE and immune scores (Figure 3C), immune infiltration (Figure 3D), as well as the key immune checkpoints expression (Figure 4) were all highly in accordance with the differences among the METAclusters, genomically verifying three distinct metabolism-associated patterns in CESC.

Given the individual complexity and heterogeneity of metabolic modification, we used the PCA algorithm to construct the METAscore to quantify metabolic patterns in CESC patients. Two aggregate score groups (high- and low- METAscore groups) were generated by the sum of relevant scores, and GSVA analysis uncovered that the METAscore was related to the immune signaling pathways, cancer pathways (Figure 5A), and key metabolic pathways (Figure 5B).

Then we evaluated the potential of the METAscore to predict CESC survival. Univariate and multivariate Cox regression model analysis, which considered including patient age, stage, grade, pregnancy count, radiation therapy and targeted therapy, confirmed that the METAscore was an independent and reliable prognostic biomarker (Figures 5C,D). Besides, the prognostic significance of the METAscore was measured in four independent gynecological cancers including CESC, breast cancer (BRCA), ovarian cancer (OV) and endometrial cancer (UCEC) (Figure 5E). Notable OS differences emerged between the high- and low-METAscore groups in BRCA, OV and CESC which were all recognized as hot tumors with distinct T-cell invasion (Figures 5F–I).

Accordingly, the next evaluation was concerned in immune checkpoints expression between two METAscore groups. Robust correlation between METAscore and different response of immune checkpoints including receptor, ligand, cell adhesion, co-inhibitor, antigen present and other checkpoints was demonstrated, indicating the guiding role in immunotherapy of METAscore (Figure 6).

To determine the difference in somatic mutation frequencies between the high- and low- METAscore groups, we analyzed the TCGA genomic files. The consequence revealed that the high- and low- METAscore groups exhibited distinct mutation characteristics and the genes with a high mutation frequency in TTN, MUC4, PIK3CA, and MUC16 which all correlated with EMT (Chen et al., 2018) and critical cancer pathways including PI3K/AKT (Razia et al., 2019) and JAK2/STAT3 (Shen et al., 2020) (Figures 7A,B) in CESC. Somatic mutations in the PIK3CA denoted a late event during cervical carcinogenesis, and have been detected in multiple cervical carcinoma subgroups (Verlaat et al., 2015). Besides, MUC4 and PIK3CA were frequently mutated in the HPV16-KRT, a HPV16 subtype typified by increased expression of keratinization genes, biological oxidation and Wnt pathway signaling (Lu et al., 2019). Similarly, regarding altered somatic copy number variation (CNV), remarkable differences in driver gene amplification and deletion were emerged between the METAscore groups (Supplementary Figure S4). These analyses revealed a high genomic heterogeneity and altered gene expression landscape during the METAscore groups, supporting the association between the METAscore and genomic alterations.

Immune checkpoint blockade (ICB) therapy is a revolutionary immune-based treatment in cancers including CESC. Inhibition of the immune checkpoints using monoclonal antibodies that block the T-cell molecules PD-1, PD-L1, as well as CTLA4 has emerged as a novel anti-cancer treatment with extraordinary survival advantages (Curran et al., 2010). Considering the correlation between the METAscore and immune characteristics, we elaborated the predictive potential of the METAscore to estimate ICB therapeutic value by using the melanoma GSE78220 cohort. TIDE algorithm is a method of modeling two primary mechanisms of tumor immune infiltration: the stimulation of T-cell dysfunction companying with high cytotoxic T-lymphocytes (CTL) infiltration, and the prevention of T-cell infiltration with low CTL levels, which estimates potential response to immunotherapy (Wang et al., 2020b). We conducted TIDE algorithm and obtained that patients in high- METAscore group tended to respond to immunotherapy, prompting CESC patients with high- METAscore might more likely benefit from immunotherapy (Figures 7C,D). Combined with prediction of survival outcomes in CESC (Figure 7E), we assured the guiding value of METAscore in immunotherapy.