Identification and Characterization of Cancer-Related Risk Metabolic Subpathways Reveal Their Functional Significance in Cancer

Article type: Research Article

Keywords: subpathways, biomarkers, metabolism

Authors: Hongying Zhao, Jinxing Yan, Ming Wu, Shiyi Li, Weiming He, Xiangzhe Yin, Wangyang Liu, Ying Liu, Meiting Fei, Wan Li ⚬, Junjie Lv, Lina Chen, Li Wang ⚬

Affiliations: College of Bioinformatics Science and Technology, Harbin Medical University, Harbin 150081, China; zhaohongying@hrbmu.edu.cn (H.Z.); yanjinxing0712@163.com (J.Y.); wm734393188@163.com (M.W.); lsyshea@163.com (S.L.); yxzswyxgc@163.com (X.Y.); liuwangyang2022@163.com (W.L.); liuying2022020510@163.com (Y.L.); fmt2844143432@163.com (M.F.);; Institute of Opto-Electronics, Harbin Institute of Technology, Harbin 150000, China; hewm@hit.edu.cn

License: © 2026 by the authors. CC BY 4.0 Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license.

Article links: DOI: 10.3390/ijms27104246  |  PMC: PMC13207224

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (7.7 MB)

1. Introduction

The reprogramming of metabolism, the process of tumor cells reprogramming metabolic pathways to meet the energy and biosynthetic demands of tumor cells in order to adapt to hypoxia and nutrient deficiency, is a hallmark of malignancy with great potential for prognosis and targeted therapy [ref. 1]. For example, the abnormal activation and reprogramming of multiple key metabolic pathways, such as glycolysis, fatty acid synthesis, and amino acid metabolism, are considered to be one of the key factors promoting tumor progression in breast cancer [ref. 2]. Recently, an increasing number of studies have reported that abnormalities in mRNA, proteins, and metabolites play important roles in metabolic reprogramming in cancer. For instance, the study by Torresano et al. identified the up-regulation of the glucose transporter GLUT1 in cancer as a mechanism to increase glucose uptake to support tumor growth. Consistently, the overexpression of GLUT1 is associated with increased tumor aggressiveness and poor survival in patients with colon, hepatocellular, esophageal squamous cell carcinoma, and non-small lung cancer [ref. 3]. The up-regulation of the IGF2BP3 mRNA and protein expression regulates nicotinamide metabolism, which could alter oxidative phosphorylation and promote acquired resistance to EGFR inhibitors [ref. 4]. In addition, cancer metabolic reprogramming not only promotes tumor growth and metastasis but is also closely related to the immune response in the tumor microenvironment, which affects the progression of tumors and the prognosis of patients. Studies have shown that abnormal activation of glucose metabolism promotes tumor cell survival and proliferation through the Warburg effect and activates M2-like macrophages through lactate production, leading to tumor immunosuppression [ref. 5,ref. 6,ref. 7]. Therefore, understanding how reprogrammed metabolism supports the growth and evolution of specific tumor cells and identifying which reprogrammed metabolic activities are most relevant to therapeutic vulnerabilities can help reveal the mechanism of cancer progression.

Traditional pathway identification methods, using databases such as KEGG, Reactome, Pathway Commons, and WikiPathways, rely on abstract concepts built on whole biological concepts, such as apoptosis or the cell cycle. Moreover, at the genomic and transcriptomic levels, no unique mechanistic results are provided for the activity of the whole pathway, and the experimental results are indicative and of limited utility. In the context of cancer, it is not a universal dysfunction of all genes within a pathway; rather, it is typically a subset of genes that exhibits aberrant functions within specific segments of the pathway. These specific parts, known as subpathways, which are local gene subregions within biological pathways, have drawn significant attention. Metabolic subpathway identification focuses on specific branches or submodules in a pathway, which can more accurately locate metabolic processes related to specific biological functions, and decompose the complex metabolic network into simpler and clearer subunits [ref. 8]. Abnormalities in metabolic subpathways have been reported to be associated with the occurrence, progression, and prognostic outcomes of cancer [ref. 9]. For example, (S)-2,3-Epoxysqualene, located in the sterol biosynthesis subpathway, which holds a central and crucial position, is likely to have a more significant impact on the entire pathway and thus potentially has a high association with glioblastoma [ref. 10]. The identification of cancer-related subpathways not only benefits the understanding of carcinogenesis and cancer progression, but also helps to identify more precise biomarkers for cancer patients. Cancers were highly heterogeneous diseases that had great differences in transcriptome, survival prognosis and other aspects. How to effectively integrate metabolic pathways and multiple omics data to accurately identify abnormal metabolic subpathways in cancer is a challenge. In this paper, we proposed a computational approach to identify cancer-related risk metabolic subpathways based on the topology of large-scale metabolic networks and cancer-specific multi-omics. Moreover, we applied this algorithm to identify risk metabolic subpathways in multiple cancer types. We comprehensively characterized cancer-related metabolic subpathways and exhibited expression perturbations in cancer, which were significantly correlated with immune cell infiltration and prognosis. Thereby, the metabolic subpathways identified by our method are crucial for cancer, providing novel insights for immunotherapy, drug response, and prognosis assessment.

2. Results

2.1. Identification of Cancer-Related Risk Metabolic Subpathways Across Cancer Types

To systematically identify subpathway regions of abnormal metabolic pathways that may contribute to complex human diseases, we developed a novel computational approach, CMSubpathway, to identify cancer-related risk metabolic subpathways (Figure 1A). First, 84 metabolism-related pathways were downloaded from the KEGG database and then transformed into interaction networks. The metabolic subpathways were identified based on the LTE algorithm, which is a fast local expansion algorithm for discovering communities in large-scale networks. Second, genes that were differentially expressed between cancer samples and normal samples (|logFC| ≥ 1 and p < 0.05) were identified, and standardized enrichment scores (NESs) for each metabolic subpathway were calculated based on gene set enrichment analysis (GSEA). The subpathways with an absolute NES > 1, p < 0.05 and FDR < 0.25 were considered dysregulated metabolic subpathways in cancer. Third, candidate subpathways were identified as prognostic markers for each cancer. The GSVA score for each subpathway in every sample was calculated, and the samples were categorized into high-risk and low-risk groups based on their GSVA scores. The log-rank test was applied to evaluate the differences between the two groups, with p < 0.05 considered significant. Fourth, the classification performance of candidate subpathways was evaluated. An SVM model was constructed to assess subpathway classification efficiency by using the GSVA score. An AUC > 0.75 was considered efficient. Finally, we defined cancer-related metabolic subpathways that are dysregulated in cancer, correlate with prognosis, and facilitate the classification of normal and cancer samples.

