BRCA gene expression profiles were downloaded from five independent cohorts of patients, including TCGA-BRCA, GSE25066, GSE21653, GSE103091and METABRIC. Only samples from TNBC were reserved in all cohorts. Survival analysis only considered overall survival (OS) and disease-free survival (DFS). In the above five cohorts, METABRIC had no patient prognostic information. The remaining histological data were obtained from the TCGA-BRCA cohort, including copy number variant data obtained via firehose, and mutation MAF files obtained from the cBioPortal Pancancer Project. TCGA-BRCA partial samples of the predicted neoantigen numbers were obtained from published literature (Rooney et al., 2015). The metabolic gene file used for clustering (Possemato et al., 2011), the metabolism signatures (Rosario et al., 2018), the immune pathway signatures (Bindea et al., 2013) and the oncogenetic signature.txt (Sanchez-Vega et al., 2018) from different published literature. The drug information is from the GDSC database involving the drug’s R package “pRRophetic” for use in predicting the drug’s IC50. The external datasets were used to determine whether the defined subtypes are likely to respond to immunotherapy (Roh et al., 2017).

Limma package used for identifying DEGs (|log2FC| > 1 and p < 0.01). Genetic feature set files “c2.cp.kegg.v6.2.symbols.gmt” and “h,all,v60.2.symbols” were obtained from the Molecular Signature Database (MSigDB). Then, Clusterprofiler R package was used for Gene Set Enrichment Analysis (GSEA) (Yu et al., 2012).

To identify Aldo-Keto reductase family 1 member B10 (AKR1B10), we collected peripheral blood samples from 30 TNBC patients and 30 healthy individuals as controls from Tianjin Medical University Cancer Institute and Hospital (Tianjin, China) in 2022 for RT-qPCR, and their paraffin-embedded tissues for IHC. All patients were female patients who were recently admitted and had not undergone radiotherapy, chemotherapy, or surgery. Control group was determined to be free from TNBC and other malignant tumors.

Because all data used in this study were derived from five platforms, and some of the data were normalized, we combined the data after normalizing each data using z-score to eliminate potential batch effects. We performed consensus NMF (Possemato et al., 2011) with 2-5 cluster numbers using TCGA data expression profiles and calculated the covariance coefficients for each decomposition. The MOVICS package (Lu et al., 2021) was used for differential expression analysis of these two subtypes, while the top 50 most significantly upregulated genes in each subtype were used as biomarkers for the different subtypes (p < 0.05, FDR < 0.25). In addition, we constructed a template using MOVICS.

Genomic variation analysis (GSVA) is a method genome augmentation, which calculates the characteristics of certain pathways or different populations based on expression spectra. The differences of gene sets between samples were investigated by GSVA R software package from relevant metabolic pathway gene sets (Rosario et al., 2018). Then, the limma package (Liu et al., 2019) was used to obtain the substitution gene scores, analyze the differences, and screen for DEGs features.

To identify the extent to which genes are up or downregulated within a single sample, ssGSEA is used for quantifying the immune composition of tumors. Here, we assessed the enrichment fraction of gene sets representing biological processes as well as biological pathways in bulk tumors or individual cancer cells by ssGSEA.

The ESTIMATE algorithm (Yoshihara et al., 2013) can be applied to calculate the permeability and matrix content of immune cells (immune fraction) and stromal content (stromal fraction) of different subtypes, thus reflecting the microenvironmental characteristics of tumors. Microenvironmental Cell Population counter (MCPcounter) (Becht et al., 2016) was used for evaluating the penetration frequency of immune and non-immune cell populations in two subtypes.

The detection of co-occurrence and mutually exclusive mutations mainly relied on the CoMEt algorithm. Next, we predicted the different genotypes between different subtypes, including the number of neoantigens, TMB, copy number amplification, and the frequency of copy number deletions. We also performed an online analysis using GISTIC2 (Cibulskis et al., 2013) to obtain the number of amplifications and deletions for all samples and to calculate arm- and focal level somatic copy-number alterations (SCNAs) and G-scores in tumors with the input of “SNP6” files.

MD-Anderson melanoma cohort treated with anti-CTLA-4 or anti-PD-1 is considered to be used to predict immunotherapy response (Roh et al., 2017). And then, we analyzed the sub map from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Roh et al., 2017) and studied the sensitivity differences between the C1 and C2 groups after multiple drug treatments.

IHC staining was used to slice the dewaxed tissue portion of TNBC samples and cure with 3% hydrogen peroxide for a period of time. Block endogenous peroxidase for 30 min, then solidify with appropriate horseradish AKR1B10 antibody. The IHC fraction is calculated by multiplying the dyeing intensity by the percentage of cells. Definition of intensity: 0 (unstained), 1 (soft), 2 (medium), 3 (strong). Definition of percentage of cells: 1 (25%), 2 (26%–50%), 3 (51%–75%), and 4 (>75%). More than 3 was defined as positive, while less than or equal to 3 as negative.

