RNA-sequencing and matched clinicopathological data (including histological grade, sex, stage, age, and survival data) were retrieved from the Cancer Genome Atlas (TCGA) database, the Cancer Genomics Browser of University of California Santa Cruz (https://genomecancer.ucsc.edu/).

The feature gene panels of cholesterol catabolism process, fatty acid metabolism, glutamine metabolic process, glycolysis, one carbon metabolic process, and pentose phosphate process were obtained from the molecular signature database (MsigDB). The enrichment scores of the six metabolic processes of each patient were analyzed by single sample gene set enrichment analysis (ssGSEA). ssGSEA estimates the relative enrichment of the six metabolic process gene sets in each sample by comparing the gene expression data of each sample with a specific gene set, which was downloaded from MsigDB.

The clustering method that we utilized in this research is called consensus clustering analysis. It is an unsupervised clustering, which use resampling method to extract data. For each sampling, the number of clusters was specified. Then, the process calculates the rationality of different cluster numbers using PAC method and identify and produce the most reasonable clustering result. In our research, the consensus clustering analysis was performed to classify patients into several bio-similar groups by multi-dimensional scaling, which classify samples into several subtypes according to six metabolic processes enrichment scores, so as to detect new disease subtypes with similar metabolic features.

Mutant allele tumor heterogeneity (MATH) algorithm is a reliable and applicable method to measure intra-tumor heterogeneity, which has been validated and reported in gastric cancer (13). The calculation method of MATH was identified at the Broad Institute of MIT and Harvard. The calculation method includes three steps. First, each difference value of the mutant-allele fraction (MAF) was obtained from the median difference value. Second, the median absolute deviation (MAD) in R was calculated as values scaled by a factor (1.4826) to make the expected MAD of a sample from a normal distribution equals the standard deviation. Third, MATH was figured up as MATH = 100 * MAD/median. Maftools package in R was utilized to calculate MATH. Maftools package includes clustering algorithm to improve the accuracy of genomic pattern.

Gene set enrichment analysis (GSEA) was used to analyze signaling pathway enrichment among different metabolic bio-similar groups (14). The GSEA analysis was performed according to the following standards: 1) the nominal p-value is lower than 0.05; 2) the value of NES is 1 or more than 1; 3) false discovery rate (FDR) q-value is lower than 0.25. Pathways that meet the above three criteria were considered statistically significant. The gene sets of signaling pathways were downloaded from MSigDB database (14).

The 28 immune pathway activations of patients were analyzed by single-sample gene set enrichment analysis (ssGSEA). The activation of the immune pathways were identified by a feature panel of genes overexpressed in each kind of immune cell (15, 16). The 28 immune pathway scores include activated CD4 T cell, activated CD8 T cell, central memory CD4 T cell, central memory CD8 T cell, effector memory CD4 T cell, effector memory CD8 T cell, type 1 T-helper cell, activated dendritic cell, natural killer cell, regulatory T cell, type 2 T helper cell, macrophage, myeloid-derived suppressor cells (MDSC), activated B cell, gamma delta cell, immature B cell, T follicular helper cell, eosinophil, mast cell, and monocyte. The differences in immune pathway activation among the metabolic bio-similar groups were analyzed.