Taking advantage of the expression and clinical data in The Cancer Genome Atlas (TCGA) and metabolic pathways in KEGG, the CMSubpathway was carried out to extensively identify cancer-related risk metabolic subpathways across 21 different cancer types. We identified 12 cancer-related metabolic subpathways across six types of cancer: breast cancer (BRCA), kidney renal clear cell carcinoma (KIRC), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), pheochromocytomas and paragangliomas (PCPG), and thyroid carcinoma (THCA) (Figure 1B). Specifically, we identified three BRCA-related risk metabolic subpathways, including hsa00620_2 (pyruvate metabolism pathway), hsa00330_4 (Arginine and proline metabolism pathway), and hsa00071_2 (fatty acid degradation; Figure 1C, Table S1). The subpathway hsa00620_2 was shared by BRCA and KIRC, while hsa00071_2 was shared by BRCA and THCA. For example, the low activities of the risk metabolic subpathways hsa00620_2, hsa00071_2, and hsa00330_4 were associated with a poor prognosis in BRCA (for hsa00620_2, p = 4.9 × 10⁻²; for hsa00071_2, p = 4.7 × 10⁻³; for hsa00330_4, p < 1.0 × 10⁻⁴; Figure S1A). It is consistent with our findings that the inhibition of pyruvate metabolism was sufficient to impair collagen hydroxylation and consequently the growth of breast cancer-derived lung metastases in different mouse models [ref. 11]. Decreased activity of fatty acid metabolism genes was observed in RCC tumors and was described as a feature of ccRCC [ref. 12]. In conclusion, CMSubpathway is capable of effectively identifying cancer-related risk metabolic subpathways and uncovering the commonalities and specificities of cancer metabolic reprogramming, thereby offering novel insights for further clinical research.

2.2. Independent Dataset Validation

Since breast cancer has the largest number of risk metabolic subpathways and shares a relatively large number of subpathways with other cancers, we take breast cancer as an example to verify the effectiveness of the method. An independent validation dataset, which contains matched RNA-seq and survival data, was obtained from the GSE42568 dataset, consisting of 104 breast tumor samples and 17 normal samples. The same cancer-related risk metabolic subpathway identification method was applied to the validation dataset. In line with the outcomes from the TCGA-BRCA dataset, the down-regulation of two out of the three risk metabolic subpathways (hsa00071-2 and hsa00620-2) within tumor samples exhibited a significant correlation with survival prognosis (log-rank p = 0.017 and 0.016, Figure 2A, Table S2). Elevated expression levels of these subpathways were associated with a decreased risk of adverse outcomes. Furthermore, these subpathways exhibited fine classification performance in distinguishing between normal and diseased samples (AUC = 0.76 and 0.72, Figure 2A, Table S2). We regarded the subpathways hsa00071-2 and hsa00620-2 that were consistent between the two datasets as the breast cancer risk metabolic subpathways. Moreover, to enhance the confidence of our method, we utilized the Subpathway-CorSP algorithm [ref. 13] instead of the LTE algorithm in our approach to validate the identified metabolic subpathways within the TCGA-BRCA and GSE42568 datasets. The TCGA-BRCA dataset and GSE42568 dataset consistently identified one risk metabolic subpathway (hsa00830_1) that was significantly down-regulated and associated with a high risk due to low expression (p < 0.05, Tables S2 and S3). The subpathway hsa00830_1 was significantly down-regulated in patients and showed significant prognostic value (log-rank p = 1.5 × 10⁻²) in both the TCGA-BRCA and GSE42568 datasets, effectively distinguishing between normal and cancer samples (AUC = 0.858 in TCGA and AUC = 0.781 in GSE42568; Figure 2B,C and Figure S1B, Tables S3 and S4). The subpathway hsa00830_1 has a significant intersection of metabolic genes with breast cancer risk subpathways hsa00071-2 and hsa00620-2 (hypergeometric test, p = 3.17 × 10⁻¹⁰ for hsa00071-2 and p = 4.04 × 10⁻⁸ for hsa00620-2) by our method using the LTE algorithm, including shared genes ADH1A, ADH1B, ADH1C, ADH4, ADH5, ADH6, and ADH7 (Figure S1C, Tables S2 and S3). Furthermore, to assess the performance of CMSubpathway, we performed a comparison with several established subpathway analysis methods, including Subpathway-CorSP, LTE, and SubpathwayGMir. Predictive accuracy and robustness were assessed across the aforementioned 11 independent validation datasets. The results demonstrated that our methodology achieved substantially higher predictive accuracy (mean AUC = 0.95) and greater algorithmic robustness than LTE (mean AUC = 0.917), Subpathway-CorSP (mean AUC = 0.924) and SubpathwayGMir (mean AUC = 0.898) across all tested datasets (Figure S1D). These results further validated the feasibility and reliability of our proposed method.

2.3. Metabolic Crosstalk of the Risk Metabolic Subpathway Driven by the Core Module in BRCA

2.3.1. Consistent Dysregulation of the Core Metabolic Module at Multiple Levels

Based on the significant intersection of metabolic genes in the breast cancer risk metabolic subpathways hsa00071-2 and hsa00620-2 identified above (p = 3.12 × 10⁻²³), we identified a set of 12 overlapping genes. Biologically, these genes (predominantly encoding ADH and ALDH family enzymes) reside at the critical topological crossroads linking fatty acid degradation and pyruvate metabolism, facilitating essential metabolic crosstalk. Given their key roles in bridging these distinct metabolic processes, we defined this 12-gene set as a breast cancer-related core metabolic module (Figure S1E). We found that the core metabolic module genes exhibit consistent down-regulation at the transcriptome and proteome levels in breast cancer patients. Specifically, at the transcriptomic level, through analyzing the results of differentially expressed genes from TCGA, GSE42568 and the gene expression data of breast cancer from the cancer atlas of metabolic profiles (CAMP), the core metabolic module genes consistently showed significant down-regulation results in tumor patients (Figure S2A,B). Additionally, we collected two independent breast cancer clinical datasets from the Expression Atlas database. The first dataset (E-GEOD-45581) comprises 40 untreated breast cancer patient samples and 5 normal controls. The second dataset (E-MTAB-779) includes 20 untreated breast carcinoma patient samples and 22 normal controls. These two breast cancer clinical datasets were used to validate the expression differences in core metabolic module genes. The results showed that most core metabolic module genes were down-regulated in breast cancer tissues relative to normal tissues (Figure S2C). Furthermore, we also investigated the expression of the genes associated with this core module within the peripheral blood cohort of breast cancer patients and discovered that they were down-regulated across multiple sets of blood samples. For example, ALDH3A2, ALDH2, ALDH9A1, ADH5, and ALDH7A1 had significantly lower expression in multiple cancer blood sample cohorts (Figure S2D, Table S4). At the proteomic level, we obtained the protein expression data from the CPTAC dataset, which includes 18 normal and 218 BRCA samples, and conducted a differential analysis. Consistent with the RNA-seq results, these two risk metabolic subpathways and their core metabolic module were significantly down-regulated at the protein level in breast cancer samples (Figure 3A,B), which revealed a congruent down-regulation, corroborating the RNA-seq level findings. These 12 key genes also exhibited a correlation pattern at the protein level that mirrored the RNA-seq data, with ADH1A, ADH1B, and ADH1C demonstrating pronounced correlations in both tumor and normal samples (Figure 3C). Additionally, the protein expression levels of the core metabolic module genes were significantly down-regulated in tumor tissue samples, which is consistent with the findings at the RNA-seq level (Figure 3D). To further validate our results, we searched the HPA database for the IHC data of the proteins encoded by these 12 key genes bridging the risk subpathways in both normal and tumor tissue sections. Through the analysis using the HPA database, it was found that, among the 12 key genes, the protein expression levels of 9 were decreased in breast cancer tissues relative to normal tissues (Figure 3E). In conclusion, we observed the dysregulation of genes within this core metabolic module at multiple levels, which might contribute to the metabolic reprogramming in breast cancer.