Triazole solution (AC0101-B; SparkJade, China) was used for extracting RNA from blood and tissues samples. 2×HQ SYBR qPCR Mix (ZF501; ZOMANBIO; Beijing, China) was used for PCR reaction. Primer sequences were listed in Supplementary Table S1. The levels of AKR1B10 expression were calculated by the method of 2−∆∆Cq.

R software (version 4.0.2) was used to process all data. Contingency table (χ2) variables used the chi‐square test and Fisher’s precision probability test for statistical significance. Kaplan-Meier method was used for survival analysis and compared the results by the log-rank test. Z test was used to assess whether there was a significant difference between the two groups. Univariate Cox proportional hazards regression models were used to assess the risk ratio for univariate analysis. A two-tailed p-value < 0.05 was considered statistically significant.

Based on TCGA database and NMF algorithm analysis, we divided TNBC into two subtypes with different metabolic characteristics. In this study, 115 cases of TCGA-TNBC patients were screened. Supplementary Table S2 showed clinical characteristics of TNBC patients. Before analyzing the TNBC NMF algorithm, we used the ComBat algorithm to eliminate batch processing effects in the TNBC queue. And after deleting the batch processing effect, draw a key element analysis diagram (Figure 1A). Previously, a total of 2,752 reported metabolic related genes (Possemato et al., 2011) were screened and downloaded as the basis for analyzing metabolomics in our study (Supplementary Table S3).

To identify subtypes in TNBC, Cox regression was used. A total of 637 prognostic genes were obtained (Supplementary Table S4). After a further adjusted p-value (p < 0.05), 277 candidate genes were identified (Supplementary Table S5). We then used the NMF algorithm to cluster the 277 candidate genes and drew the NMF with two to five sets (Figure 1B). After comprehensive consideration, the optimal clustering number (k value) was 2, defining two subtypes C1 (n = 80) and C2 (n = 35). To verify the consistency between subtype designations and two-dimensional distribution patterns, we reduced the PCA dimension (Figure 1C). Subsequently, the same conclusion was validated in the validation set (GSE25066, GSE21653, and GSE103091, excluding METABRIC with a lack of patient OS).

Finally, two TNBC molecular subtypes were established. We also used the survival information in the four queues to analyze the subtype survival of TNBC subsets. The OS of C1 was verified better than that of C2 in TCGA-TNBC patients (p = 0.014, Figure 1D) and other datasets patients (GSE25066 and GSE103091) (p = 0.014 and p = 0.017, respectively) (Figures 1E, G). No significant difference was observed in datasets (GSE21653) (p = 0.15) (Figure 1F).

In this study, we analyzed whether different TNBC subtypes have their own characteristics in distinct metabolic pathways. Firstly, we used the GSVA R package to score metabolic pathways (Rosario et al., 2018) (Supplementary Table S6). Limma difference test cross group was performed to confirm subtype-specific differential metabolic pathways, and heatmaps were constructed for visualization (Supplementary Figure S1).

Furthermore, the DEGs between two groups were detected by GSVA enrichment again, and it was found that C1 and C2 had specific metabolic characteristics (Figure 2A). There were 17 metabolism-related pathways in C2 that were significantly upregulated, mainly involving pentose and glucuronate interconversion, steroid hormone biosynthesis, tyrosine metabolism, oxidative phosphorylation, ketone biosynthesis and metabolism. Similar outcomes that 17 metabolism-related pathways were activated in C2 were observed in the validation datasets (GSE25066, GSE21653, GSE103091 and METABRIC) (Supplementary Figures S2, S3). To determine the different activities of metabolic pathways, we represented the two subtypes in the TCGA-TNBC cohorts and validation cohorts and revealed that C2 contained the highest activation of metabolic pathways in all cohorts (Supplementary Figure S4).

In analyze the differences between the two subtypes in carcinogenesis-related pathways, we counted the GSVA enrichment points and plotted box-line plots. Eleven carcinogenesis-related pathways were selected and quantified. The results showed that different subgroups were closely related to the activation of different carcinogenic signaling pathways, which mainly involved cell cycle, PI3K, RTK-RAS, and angiogenesis (Figure 2B). C2 displayed significant activation of the PI3K and angiogenesis pathways. C1 had a stronger cell cycle, HIPPO, RTK-RAS and WNT signature than C2. These differences in carcinogenic pathway activity may affect their prognosis. After evaluating whether the subtype was related to the tumor microenvironment, it was found that the immune score of C1 was higher than that of C2 (p = 0.046), and the stromal score had no significant difference (Figure 2C).