RNA-seq and clinical data were retrieved from GEO database. GSE15459 dataset was downloaded from the GEO (https://www.ncbi.nlm.nih.gov/gds/) on the GPL570 platform. The classification of patients by metabolic bio-similar features was also performed.

Statistical methods were all analyzed using R version 3.5.1 (http://cran.r-project.org) and Stata statistical software, version 12.0 (StataCorp, College Station, TX). Patient baseline characteristics were compared among the different metabolic bio-similar groups via Wilcoxon rank-sum test, while multiple comparisons were performed using Kruskal–Wallis test. The association between patient characteristics and metabolic clusters was analyzed using the Spearman test. Kaplan–Meier curves were presented between different metabolic bio-similar groups for OS analysis. Log-rank test was utilized to compare the overall survival between different groups. Multivariate analysis of prognostic predictors was carried out using a Cox proportional hazards model. The differences in immune pathway activation among the metabolic bio-similar groups were analyzed using Kruskal–Wallis test. The difference in MATH among the metabolic bio-similar groups were analyzed using Wilcoxon rank sum test. A two-sided p-value <0.05 was considered significant.

A total of 370 GC patients with complete information from TCGA database were enrolled in the study as training cohort. The clinical and pathological characteristics are illustrated in Table 1. The enrichment scores of the six metabolic processes of each sample were analyzed by single sample gene set enrichment analysis (ssGSEA). As we have calculated the enrichment score of the six metabolic processes for each sample, then the consensus clustering analysis, an unsupervised analysis, was performed to classify patients into several bio-similar clusters according to six metabolic processes enrichment scores, so as to detect new clusters with similar metabolic features. In our study, four clusters were identified. Figure 1 illustrated the six metabolic process enrichment of each patient and the four clustering of the patients according to metabolic process enrichment. Patients in different clusters have prominent characteristics. In cluster 3, except for cholesterol catabolism, all metabolism processes were activated. In cluster 2, all metabolism processes were activated, and cholesterol catabolism was activated only in cluster 2. In cluster 4, all metabolism processes were in low expression levels. In cluster 1, none of the metabolic processes are specifically activated or silenced. Then, we named the four clusters according to metabolic characteristics. Cluster 1 was named as metabolic non-specificity cluster. Cluster 2 was named as metabolic-active cluster. Cluster 3 was named as cholesterol-silence cluster. Cluster 4 was named as metabolic-silence cluster.

There were no significant relationships between metabolic subtypes and patient characteristics, including sex, age, stage, and histological grade (Table 1).

Survival curves were estimated by Kaplan–Meier method and analyzed by log-rank test to assess the survival difference between the four metabolic subtypes. In the total of 370 patients, patients in metabolic non-specificity cluster and cholesterol-silence cluster had significantly longer OS than patients in metabolic-active cluster and metabolic-silence cluster (p = 0.031, Figure 2A). The median OS in the four metabolic subtypes were 57.4, 21.0, 55.4, and 28.7 months, respectively. We identified metabolic non-specificity cluster and cholesterol-silence cluster as good-prognosis group, while metabolic-active cluster and metabolic-silence cluster were identified as poor-prognosis group. Patients in good-prognosis group had significantly longer OS than patients in good-prognosis group (p = 0.007, Figure 2B). The multivariate Cox proportional hazards model showed that metabolic bio-similar clustering, which was categorized as good- and poor-prognosis groups, was an independent prognostic factor for OS (hazard ratio, 0.622; 95% CI, 0.441–0.879; p = 0.007, Table 2), after adjusting for clinicopathological characteristics, including age, sex, pathological stage, and histological grade.

Intra-tumor heterogeneity was measured using mutant allele tumor heterogeneity (MATH) algorithm. MATH score is significantly associated with overall survival in gastric cancer patients, of which the result has been published. The Wilcoxon rank sum test showed that patients in metabolic-active cluster and metabolic-silence cluster had significantly higher MATH score than those in metabolic non-specificity cluster and cholesterol-silence cluster (Figure 2C, p = 0.032).

We think it is necessary to analyzed the metabolic feature of normal gastric tissues using the same method, in order to verify that tumor do affect the metabolism of gastric cells. Then, we analyzed 35 normal gastric tissue from TCGA database (Figure 2D). None of the metabolic processes were significantly activated or inhibited in normal gastric tissue. There was no specific metabolic characteristics in normal gastric tissues.

In metabolic-active cluster and cholesterol-silence cluster, the patients’ outcome was different, while the metabolic characteristic was in fact similar, except for cholesterol catabolism. Patients in cholesterol-silence cluster had significantly longer OS, which was of most concern in further analysis. The difference in signaling pathways between different clusters was analyzed using GSEA. Three cell-cycle-related pathways are significantly downregulated in cholesterol-silence cluster than metabolic-active cluster, which were G2M checkpoints, E2F targets, and MYC targets (Figure 2E).

The 28 immune pathway activations of patients were analyzed by ssGSEA. Patients in cholesterol-silence cluster had significantly longer OS than patients in metabolic-active cluster and metabolic-silence cluster. The Wilcoxon rank sum test showed that several anti-tumor immune pathway activation was significantly higher in cholesterol-silence cluster than in metabolic-active cluster, including activated CD4 T cell, type17 T-helper cell, and several pro-tumor immune pathway activations were significantly higher in metabolic-active cluster than in cholesterol-silence cluster, including macrophage, MDSC, and plasmacytoid dendritic cell (Figure 2F).

We enrolled 192 gastric cancer patients with complete RNA-seq and clinical information from GSE15459 dataset to establish validation cohort. The 192 patients were classified into four clusters according to metabolic bio-similar clustering in training cohort. The enrichment of the six metabolic processes of the patients and the four clusters are illustrated in Figure 3A. It is obvious that the metabolic features of the four clusters were consistent with that in the training cohort.

Survival curves were also estimated by Kaplan–Meier method and analyzed by log-rank test between the four metabolic subtypes. In a total of 192 patients, patients in metabolic non-specificity cluster and cluster 3 had significantly longer OS than patients in cluster 2 and cluster 4 (p = 0.030, Figure 3B). The median OS in the four metabolic subtypes were 99.4 months, 23.6 months, not reached, and 31.2 months, respectively. We identified metabolic non-specificity cluster and cholesterol-silence cluster as good-prognosis group, while metabolic-active cluster and metabolic-silence cluster were identified as poor-prognosis group. Patients in good-prognosis group had significantly longer OS than patients in good-prognosis group (p = 0.007, Figure 3C). The multivariate Cox proportional hazards model showed that metabolic bio-similar clustering, which was categorized as good- and poor-prognosis groups, was an independent prognostic factor for OS (hazard ratio, 0.468; 95% CI, 0.291–0.559; p = 0.002, Table 3), after adjusting for clinicopathological characteristics, including age, sex, and pathological stage.