2.3.2. Dysregulation of the Core Module Is Associated with Metabolic Reprogramming

These core metabolic module genes jointly regulate pyruvate metabolism, glycolysis/gluconeogenesis, and fatty acid degradation, affecting metabolic crosstalk and reprogramming (Figure 4A and Figures S3–S5). The genes within this module can be divided into ALDH and ADH families according to their regulatory roles in metabolism, which are derived from the aldehyde dehydrogenase and alcohol dehydrogenase gene families and exhibit strong internal correlation (Figure 4B and Figure S6A). To further investigate whether these key genes play a central role in subpathway crosstalk, we carried out a functional enrichment analysis of the differentially expressed gene sets in breast cancer, with and without their being considered separately. The results of our analysis demonstrated that only the gene sets that included this core metabolic module were significantly enriched in the pyruvate metabolism pathway (p = 2.63 × 10⁻⁴) and the fatty acid degradation pathway (p = 1.33 × 10⁻⁴), whereas the gene sets without these key genes did not show such significant enrichment (Figure S6B). By utilizing the gene expression and metabolites data of breast cancer, which encompassed 61 cancer samples and 47 normal samples from the CAMP [ref. 14], we found that the expression levels of ADH1A, ADH1B, and ADH1C were significantly down-regulated and positively correlated with glucose levels in cancer samples (Pearson correlation test, with p < 0.05 and Pearson correlation coefficient > 0). However, such a correlation was not detected in normal samples (Figure 4C and Figure S2B). This finding implies that the ADH genes within this module might be related to the high glucose consumption trait in breast cancer. Glucose serves as the principal source of lactate in the body. Glycolysis triggers the production of lactate, and lactate acts as the main fuel for the tricarboxylic acid cycle (TCA cycle), supplying energy for tumor cell respiration [ref. 15,ref. 16]. We hypothesized that the down-regulation of the ADH and ALDH genes of the core metabolic module inhibits lactate metabolism, thereby maintaining a state of acidosis, while ensuring that lactate, with the catalysis by the up-regulation of ACSS1, generates acetyl-CoA to participate in the TCA cycle (Figure 4B). Indeed, correlation analysis indicated a significant negative correlation between the expression of ADH1A, ADH1B, ADH1C genes and lactate levels in breast cancer (Figure 4D). At the same time, the differential results of metabolite levels showed that the lactate levels in cancer samples were significantly higher than those in normal samples (p = 4.86 × 10⁻⁸; Figure 4E). The lactate accumulates as a reducing product of glucose in the tumor microenvironment during glycolysis based on the Warburg effect [ref. 17]. Lactate can serve as a metabolic substrate and signaling molecule, playing a predominant role in metabolic crosstalk between cancer cells and the tumor microenvironment (TME) [ref. 18,ref. 19]. Additionally, the core metabolic module genes were involved in the processes in which alcohol can also be catalytically converted into aldehyde and in the metabolism of fatty acids. In cancer progression, fatty acid degradation, along with lactate metabolism, represents an effective energy acquisition pathway. The differentially expressed ADH impacts the oxidation of 1-alcohol to produce aldehyde and, under the dysregulated ALDH regulation, participates in the synthesis of fatty acids, thereby exerting significant influences on cancer metabolism (Figure 4B). The combined action of ADH and ALDH, through their effects on aldehyde production and fatty acid synthesis, significantly impacts cancer metabolism, potentially altering the energy dynamics and metabolic reprogramming that are characteristic of cancer cells. Indeed, a total of 17 RNA samples were isolated from MCF7 xenograft tumors, which were co-implanted with either UCP1-CRISRPa-modulated or dCas9-VP64 (control)-treated adipose organoids in SCID mice fed with chow, high-fat diet, or 15% glucose water (GSE246231). These samples were obtained to explore the relationship between fat interference and the expression of this core metabolic module. The differential analysis revealed that ALDH2 (p = 0.041), ALDH9A1 (p = 0.045), and ADH5 (p = 0.041) were significantly lower expressed in high-fat diet samples than in chow diet samples (Figure S6C). To further support our hypothesis about the essentiality of these core metabolic module, we utilized publicly available CRISPR-Cas9 knockout screening data from Depmap (Dependency Map; https://depmap.org/portal/, accessed on 8 April 2026). The results revealed that, upon CRISPR-mediated knockout of ADH5 in breast cancer cell lines, the CERESs were <0 across all evaluated cell lines. Similarly, following the knockout of ALDH1B1 and ALDH7A1, over 75% of the cell lines exhibited CERESs < 0 (Figure S7A). These findings substantiate that ADH5, ALDH1B1, and ALDH7A1 play critical roles in the proliferation and survival of breast cancer cells. Additionally, ADH plays an essential role in drug metabolism involving cytochrome P450 (CYPs), such as drug metabolism-cytochrome P450 and the metabolism of xenobiotics by cytochrome P450 (Figure S5B,C). Research has elucidated the roles of ADH and CYP in ethanol pharmacokinetics through pharmacokinetic analysis [ref. 20]. Thus, the dysregulated core metabolic module genes, serving as crucial regulatory enzymes in multiple metabolic pathways, might disrupt the metabolism subpathway for glucose and fatty acids, consequently being closely associated with the energy reprogramming related to the tumor cell growth of breast cancer.

2.3.3. Functional Characteristics of the Core Module Linking Risk Metabolic Subpathway

Considering the crucial role of the core metabolic module in energy and drug metabolism, we will further investigate their regulation on the progression of breast cancer. By using the GSVA algorithm to calculate the functional scores of the core metabolic module genes and dividing patients into high and low metabolic activity groups, based on these scores, significant differences in survival could also be observed between these two groups. Given that patients may receive different treatment regimens, we evaluated the prognostic value of the metabolic signature within the major treatment subgroups of the TCGA-BRCA cohort, in which 49.7% of patients were chemotherapy-naive and 22.9% received anthracycline–taxane combination regimen (AC-T) chemotherapy. For each subgroup, we stratified patients into high- and low-score groups based on the key subpathway score and performed survival analyses. As shown in Figure S7B, the metabolic signature remained prognostic in both the AC-T subgroup (n = 248, p =3.6 × 10⁻²) and the treatment-naive subgroup (n = 554, p = 3.9 × 10⁻³). Furthermore, we analyzed the differential genes between two groups, including 258 down-regulated genes and 3 up-regulated genes. The under-expressed genes in the high metabolic activity group exhibited an overlap with the key genes (ADH1A, ADH1B, ADH1C, and ADH4) of ADH in breast cancer (Figure S8A). Studies have found that ADH1B may inhibit the proliferation, invasion and migration of BRCA cells by inactivating the MAPK signaling pathway, making it a promising target for the clinical treatment of BRCA [ref. 21]. The association between ADH1C polymorphisms and the risk of breast cancer was significant in a meta-analysis of a case–control study [ref. 22]. The functional pathways enriched in the high metabolic activity group include pyruvate metabolism (p = 9.42 × 10⁻⁴), fatty acid degradation (p = 6.23 × 10⁻⁴), glycolysis/gluconeogenesis (p = 7.38 × 10⁻⁴), drug metabolism-cytochrome P450 (p = 1.16 × 10⁻³), and the AMPK signaling pathway (p < 0.05) (Figure S8B). The AMPK signaling pathway, involved in cellular energy metabolism, impacts breast cancer progression by enhancing glucose uptake and aerobic glycolysis, leading to lactate accumulation [ref. 23]. These findings are highly consistent with the regulation of metabolic pathways by the core metabolic module, further highlighting their important role in regulating metabolism in BRCA. In terms of biological processes, they are related to alcohol metabolism and immune-related functions. Therefore, in-depth exploration of this core metabolic module in the metabolic mechanisms of breast cancer holds significant clinical significance and promising potential applications.

2.3.4. Subpathway Crosstalk Mediated by the Core Module Correlates with Immune Cell Infiltration in BRCA

Given the enrichment of immune-related functions in the high metabolic activity group, we further investigated the impact of the core metabolic module on immune cell infiltration. The results showed that the immune infiltration levels of multiple immune cells, including B cells naive, CD8 T cells, resting memory CD4 T cells, and monocytes, are significantly higher in the high metabolic activity group, whereas the levels of M0 macrophages, M2 macrophages, and resting mast cells are significantly lower than those in the low metabolic activity group (Figure 4F). In the GSEA enrichment results of the c7 immune-related gene set, the high metabolic activity group was significantly enriched in multiple T cell-related immune gene sets, such as GSE19198_TREATED_TCELL_UP, GSE3039_CD8_TCELL_DN and GSE36476_MEMORY_CD4_TCELL, indicating that the core metabolic module may play a crucial role in T cell immune infiltration (p < 0.01; Figure S8C,D). More and more evidence suggests that the immune response of tumor cells is influenced by the acidification of the tumor microenvironment (TME) caused by lactate accumulation. It is well-known that TME acidification is an immunosuppressive factor that hinders tumor immune responses [ref. 24]. Activated T cells are primarily characterized by glycolysis, which involves glucose consumption and lactate production, further contributing to lactate accumulation [ref. 25]. Therefore, we have reason to believe that the dysregulation of the core metabolic module is potentially linked to lactate accumulation and the acidification of TME, which may contribute to an immunosuppressive environment and consequently correlate with altered immune cell infiltration in breast cancer.

In conclusion, by studying the regulatory mechanisms of this core metabolic module, we would gain a deeper understanding of their roles in metabolic crosstalk and reprogramming, and further elucidate their impact on breast cancer progression and immune cell infiltration. The study of this core metabolic module in breast cancer holds great promise, providing new insights and approaches for the diagnosis and treatment of breast cancer.

2.3.5. Metabolic Reprogramming of the Core Module Linking Risk Subpathways Associated with the CD8 T Cell Differentiation and Exhaustion

To explore the mechanism of action of the core metabolic module in the breast cancer tumor microenvironment, we collected a BRCA scRNA-seq cohort (GSE195861) [ref. 26]. After quality control, 51,002 cells were contained and the cells were divided into five major cell types, including T cells (n = 12,342) identified by the expression of CD3D and CD3E, B cells (n = 11,175) which expressed MS4A1, CD79A and CD79B, fibroblasts (n = 97) marked by PDGFRB, COL1A1 and COL1A2, myeloid cells (n = 3331) identified by CD68, CD14, FCGR3A and LYZ expression, and epithelial cells (n = 24,057) which were positive for EPCAM, CDH1, KRT8 and KRT18 expression (Figure 5A,B). Next, to confirm the malignant cell subtypes, epithelial cells were further reclustered into 14 subtypes, including 7 malignant cell subtypes (0, 1, 4, 6, 7, 9, 12; n = 15,515) and 7 normal epithelial cell subtypes (2, 3, 5, 8, 10, 11, 13; n = 8542) using CNV analysis with immune cells as references [ref. 27] (Figure S9). As metabolic pathways are dysregulated in tumor cells as novel cancer markers, we used SiPSiC to establish core metabolic module scores in the breast cancer tumor microenvironment. First, analysis of the difference in these core metabolic module scores between malignant and normal cells showed that the core metabolic module scores were significantly lower in malignant cells (p < 2.2 × 10⁻¹⁶; Figure 5C). This suggests that this core metabolic module is less active in malignant cells, which is consistent with findings in the TCGA-BRCA and GSE42568 cohorts (TCGA-BRCA: p < 2.2 × 10⁻¹⁶; GSE42568: p < 3.6 × 10⁻⁹; Figure 5D). Moreover, this core metabolic module score using ssGSEA was significantly negatively correlated with tumor cells in the TCGA-BRCA, GSE42568 and five breast cancer peripheral blood cohorts (GSE67939, GSE51827, GSE111842, GSE75367 and GSE86978; R < 0, p < 0.05; Figure S10; Table S5). Next, investigating the genes within this core metabolic module is very important for us to further understand the role in mediating subpathway crosstalk. ALDH3A2, ALDH7A1, ALDH9A1, ADH5 and ALDH2 exhibited a significantly lower expression level in malignant cells (p < 2 × 10⁻¹⁶), and consistent results were also observed in cancer samples from the TCGA-BRCA and GSE42568 cohorts (p < 0.05; Figure 5E). Results illustrated that this core metabolic module was significantly dysregulated and its related genes had significantly lower expression levels in cancer cells and samples through the integration of single-cell and bulk data.

Since the core metabolic module was correlated with the immune cell infiltration of CD8 T cells, CD8 T cells were subclassified through clustering analysis into eight subtypes (Figure 6A). To reveal the origin of CD8 T cells at the single-cell level, CytoTRACE analysis was performed to decipher the differentiation state of CD8 T cell subtypes. A higher CytoTRACE score was predicted in CD8 T cell subtype 0 than in other subtypes (Figure 6A,B). The T cell subtype 0 had the highest cell stemness, representing the lowest degree of differentiation. Notably, we found that CD8 T cell subtype 0 exhibited a higher expression of exhausted T cell marker genes (PDCD1, HAVCR2, LAG3, and TOX) and a lower expression of transcription factor TCF7 (Figure S11A), which was consistent with stem-like CD8 T cells in a previous study [ref. 28]. Moreover, we observed a significant positive correlation between core metabolic module scores and CytoTRACE score in CD8 T cells (R = 0.39; p < 2.2 × 10⁻¹⁶; Figure 6C). Additionally, CD8 T cell subtype 0 with high cell stemness exhibited the highest core metabolic module scores (p < 2.0 × 10⁻¹⁶; Figure 6D and Figure S11B). Especially, we also observed a significant positive correlation between the core metabolic module scores and CytoTRACE scores in CD8 T cell subtype 0 (R = 0.23; p < 3.1 × 10⁻¹⁴; Figure 6E). Using ssGSEA, the core metabolic module was found to a significant positive correlation with CD8 T cell subtype 0 in the TCGA-BRCA cohort (R = 0.24; p < 3.1 × 10⁻¹⁵; Figure 6F). These results suggest that the core metabolic module in the tumor microenvironment may be involved in the differentiation process of CD8 T cells and affect the activation and exhaustion of CD8 T cells, which is consistent with previous studies [ref. 29,ref. 30]. Furthermore, the infiltration score of CD8 T cell subtype 0 was significantly negatively correlated with tumor cells in TCGA-BRCA (R = −0.097; p = 0.0016; Figure S11C) and GSE42568 (R = −0.22; p = 0.023) cohorts (Figure S11D). This is consistent with previous studies that stem-like CD8 T cells are critical for maintaining anti-tumor response and immunotherapy response [ref. 28,ref. 31,ref. 32]. Taken together, these results suggested that metabolic reprogramming of the core metabolic module is highly correlated with the differentiation and exhaustion of CD8 T cells. This association implies a potential impact on anti-tumor immunity, highlighting the core metabolic module as a promising therapeutic target for BRCA immunotherapy.

3. Discussion

Subpathways, as the crucial components of pathways, had the potential to be biomarkers for clinical applications. In this study, we developed a novel computational approach CMSubpathway to identify cancer-related risk metabolic subpathways that were dysregulated in cancers and possessed the capabilities of cancer prognosis and classification. Firstly, differential expression analysis showed that 33.40% of the 156 metabolic subpathways were dysregulated in cancer (p < 0.05, Figure S12A), and 66.04% of the dysregulated subpathways were present in multiple cancer types. For instance, the subpathways hsa00620_2, hsa00071_2, and hsa00982_1 (drug metabolism-cytochrome P450), were dysregulated in more than five cancer types (Figure S12B). Further analysis revealed that a minority of dysregulated subpathways were associated with prognostic and classification efficacy, particularly in BRCA and PCPG (Figure S12C). Taking breast cancer as an example, an independent validation dataset was used to verify the robustness of the method. The discovery set and the validation set consistently identified two risk subpathways, hsa00071-2 and hsa00620-2, as breast cancer risk metabolic subpathways. Furthermore, another subpathway identification method was used to verify the robustness of the method. Consistently, one risk subpathway was identified in both the discovery set and the validation set, which had a significant intersection with the breast cancer risk subpathways. These results validated the feasibility and reliability of CMSubpathway to identify cancer-related risk metabolic subpathways.

Furthermore, we found a significant gene intersection between the two risk metabolic subpathways of breast cancer, so we regarded the intersection genes as a core metabolic module linking these subpathways. Using multiple sets of expression data and protein datasets, the genes of this core metabolic module showed consistent down-regulation at both the transcriptome and proteome levels in breast cancer patients. The core metabolic module gene ADH1A exhibited a significant correlation with breast cancer prognosis (Figure S11D; log-rank p = 0.032 for disease-free survival and p = 0.046 for overall survival). Functional analysis showed that these 12 key genes jointly influenced pyruvate metabolism, glycolysis/gluconeogenesis, and fatty acid degradation. Using paired expression and metabolite data, it was found that the metabolite lactate was significantly increased in cancer and was significantly negatively correlated with the expression of the core metabolic module genes, including ADH1A, ADH1B, and ADH1C. It is consistent with existing research showing increased lactate levels may be a marker for tumor aggressiveness, and patients with high lactate scores contribute to immune evasion and a decrease in the function of cytotoxic T cells in breast cancer [ref. 33,ref. 34]. Indeed, immune infiltration analysis showed that multiple immune cells were significantly different in the high and low metabolic activity groups. For example, the infiltration level of CD8 and CD4 T cells was significantly lower in the low metabolic activity group. Using the scRNA-seq cohort, the activity of the core metabolic module was significantly lower in malignant cells than in normal cells. Genes within this module, such as ALDH3A2, ALDH7A1, ALDH9A1, ADH5 and ALDH2 exhibited a significantly lower expression level in malignant cells and consistent results were also observed in cancer samples. The activity of the core metabolic module was significantly positively correlated with the degree of differentiation of CD8 T cells. Specifically, differential activity of this metabolic module was observed in CD8 T cell subtypes, with the stem-like CD8 T cell subtype showing the highest metabolic activity. The stem-like CD8 T cell subtype exhibited a higher expression of exhausted T cell marker genes and a lower the degree of cell differentiation. Research showed that metabolic limitations in the tumor microenvironment may exacerbate T cell exhaustion and lead to impaired anti-tumor immune responses [ref. 35]. In contrast to existing subpathway methods that rely solely on differential expression analysis or network topological structure, our CMSubpathway framework introduces a novel multi-dimensional integrative strategy that systematically identifies cancer-specific risk metabolic subpathways by integrating modeling metabolic network topology, expression dysregulation, prognostic association, and disease classification power. Beyond methodological improvements, we further revealed a novel immune–metabolic mechanism by which the core metabolic module is linked to lactate accumulation and is associated with CD8+ T cell differentiation and stem-like exhaustion. Therefore, we developed a novel method for identifying cancer risk metabolic subpathways. The identified metabolic reprogramming events were significantly associated with tumor immune cell infiltration and could act as potential biomarkers, which was conducive to the development of more effective cancer treatment strategies.

Finally, it is important to note that our findings regarding metabolic reprogramming and immune regulation are primarily based on correlative computational analyses. While these multi-omics associations are robust and consistently observed across single-cell and bulk datasets, they do not strictly establish direct causation. However, these theoretical findings can serve as a valuable theoretical foundation and be highly helpful in planning future laboratory analyses. Future in vitro and in vivo functional experiments are warranted to mechanistically validate the direct causal roles of this core module in driving metabolic reprogramming and shaping the immunosuppressive tumor microenvironment.

4. Materials and Methods

4.1. Data and Preprocessing

We downloaded the human metabolism Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway KGML files from the KEGG database. A total of 84 metabolism pathway KGML files were downloaded, and the R package “KEGGgraph”(v 1.72.0) [ref. 36] was used to import and convert them into the form of interactive networks. In total, 21 different The Cancer Genome Atlas Program (TCGA) projects, each representing a specific cancer type, were analyzed. RNA-seq-based gene expression profile data were downloaded from the TCGA project. The clinical information of tumor patients, including the survival status, stage, grade, and survival time, was also downloaded from the TCGA project. The scRNA-seq data of breast cancer were downloaded from GSE195861 via the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/, accessed on 14 March 2024) [ref. 26]. The bulk RNA-seq of breast cancer fat metabolism data were downloaded from GSE246231 via the Gene Expression Omnibus database. The bulk RNA-seq datasets of 9 breast cancer peripheral blood cohorts (including GSE111842, GSE109761, GSE111065, GSE51827, GSE55807, GSE67939, GSE75367, GSE86978, GSE41245) were collected from LncTarD 2.0 and Gene Expression Omnibus database [ref. 37]. The Seurat package (v 5.1.0) was used to load the 10X Genomics data from each sample into the R software (v 4.3.2). Using the R Harmony package (version 1.0), batch effects were eliminated based on the top 50 PCA components [ref. 38]. Based on harmony-corrected data, k-nearest neighbors (KNNs) were calculated, and a shared nearest neighbor (SNN) graph was constructed. The modular function was then modified based on the clustering algorithm to accomplish cluster recognition. The identified clusters were presented on the 2D map made with the t-distributed stochastic neighbor embedding (tSNE) or uniform manifold approximation and projection (UMAP) for the dimension reduction method [ref. 39]. For the TCGA cohort, gene expression levels were quantified as Transcripts Per Million (TPMs) and log2-transformed. For the GEO datasets, gene expression matrices were normalized using the Robust Multi-array Average (RMA) method and log2-transformed. This approach is suitable and widely applied to cross-sample comparisons, and differential expression analysis [ref. 40]. As a representative example, the TCGA-BRCA and GSE42568 cohorts were analyzed separately, and the consistently identified core metabolic module was used for subsequent analyses. Furthermore, multiple GEO datasets mentioned above were employed as independent validation sets for verification.