To evaluate the immune status of two subtypes, the MCP counter and ssGSEA algorithm were used to estimate the abundance of immune cells (Figure 3A). The results showed significant differences between different immune cell groups between the two subtypes (Figures 3B, C). Especially, the immune value of C2 in most immune cells was obvious lower than that of C1, except for neutrophils, fibroblasts and Th17 cells. According to this study, C1 was rich in more immune cells and had the highest immune score, which indicated that differences in the distribution of different immune cells may be the reason for the poorer prognosis of C2 than C1.

To explore the relationship between these subtypes and clinical features, we analyzed the clinicopathological parameters between the two subtypes and constructed a clinical information heatmap of subtypes (Figure 3A). The results revealed that larger tumor size (p = 0.007) and advanced pathologic stage (TNM III/IV stage) (p = 0.001) were related to the C2 subtype (Supplementary Table S7). We also constructed a clinicopathological variables heatmap of subtypes in the validation cohorts and presented detailed data (Supplementary Figure S2; Supplementary Table S2). It is well known that larger tumor size and advanced TNM stage represent shorter survival in TNBC (Johansson et al., 2021).

Breast cancer has been closely related to many genomic mutations in the body (Kim et al., 2021). To investigate the difference of somatic mutations frequency between TNBC subtypes, we applied specific driver mutations for breast cancer (Bailey et al., 2018) to estimate gene mutations and draw waterfall map. High mutation frequencies of TP53, BRCA1, PIK3CA, PTEN, FBXW7, NF1, RB1, KMT2C, and PTPRD in both TNBC subtypes were observed (Figure 4A; Supplementary Table S8). We found that C2 exhibited different mutation characteristics from C1. Specifically, C2 has a higher mutation frequency, such as PIK3CA, KMT2D, KMT2C, and so on (Figure 5A; Supplementary Table S9). We calculated the TMB for each metabolic subtype (Figure 4B). Although there was no difference (p = 0.16), a trend showed that the TMB of C1 was higher than that of C2. We also analyzed the total number of mutations and expected neoantigens (Figure 4C) and observed a significant difference between them (p = 0.0033). Subsequently, the frequency of amplification (Figure 4D) and deletion (Figure 4E) was showed and found that patients within C1 only showed higher amplification than C2 (p = 0.043).

Finally, we mapped a cell column to change the number of copies of each group by performing online GISTIC2.0 analysis, in which red represented gains and blue represented losses (Figures 5B, C). Both C1 and C2 observed copy number alterations in chromosome regions, including amplification at 11p13 and deletion at 8p23.2, 9p21.3, 13q14.2, and 19p13.3. In contrast to C1, C2 has significant amplification at 1q23.3 and 7p11.2. These differences could also explain that C2 has a better prognosis than C1. Therefore, changes in copy quantity might be the main mechanism behind the differences in metabolism and prognosis between the two group.

The difference in sensitivity to chemotherapeutic drugs and targeted drugs between two groups was analyzed by using the GDSC drug sensitivity database. The top 12 drugs with differential responses were plotted and listed (Figure 6A). After estimating the IC50 value, we found that C2 may be less sensitive to chemotherapy, including bleomycin, vinorelbine and doxorubicin (all p < 0.001).

The different immune infiltration patterns among TNBC subtypes suggested that further research on the response of immunotherapy was needed. To this aim, we matched the expression spectra of two subspecies to determine the similarity of the TCGA reaction spectra (Figure 6B). The results indicate that C1 may be more sensitive to PD1 treatment (p = 0.02), and C2 may have a better therapeutic effect on anti‐CTLA4.

And then, we extracted the top 50 genes of each metabolic specificity as biomarkers and constructed clinical models, and plotted correlation heatmaps using MOVICS (Bailey et al., 2018) package analysis. The classifier based on 100 genes was generated and visualized by heatmap (Figure 7A; Supplementary Table S10). In order to predict the identification of metabolic subtypes in each sample, we conducted consistency testing on the results of the two subtypes using the NTP algorithm and indicated that the characteristics of these genes can be replicated to determine the TNBC type (Figure 7B).

To better distinguish the two subtypes, we assume that AKR1B10 was an effective biomarker for C2. RT-qPCR and IHC staining were used to preliminarily verify this hypothesis. AKR1B10 were overexpressed in peripheral blood of TNBC patients than in healthy control (Figure 7C). IHC showed that AKR1B10 were positive in 11 cases of TNBC, with a positive rate of 36.7% (Figures 7D, G). The average tumor size of AKR1B10 positive group was 2.2 cm from 0.8 to 3.5 cm, which was higher than that of negative group (Figure 7E). In addition, the number of lymph node metastasis in AKR1B10 positive group were more than that in negative group (Figure 7F). Large tumor and many lymph node metastases often indicate poor prognosis of TNBC, which was consistent with C2 subtype. These results were listed in Supplementary Table S11.