We downloaded LGG patients’ gene transcriptome profiles and corresponding RNA-seq data from TCGA database (https://portal.gdc.cancer.gov/) and the Genomic Data Commons (GDC) tool. 509 LGG samples logged in TCGA were available for analysis. Additionally, mRNA-seq 693 (batch 1) and mRNA_seq325 (batch 1) datasets were downloaded from CGGA (http://www.cgcg.org.cn/). Next, these two RNA-seq datasets were merged into one metadata set; 408 samples were ultimately included in this study.

Previous studies by other researchers screened a set of metabolic genes relevant to malignant tumor activity and stemness properties required for tumorigenesis (Possemato et al., 2011). A total of 2,752 metabolism-relevant genes which encoded metabolic enzymes and transporters were selected.

Poor-quality samples were excluded before data preprocessing. Samples without clinical data and with more than 50% missing were removed. In the end, we obtained an expression profile dataset containing 20,485 genes profiles for subsequent analysis.

ConsensusClusterPlus implements the consensus clustering (CC) method, which facilitates more quantitative stability evidence in unsupervised class discovery (Wilkerson and Hayes, 2010). We used the Normalized Enrichment Score (NES) to measure the gene-sets enrichment (Wilkerson and Hayes, 2010; Yang et al., 2020), and obtained 113 normalized enrichment scores of metabolism-relevant gene signatures. Each of the 113 metabolism-related signatures had a class of gene sets and contained multiple genes. We adopted the “PAM” algorithm along with “Canberra” as a measure of distance, and performed 500 bootstraps, each involving 80% of the patients in the training set. The clustering number was set as 2–10, and the consistency matrix and consistency cumulative distribution function were calculated to determine the best classification.

To infer the relative abundance of 22 types of tumor-infiltrating immune cells and non-immune cells in the tumor microenvironment, we applied the CIBERSORT algorithm via the “CIBERSORT” R package (CIBERSORT R script v1.03; http://cibersort.stanford.edu/). We adopted the ESTIMATE package to determine the presence of infiltrating immune cells (Yu et al., 2019), using the ImmuneSignature gene set based on LGG RNA-seq data (Jusakul et al., 2017).

To further explore the differences between different clusters in biological processes, we performed the gene set variation analysis (GSVA) with the R package “GSVA” to estimate pathway enrichment for different clusters. In total, 113 metabolism-associated pathways were included in GSVA analysis genesets. R package clusterProfiler (https://guangchuangyu.github.io/software/clusterProfiler) was used to process the GSEA analysis.

Considering that different metabolic characteristics existed in different metabolic subtypes, we applied the principal components analysis (PCA). And we established subtype classification scores better to quantify metabolic characteristics of patients in different sample cohorts. Specifically, we used 113 metabolic characteristics for PCA analysis, and the first two-component scores. Then, we calculated the metabolic subtype characteristic score of each sample and calculate the formula MRGs-score = ∑(PC1 I + PC2 I), where I represented metabolic characteristics.

To identify prognosis-related metabolic signatures, we calculated enrichment score and prognosis features of 113 metabolic genes from LGG patients in the TCGA and CGGA cohorts. The results indicated that there were 69 prognosis-related metabolic signatures in the TCGA database, and 73 prognosis-associated metabolic genes (Figure 1A). Furthermore, we adopted ConsensusClusterPlus to explore reasonable classifications according to the characteristics of the 509 LGG sample, and classified these similar characteristics genes into one category. The corresponding cumulative distribution function (CDF) curve and the delta area plot indicated that the optimal choice was K = 3 (Figures 1B,C). Finally, we obtained three different metabolic clusters (MC) which were characterized by consistent metabolic subtype-related specificities (Figure 1D). We further investigated the prognosis linked to these three clusters. Results for overall survival (OS) indicated that MC1 predicted the shortest survival time compared with MC2 and MC3 in the TCGA dataset (Figure 1E). The progression-free survival (PFS) outcomes showed that the MC1 subgroup had the poorest survival rates in both the TCGA cohort and the CGGA cohort (Figures 1F,G).

To clarify the clinical signatures between the three clusters, the results demonstrated that there was no significant difference in occurrence of any signature based on age or gender (Supplementary Figures S1A,B). There was a remarkable variation in IDH mutation in the TCGA database, and IDH mutation in the MC2 and MC3 subgroups was significantly higher than that in the MC1 subgroup (Supplementary Figure S1C). The probabilistic infiltration weighted gradient maps also demonstrated that the chromosome 1p and 19q (1p/19q) codeletion status was significantly higher compared with MC2 and MC1 in the TCGA dataset (Supplementary Figure S1D). The methylguanine-DNA methyl-transferase (MGMT) expression in MC1 was obviously decreased compared with the MC2 and MC3 subgroups (Supplementary Figure S1E). At the same time, we explored these features in the CGGA database, and noticed that consistent with the prior observations, the three clusters did not differ significantly based on age or gender (Supplementary Figures S1F,G). IDH mutation and 1p19q co-deletion were significantly higher in MC2 and MC3 compared with MC1 (Supplementary Figures S1H,I). MGMT expression was not significantly different in the three subgroups (Supplementary Figure S1J).

To thoroughly investigate the metabolic characteristics, we applied differential expression analysis to identify subtype-specific metabolic characteristics based on the GSVA scoring system. According to the results, 47 specific metabolic characteristics were present in MC1, 6 in MC2, and 39 in MC3, as visualized by the heatmap (Supplementary Figure S2). These genes were applied to further distinguish between metabolic subgroups.

We also analyzed the differences in genomic changes among these three metabolic subtypes in the TCGA cohort. The MC3 subtype exhibited a lower aneuploidy score (Figure 2A), fraction altered (Figure 2B), number of segments (Figure 2C), and number of homologous recombination defects (Figure 2D). Tumor mutation burden was similar in all three subgroups (Figure 2E). In addition, we analyzed the correlation between gene mutations and metabolic subtypes, and found metabolic subtypes. There was a significant correlation with gene mutation. TP53, IDH1, ATRX, and EGFR genes showed extensive somatic mutations in LGG, among which the IDH1 gene had a higher mutation frequency in MC2 and MC3 subtypes, and patients with tumors with IDH1/2 mutations had favorable prognoses. The mutation frequency of the ATRX gene in the MC2 subtype was the highest, followed by MC1. EGFR had a higher mutation frequency in the MC1 subtype, but a lower frequency in both MC2 and MC3 subtypes. In terms of copy number variation, the MC1 subtype had a wide range of copy number amplification and deletion frequencies (Figure 2F).

To comprehensively investigate the tumor heterogeneity among different metabolic subtypes, we obtained genomic characteristics including tumor purity, ploidy, amplification, and intertumoral heterogeneity from known studies (Thorsson et al., 2018). They indicated that the purity, ploidy, and intratumor Heterogeneity of MC1 and MC3 were significantly lower than those of MC2 (Supplementary Figures S3A–C). MC1 showed the highest proliferation score, followed by MC2 and MC3 (Supplementary Figure S3D). Studies have demonstrated that LGG patients with higher mRNAsi indices have a better prognosis than those with low mRNAsi indices (Tan et al., 2021) (Supplementary Figure S3E). We obtained and analyzed mRNAsi differences in LGG patients from our data, and found that the mRNAsi indices of MC2 and MC3 were significantly higher than that of MC1. We also analyzed epigenetic regulation based on index (EREGmRNAsi) differences among different metabolic subtypes, and found that there was no significant difference between EREGmRNAsi subtypes (Supplementary Figure S3F).

To compare the differences between the three metabolic subtypes identified by us and the seven metabolic subtypes reported in previous studies (Ceccarelli et al., 2016), we determined that Codel and G-CIMP-high were higher in MC2 and MC3 subtypes. The percentages of patients with Codel and G-CIMP-high were 35 and 57%, respectively, in the MC2 subtype, and 55 and 36%, respectively, in the MC3 subtype. Patients with Codel and G-CIMP-high had a better prognosis (Figure 3A). We also compared the metabolic subtypes with the six previously reported immunomolecular subtypes (Thorsson et al., 2019), and found C3 to have the best prognosis and C4 and C6 subtypes to have the worst. The percentages of C5 subtype in MC2 and MC3 subtypes were 77 and 86% respectively, while the percentage of C4 subtype in the MC1 subtype was 62%, higher than those in MC2 and MC3 subtypes (Figure 3B). Immuno-infiltration analysis showed that MC1 had the highest immuno-microenvironment infiltration in the TCGA cohort, followed by MC2, and MC3 had the lowest immune-infiltration score (Figure 3C). Consistent findings were also observed in the CGGA cohort: MC1 subtype had the highest immune infiltration, significantly higher than that of MC2 and MC3 metabolic subtypes (Figure 3D). We then applied the CIBERSORT method to investigate the components of immune cells in each metabolic subtype. The results demonstrated that there were significant differences among different subtypes of immune cells. In the TCGA cohort, the Macrophages_M2 was significantly enriched in the MC1 subtype, and the CD8⁺ T cells, Naïve CD4⁺ T cells, activate CD4⁺ memory T cells, follicular helper T cells, Treg, gamma delta T cells, M1 macrophages, activated dendritic cells, mast cells, eosinophils, and neutrophils were differently expressed in the three clusters (Figure 3E). In the CGGA cohort, Macrophages_M2 in the MC1 subtype was significantly higher than that in MC2 and MC3, while naïve B cells, memory B cells, plasma cells, CD8⁺ T cells, resting memory CD4⁺ T cells, gamma delta T cells, resting NK cells, monocytes, M1 macrophages, and neutrophils were differently expressed in these classifications (Figure 3F). In addition, we also applied the xCell, EPIC, and MCP-counter to comprehensively evaluate the tumor-infiltrating immune cell score in the TCGA database and CGGA dataset. And results demonstrated that the MC1 subtype had the highest immune infiltration, while the MC2 and MC3 metabolic subtypes were significantly lower than the MC1 subtype. In addition, most immune cells had significant differences between subtypes in the TCGA database and CGGA database Supplementary Figure S4. These findings suggested that our metabolism-related classification was closely related to immune cell infiltration, which may be an effective evaluation method in LGG immune evaluation.

To evaluate the different immunotherapy and potential clinical effects of metabolism-associated subtypes, we adopted TIDE (http://tide.dfci.harvard.edu/) software and determined that MC had the highest TIDE score in the TCGA cohort, which suggested a higher possibility of MC2 immune escape and a lower possibility of benefit from immunotherapy (Figure 4A). We also compared the differences in T cell dysfunction scores and T cell rejection scores among different molecular metabolic subtypes in the TCGA cohort. The results indicated that the MC2 subtype had the lowest T cell dysfunction scores but the highest T cell rejection scores. The T cell dysfunction scores of MC1 and MC3 subtypes were higher than that of the MC2 subtype, but the T cell rejection scores were lower. However, there was no significant difference in T cell dysfunction scores and T cell rejection scores between the MC1 and MC3 subtypes (Figures 4B,C). Similar results were also observed in the CGGA cohort (Figures 4E–G). In addition, we investigated the differences in the predicted immunotherapy response for different metabolic subtypes, and the results showed significant differences in immunotherapy response status between the MC2 and MC3 subtypes. In the CGGA cohort, there was a significant difference between the MC1 subtype and the MC2 and MC3 subtypes (Figures 4D,H). Furthermore, we carried out subclass mapping and compared the immunotherapy data on the three metabolic subtypes; we found that the MC1 subtype in the TCGA and CGGA cohorts was sensitive to anti-PD1 therapy, while the MC3 subtype was sensitive to CTLA4-R (Figures 4I,J).

Because different metabolic subtypes have different metabolic characteristics, we performed PCA to construct a subtype classification feature score and thus better quantify the metabolism-related characteristics of each sample. The PCA indicated that PC1 and PC2 could successfully discriminate and classify LGG samples according to the metabolic subtypes classification features in the TCGA database (Figure 5A). The results revealed that the three metabolic subtypes had significantly different microenvironment-related genes scores (Figure 5B). The receiver-operating characteristic (ROC) curve, combining sensitivity and specificity, showed the performance of different metabolic subtype feature scores of different clusters. The multiclass area under the curve (AUC) was 0.86, which indicated that this signature score model had excellent predictive power (Figure 5C). We also observed similar features in the CGGA database (Figures 5D–F), in which the multiclass AUC was 0.85.

To assess the correlation between metabolic subtype characteristics scores and immune cell characteristics, we performed Pearson’s correlation analysis, and found that metabolic subtypes signature scores had no significant correlation with immune cells, with the exception of M2 macrophages and eosinophils.

These results highlighted the fact that M2 macrophages had a remarkably strong negative correlation with metabolism-associated scores, and eosinophils had a positive correlation with these scores (Supplementary Figure S5A). The results for M2-macrophage infiltration in the CGGA cohort agreed with those from the TCGA cohort (Supplementary Figure S5B). Furthermore, to better describe the correlation analyses between metabolic subtypes, metabolic score, and immune checkpoint molecules, we adopted correlation row analysis, and we found that the metabolic subtypes, metabolic score and immune checkpoint molecules have significant correlation with each other (Supplementary Figure S6). These findings suggested that the metabolic signature had good cooperativity for not only metabolic score but also immune associated molecules.

To further identify metabolism subtypes associated with co-expressed gene models, we performed a weighted gene co-expression network analysis (WGCNA) to identify modules connected with a variety of LGG metabolism signature-based subtypes. The genes, whose MAD (median absolute difference) were more than 50%, were selected for further WGCNA analysis in the TCGA gene expression profile. The sample clusters are illustrated in Figure 6A. To ensure a scale-free network, a power of β = 8 (scale-free R2) was selected as the soft-thresholding parameter (Figures 6B,C). Furthermore, similar clusters were merged into new modules using the following settings: height = 0.25, deepSplit = 2, and min ModuleSize = 30. This produced 18 gene molecules (Figure 6D). The 18 gene numbers of each co-expression network are illustrated in Figure 6E. We further analyzed the correlation between each co-expressed gene module and Age, Gender, IDH1/IDH2 mutation status, MGMT methylation status, 1P/19Q CODEL status, and MC1, MC2, and MC3. The results showed that the green module had the highest correlation with MC1, the Cyan module was significantly positively correlated with MC2, and the blue module was significantly positively correlated with MC3 (Figure 6F). The Green module had a significant correlation with MC1 (Figure 6G, cor = 0.69, p = 2.4e−82). The Cyan module had a significant correlation with MC2 (Figure 6H, cor = 0.47, p = 7.8e−11). The Blue module had a significant correlation with MC3 (Figure 6I, cor = 0.62, p = 1.2e−116).

We calculated the correlation between the feature vectors of the 17 modules (excluding the Grey module) and the metabolic subtype feature index, as shown in Figure 7A. The Green and Cyan modules had significant negative correlation with the metabolic subtype feature score (Figures 7B,C). The blue module was significantly positively correlated with the metabolic subtype characteristic score (Figure 7D). We screened the genes of the three modules for functional enrichment, with the following results. As shown in Figures 7E–G, our green module is related to leukocyte proliferation, positive regulation of leukocyte activation, T cell activation, and other immune processes. The KEGG pathway was enriched in the intestinal immune network for IgA production and Th17 cell differentiation. The Cyan module was significantly enriched to protein targeting to membrane, protein localization to endoplasmic reticulum, and other processes (Figure 7F). Furthermore, the results of significant enrichment in the blue module are shown in Figure 7G.

In addition, we identified 24, 13, and 21 key genes in the green, Cyan, and Blue modules, respectively, with correlation coefficients greater than 0.9 and significant association with prognosis (Figure 8A). We then used a Venn diagram to identify the intersection of a total of eight key genes: BTK, ALOX5, ARPC1B, REL-10, ITGB2, GPSM3, LAPTM5, and MYO1F (Figure 8B).

According to gene expression levels, patients were divided into a high-expression group and low-expression group. Kaplan-Meier prognostic analysis revealed that high expression levels of these eight genes were correlated with poor prognosis (Figure 8C). In conclusion, these eight genes are potential markers associated with the metabolic subtype characteristics score.