4.2. Identification of Cancer-Related Risk Metabolic Subpathways

To identify cancer-related risk metabolic subpathways (CMSubpathway), a step-by-step procedure was taken. First, the metabolic subpathways were located by means of the LTE (local tightness expansion) algorithm [ref. 41]. Second, dysregulated metabolic subpathways in cancer were identified. Third, candidate subpathways with prognostic efficacy were identified. Fourth, the classification performance of candidate subpathways was evaluated. In this study, we defined cancer-related risk metabolic subpathways that were dysregulated in cancer, correlated with prognosis, and facilitated the classification of normal and cancer samples. The schematic workflow is shown in Figure 1A. The detailed processes for the identification are as follows.

• Step 1: Identification of metabolic subpathways

First, we performed a fast greedy module community recognition algorithm LTE to extract the representative subpathways from KEGG metabolism pathways. For each network, this algorithm for community detection is operated by optimizing modularity, a measure of the strength of the division of each network into communities [ref. 41]. Each node started as its community in a network, and then iteratively merged communities if two communities were adjacent to increase the overall modularity. This process continues until the modularity has no further improvement. Modularity Q is defined as follows:

$A_{i,j}$ is the adjacency matrix of the network (1 if there is an edge between nodes $i$ and $j$, 0 otherwise). $n$ is the total number of edges in the network. $W_i$ and $W_j$ are the degrees of nodes $i$ and $j$.${ \delta }_{i,j}$ is a function that is 1 if nodes $i$ and $j$ belong to the same community, and 0 otherwise. In our study, the minimum number of nodes in the module is 5, and the minimum modularity is 0.3.

• Step 2: Identification of dysregulated metabolic subpathways in cancer

Second, to investigate the differences in metabolic subpathways between normal samples and cancer samples, we calculated the scores of metabolic pathways using GSEA. The R package “limma”(v 3.6.9) is used to calculate differential genes between normal and cancer samples. Differentially expressed genes are sorted in a descending rank list. For each metabolic subpathway, an enrichment score (ES) is calculated [ref. 42]. The ES reflects the degree to which the genes in the metabolic subpathway are overrepresented at the top (or bottom) of the ranked list. The formula for ES can be expressed as follows:

$r_i$ is the rank metric for the gene at position $i$ of the rank list, and $N$ is the total number of genes in which the subpathway is located. The ES is normalized to account for differences in subpathway sizes. This results in a normalized enrichment score (NES), which allows for comparison across different subpathways. To properly address multiple hypothesis testing and strictly control for false positive findings, the Benjamini–Hochberg method was applied. Furthermore, to mitigate potential biases arising from tumor heterogeneity, we implemented an iterative subsampling strategy [ref. 43]. Specifically, the process of identifying differentially expressed genes and calculating NESs was repeated 100 times. In each iteration, 70% of the total cancer and normal samples were randomly subsampled. A subpathway was ultimately defined as significantly dysregulated only if it met the criteria of abs (NES) > 1 and nominal p-value < 0.05 and false discovery rate (FDR) < 0.25 in at least 80 out of the 100 resampling iterations [ref. 44].

• Step 3: Evaluating the prognostic efficacy of metabolic subpathways in cancer

Third, to investigate the impact of subpathways on survival, we conducted a survival analysis on metabolic subpathways. We conducted gene set variation analysis (GSVA) on cancer samples. The GSVA score for each subpathway in each sample is computed using the following formula [ref. 45]:

$S_{ij}$ is the GSVA score for subpathway $j$ in sample $i$. $M_i$ is the expression level of gene $i$. $N$ is the number of genes in the subpathway. $K$ is a kernel function (Gaussian) that defines the weighting of expression differences.

GSVA scores were used to distinguish the high-score and low-score groups, using the R function “surv_cutpoint” to find the best threshold, and then compared the survival status of the two groups. To compare survival curves between the high-score and low-score groups, the log-rank test can be used to assess whether there are significant differences between them; p < 0.05 was considered significant.

• Step 4: Evaluating the classification performance of metabolic subpathways

Fourth, we constructed support vector machine (SVM) models to examine the classification performance of subpathways for normal and cancer samples. The SVM algorithm is widely used for cancer genomic classification and prediction [ref. 46]. We conducted five-fold cross-validation on a dataset of GSVA scores for each subpathway across various cancers in the pan-cancer TCGA dataset. We randomly divided the samples into five equal parts, selecting four for training (training set) and one for testing (test set). Each SVM model was built using the “e1071” (v 1.7.17) package in R [ref. 47] which provides an R interface to libsvm. The functions “svm()” and “predict()” were utilized to build each SVM model and to predict the sample type, respectively [ref. 48]. The outcome of the prediction was then assessed by area under the index curve (AUC) analysis. We considered subpathways with an AUC > 0.75 as candidates for classification, given that this value serves as an accepted benchmark for being considered a candidate subpathway for classification. The model has good discrimination performance [ref. 49].

4.3. Clustering and Cell Type Identification of scRNA-Seq

Using the “FindAllMarkers” function, we identified the marker genes for each cluster. Based on the DEGs and well-known cellular markers mentioned in CellMarker 2.0, the cell groups were annotated [ref. 50]. InferCNV (version 1.12.0) was utilized to estimate single-cell copy number variation (CNV) profiles to identify malignant cells [ref. 51]. Immune cells were used as reference cells to define a baseline of CNV estimates. Finally, compared with the reference cells, subclusters with significantly higher CNV scores were defined as malignant cells (Wilcoxon tests, p < 0.05, fold changes > 1.5).

4.4. CRISPR Knockdown

We retrieved the CRISPRGeneEffect dataset from the DepMap database (https://depmap.org/portal/, accessed on 8 April 2026). The CRISPR Essentiality Scores (CERESs) were utilized to calculate gene dependency, aiming to pinpoint genes crucial for cellular proliferation and survival. Specifically, a CERES score > 0 indicates that gene disruption enhances cell line viability, whereas a score < 0 implies that the knockout inhibits cell survival and proliferation [ref. 52].

4.5. Calculating the Activity of Key Metabolic Subpathways in Single-Cell RNA-Seq

The Single Pathway analysis in Single Cells (SiPSiC) (v 1.4.2) package was used to infer cancer key metabolic subpathway activity of the single-cell RNA-seq dataset [ref. 53]. We provided a list of genes contained by the cancer key metabolic subpathway, and then per-cell cancer key metabolic subpathway activity was calculated for all cells by the “getPathwayScores” function.

4.6. Predicting Differentiation States of T Cells from scRNA-Seq Data

CytoTRACE (version 0.3.3) is a computational framework for recreating the relative differentiation state of single-cell RNA sequencing data using the gene expression profile [ref. 54]. The differentiation state of cells may be deduced from scRNA-seq data without any previous knowledge. Initially, a KNN graph containing information on the undirected relationships between cells was built. CytoTRACE was then utilized to determine the recommended time order of cells. The KNN graph and planned time were then utilized to generate a transfer matrix, which was afterwards shown on a UMAP scatter plot.

4.7. Protein Level Verification

The CPTAC protein dataset was utilized to validate the results at the proteomic level, including 18 breast cancer samples and 118 normal samples. Protein expression patterns in both normal and tumor tissues were ascertained by retrieving images from the Human Protein Atlas (HPA, https://www.proteinatlas.org) database. The immunohistochemical (IHC) data of key genes in breast cancer and normal breast tissues were obtained through the HPA database.

References

1. B. Faubert, A. Solmonson, R.J. DeBerardinis. Metabolic reprogramming and cancer progression. Science, 2020. [DOI | PubMed]
2. Y. Zhang, M.-J. Wu, W.-C. Lu, Y.-C. Li, C.J. Chang, J.-Y. Yang. Metabolic switch regulates lineage plasticity and induces synthetic lethality in triple-negative breast cancer. Cell Metab., 2024. [DOI | PubMed]
3. L. Torresano, C. Nuevo-Tapioles, F. Santacatterina, J.M. Cuezva. Metabolic reprogramming and disease progression in cancer patients. Biochim. Biophys. Acta BBA Mol. Basis Dis., 2020. [DOI]
4. Z. Lin, J. Li, J. Zhang, W. Feng, J. Lu, X. Ma, W. Ding, S. Ouyang, J. Lu, P. Yue. Metabolic reprogramming driven by IGF2BP3 promotes acquired resistance to EGFR inhibitors in non–small cell lung cancer. Cancer Res., 2023. [DOI | PubMed]
5. X. Chen, Z. Wan, L. Yang, S. Song, Z. Fu, K. Tang, L. Chen, Y. Song. Exosomes derived from reparative M2-like macrophages prevent bone loss in murine periodontitis models via IL-10 mRNA. J. Nanobiotechnol., 2022. [DOI | PubMed]
6. Y. Chen, M. Hu, L. Wang, W. Chen. Macrophage M1/M2 polarization. Eur. J. Pharmacol., 2020. [DOI | PubMed]
7. J. Zhang, J. Muri, G. Fitzgerald, T. Gorski, R. Gianni-Barrera, E. Masschelein, G. D’hUlst, P. Gilardoni, G. Turiel, Z. Fan. Endothelial lactate controls muscle regeneration from ischemia by inducing M2-like macrophage polarization. Cell Metab., 2020. [DOI | PubMed]
8. C. Li, J. Han, Q. Yao, C. Zou, Y. Xu, C. Zhang, D. Shang, L. Zhou, C. Zou, Z. Sun. Subpathway-GM: Identification of metabolic subpathways via joint power of interesting genes and metabolites and their topologies within pathways. Nucleic Acids Res., 2013. [DOI | PubMed]
9. Y. Sheng, Y. Jiang, Y. Yang, X. Li, J. Qiu, J. Wu, L. Cheng, J. Han. CNA2Subpathway: Identification of dysregulated subpathway driven by copy number alterations in cancer. Brief. Bioinform., 2021. [DOI | PubMed]
10. X. Han, D. Wang, P. Zhao, C. Liu, Y. Hao, L. Chang, J. Zhao, W. Zhao, L. Mu, J. Wang. Inference of subpathway activity profiles reveals metabolism abnormal subpathway regions in glioblastoma multiforme. Front. Oncol., 2020. [DOI | PubMed]
11. I. Elia, M. Rossi, S. Stegen, D. Broekaert, G. Doglioni, M. van Gorsel, R. Boon, C. Escalona-Noguero, S. Torrekens, C. Verfaillie. Breast cancer cells rely on environmental pyruvate to shape the metastatic niche. Nature, 2019. [DOI | PubMed]
12. J. Hu, S.-G. Wang, Y. Hou, Z. Chen, L. Liu, R. Li, N. Li, L. Zhou, Y. Yang, L. Wang. Multi-omic profiling of clear cell renal cell carcinoma identifies metabolic reprogramming associated with disease progression. Nat. Genet., 2024. [DOI | PubMed]
13. C. Feng, J. Zhang, X. Li, B. Ai, J. Han, Q. Wang, T. Wei, Y. Xu, M. Li, S. Li. Subpathway-CorSP: Identification of metabolic subpathways via integrating expression correlations and topological features between metabolites and genes of interest within pathways. Sci. Rep., 2016. [DOI | PubMed]
14. E. Benedetti, E.M. Liu, C. Tang, F. Kuo, M. Buyukozkan, T. Park, J. Park, F. Correa, A.A. Hakimi, A.M. Intlekofer. A multimodal atlas of tumour metabolism reveals the architecture of gene–metabolite covariation. Nat. Metab., 2023. [DOI | PubMed]
15. S. Hui, J.M. Ghergurovich, R.J. Morscher, C. Jang, X. Teng, W. Lu, L.A. Esparza, T. Reya, L. Zhan, Y.J. Guo. Glucose feeds the TCA cycle via circulating lactate. Nature, 2017. [DOI | PubMed]
16. B. Faubert, K.Y. Li, L. Cai, C.T. Hensley, J. Kim, L.G. Zacharias, C. Yang, Q.N. Do, S. Doucette, D. Burguete. Lactate metabolism in human lung tumors. Cell, 2017. [DOI | PubMed]
17. M.V. Liberti, J.W. Locasale. The Warburg effect: How does it benefit cancer cells?. Trends Biochem. Sci., 2016. [DOI | PubMed]
18. H. Lu, R.A. Forbes, A. Verma. Hypoxia-inducible factor 1 activation by aerobic glycolysis implicates the Warburg effect in carcinogenesis. J. Biol. Chem., 2002. [DOI | PubMed]
19. P. Sonveaux, F. Végran, T. Schroeder, M.C. Wergin, J. Verrax, Z.N. Rabbani, C.J. De Saedeleer, K.M. Kennedy, C. Diepart, B.F. Jordan. Targeting lactate-fueled respiration selectively kills hypoxic tumor cells in mice. J. Clin. Investig., 2008. [DOI | PubMed]
20. T. Fukami, T. Yokoi, M. Nakajima. Non-P450 drug-metabolizing enzymes: Contribution to drug disposition, toxicity, and development. Annu. Rev. Pharmacol. Toxicol., 2022. [DOI | PubMed]
21. C. Jiang, R. Liu, X. Wu. Alcohol dehydrogenase-1B represses the proliferation, invasion and migration of breast cancer cells by inactivating the mitogen-activated protein kinase signalling pathway. J. Physiol. Pharmacol., 2023
22. Y. Xue, M. Wang, D. Zhong, N. Tong, H. Chu, X. Sheng, Z. Zhang. ADH1C Ile350Val polymorphism and cancer risk: Evidence from 35 case–control studies. PLoS ONE, 2012. [DOI | PubMed]
23. X. Zheng, H. Ma, J. Wang, M. Huang, D. Fu, L. Qin, Q. Yin. Energy metabolism pathways in breast cancer progression: The reprogramming, crosstalk, and potential therapeutic targets. Transl. Oncol., 2022. [DOI | PubMed]
24. K. Goetze, S. Walenta, M. Ksiazkiewicz, L.A. Kunz-Schughart, W. Mueller-Klieser. Lactate enhances motility of tumor cells and inhibits monocyte migration and cytokine release. Int. J. Oncol., 2011. [DOI | PubMed]
25. M. Sukumar, J. Liu, Y. Ji, M. Subramanian, J.G. Crompton, Z. Yu, R. Roychoudhuri, D.C. Palmer, P. Muranski, E.D. Karoly. Inhibiting glycolytic metabolism enhances CD8+ T cell memory and antitumor function. J. Clin. Investig., 2013. [DOI | PubMed]
26. M. Tokura, J. Nakayama, M. Prieto-Vila, S. Shiino, M. Yoshida, T. Yamamoto, N. Watanabe, S. Takayama, Y. Suzuki, K. Okamoto. Single-cell transcriptome profiling reveals intratumoral heterogeneity and molecular features of ductal carcinoma in situ. Cancer Res., 2022. [DOI | PubMed]
27. L. Wang, W. Liu, K. Liu, L. Wang, X. Yin, L. Bo, H. Xu, S. Lin, K. Feng, X. Zhou. The dynamic dysregulated network identifies stage-specific markers during lung adenocarcinoma malignant progression and metastasis. Mol. Ther. Nucleic Acids, 2022. [DOI | PubMed]
28. C.S. Jansen, N. Prokhnevska, V.A. Master, M.G. Sanda, J.W. Carlisle, M.A. Bilen, M. Cardenas, S. Wilkinson, R. Lake, A.G. Sowalsky. An intra-tumoral niche maintains and differentiates stem-like CD8 T cells. Nature, 2019. [DOI | PubMed]
29. F. Franco, A. Jaccard, P. Romero, Y.-R. Yu, P.-C. Ho. Metabolic and epigenetic regulation of T-cell exhaustion. Nat. Metab., 2020. [DOI | PubMed]
30. K.A. Brown. Metabolic pathways in obesity-related breast cancer. Nat. Rev. Endocrinol., 2021. [DOI | PubMed]
31. J.D. Chan, J. Lai, C.Y. Slaney, A. Kallies, P.A. Beavis, P.K. Darcy. Cellular networks controlling T cell persistence in adoptive cell therapy. Nat. Rev. Immunol., 2021. [DOI | PubMed]
32. A. Kallies, D. Zehn, D.T. Utzschneider. Precursor exhausted T cells: Key to successful immunotherapy?. Nat. Rev. Immunol., 2020. [DOI | PubMed]
33. V. Frisardi, S. Canovi, S. Vaccaro, R. Frazzi. The significance of microenvironmental and circulating lactate in breast cancer. Int. J. Mol. Sci., 2023. [DOI | PubMed]
34. G.A. Brooks. Lactate as a fulcrum of metabolism. Redox Biol., 2020. [DOI | PubMed]
35. L. Xia, L. Oyang, J. Lin, S. Tan, Y. Han, N. Wu, P. Yi, L. Tang, Q. Pan, S. Rao. The cancer metabolic reprogramming and immune response. Mol. Cancer, 2021. [DOI | PubMed]
36. J.D. Zhang, S. Wiemann. KEGGgraph: A graph approach to KEGG PATHWAY in R and bioconductor. Bioinformatics, 2009. [DOI | PubMed]
37. H. Zhao, X. Yin, H. Xu, K. Liu, W. Liu, L. Wang, C. Zhang, L. Bo, X. Lan, S. Lin. LncTarD 2.0: An updated comprehensive database for experimentally-supported functional lncRNA–target regulations in human diseases. Nucleic Acids Res., 2023. [DOI | PubMed]
38. Y. Zhou, D. Yang, Q. Yang, X. Lv, W. Huang, Z. Zhou, Y. Wang, Z. Zhang, T. Yuan, X. Ding. Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma. Nat. Commun., 2020. [DOI | PubMed]
39. Z. Zhao, D. Zhang, F. Yang, M. Xu, S. Zhao, T. Pan, C. Liu, Y. Liu, Q. Wu, Q. Tu. Evolutionarily conservative and non-conservative regulatory networks during primate interneuron development revealed by single-cell RNA and ATAC sequencing. Cell Res., 2022. [DOI | PubMed]
40. B. Li, V. Ruotti, R.M. Stewart, J.A. Thomson, C.N. Dewey. RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 2010. [DOI | PubMed]
41. J. Huang, H. Sun, Y. Liu, Q. Song, T. Weninger. Towards online multiresolution community detection in large-scale networks. PLoS ONE, 2011. [DOI | PubMed]
42. J. Candia, L. Ferrucci. Assessment of Gene Set Enrichment Analysis using curated RNA-seq-based benchmarks. PLoS ONE, 2024. [DOI | PubMed]
43. S. Acharya, P. Kossinna, Q. Zhang, J. Guo. Stabilized marker gene identification and functional annotation from single-cell transcriptomic data. PLoS Comput. Biol., 2025. [DOI | PubMed]
44. Z.X. Xiao, R. Liang, Y. Liu, C. Huang, Q. Fang, X. Hu, J. Wang, N. Olsen, D. Wu, S.G. Zheng. Regulation and Mechanism of Deleted in Breast Cancer-1 on Dendritic Cell Function in Systemic Lupus Erythematosus. MedComm, 2026. [DOI | PubMed]
45. S. Hänzelmann, R. Castelo, J. Guinney. GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinform., 2013. [DOI | PubMed]
46. S. Huang, N. Cai, P.P. Pacheco, S. Narrandes, Y. Wang, W. Xu. Applications of support vector machine (SVM) learning in cancer genomics. Cancer Genom. Proteom., 2018
47. D. Meyer, E. Dimitriadou, K. Hornik, A. Weingessel, F. Leisch. e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien. 2025
48. Z. Ning, S. Yu, Y. Zhao, X. Sun, H. Wu, X. Yu. Identification of miRNA-mediated subpathways as prostate cancer biomarkers based on topological inference in a machine learning process using integrated gene and miRNA expression data. Front. Genet., 2021. [DOI | PubMed]
49. Y. Xia, C. Fan, K.A. Hoadley, J.S. Parker, C.M. Perou. Genetic determinants of the molecular portraits of epithelial cancers. Nat. Commun., 2019. [DOI | PubMed]
50. C. Hu, T. Li, Y. Xu, X. Zhang, F. Li, J. Bai, J. Chen, W. Jiang, K. Yang, Q. Ou. CellMarker 2.0: An updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res., 2023. [DOI | PubMed]
51. A.P. Patel, I. Tirosh, J.J. Trombetta, A.K. Shalek, S.M. Gillespie, H. Wakimoto, D.P. Cahill, B.V. Nahed, W.T. Curry, R.L. Martuza. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science, 2014. [DOI | PubMed]
52. B. Li, Y. Pan, J. Wu, C. Miao, Z. Wang. Large-scale genomic-wide CRISPR screening revealed PRC1 as a tumor essential candidate in clear cell renal cell carcinoma. Int. J. Med. Sci., 2025. [DOI | PubMed]
53. D. Davis, A. Wizel, Y. Drier. Accurate estimation of pathway activity in single cells for clustering and differential analysis. Genome Res., 2024. [DOI | PubMed]
54. G.S. Gulati, S.S. Sikandar, D.J. Wesche, A. Manjunath, A. Bharadwaj, M.J. Berger, F. Ilagan, A.H. Kuo, R.W. Hsieh, S. Cai. Single-cell transcriptional diversity is a hallmark of developmental potential. Science, 2020. [DOI | PubMed]