Two BC tissue specimens were obtained from patients at the Second Hospital of Dalian Medical University. Invasive breast cancer was pathologically confirmed in all patients not on chemotherapy or radiotherapy before tissue collection. Written informed consent was obtained from patients, and the study was approved by the Ethics and Human Subject Committee of the Second Hospital of Dalian Medical University (NO.2023191). Procedures were performed according to hospital guidelines and regulations.

Gene expression matrices of enrolled patients were obtained from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases. We included 1,069 BC samples from TCGA as the training cohort. The gene-expression profiles of TCGA-BRCA in the Fragments Per Kilobase per Million (FPKM) format were obtained from the TCGA portal (http://cancergenome.nih.gov), and then the ID conversion was carried out through the operation of ENSG ID to GeneSymbol, and finally the data standardization was carried out, and the standardization method is log2 (X+1). In addition, the BC patients’ clinical data (gender, age, histological type, and survival) were downloaded from TCGA. After searching the datasets with more than 150 human breast cancer samples with complete expression profile data, we selected the GSE42568, GSE88770, GSE48390, and GSE162228 dataset from the GEO as the validation cohort. These datasets were verified using the GPL570 platform. To ensure the scientificity and accuracy of the research, we successfully removed batch effect with COMBAT when combining GEO multi-data sets ( Supplementary Figure S1 ). Additionally, clinical survival and outcome data of BC patients were also downloaded from this database.

Data analysis was performed using the “limma” R package. Fold change > 1.5, p < 0.05, and false discovery rate (FDR) < 0.05 were set as the cutoffs to screen for DEGs.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to enrich the DEGs into associated pathways using the “clusterProfiler” R package (version 3.14.3). p < 0.05 and FDR < 0.05 were considered significant.

Least absolute shrinkage and selection operator (LASSO) Cox regression analysis identified genes most correlated with OS, and 10-round cross-validation was performed to prevent overfitting. The risk score for each patient was then calculated based on the expression levels of genes. Risk score: -0.0419970982477039 * NPY1R - 0.162055812415471 * CELSR2 - 0.043004672153174 * STC2 - 0.0716026845406244 * SCUBE2 + 0.2810654696502 * GIMAP2 + 0.0773881988402307 * HLA-DPB1 - 0.0232515777318596 * CXCL14 - 0.721867840891611 * KLRB1 - 0.253187064109637 * BIRC3 - 0.0587584464454724 * IL18 - 0.242105852075788 * PSMB8 + 0.198881881356143 * CD1C + 0.0814403392760682 * TNFAIP8 + 0.076656198308623 * IRF1. According to the median risk score, BC patients were divided into high- and low-risk groups. Kaplan–Meier analysis was employed to estimate the difference in OS between the categorized patients via the R package “survival.” The prognostic capability of the risk model was validated using time- dependent receiver operating characteristic (ROC) analysis with the R package “pROC”.

Univariate and multivariable Cox regression analyses were performed to estimate whether the risk score was an independent predictor of BC prognosis. A subgroup analysis was conducted to confirm the independence of the risk model. The patients with BC in the training cohort were regrouped into new subgroups based on different clinical characteristics, and the patients in each subgroup were stratified into high- and low-risk groups, based on the median risk score.

Patient tissue specimens were fixed in 10% neutral formalin, embedded in paraffin, and sectioned into 4 µm sections before staining. Sections were deparaffinized, rehydrated, and blocked for endogenous peroxidase activity. Next, antigen retrieval was performed in citrate buffer (pH 6.0) and sections autoclaved for 90 s at 121°C. After washing in phosphate buffered saline (3 min × 3), sections were blocked in goat serum at room temperature for 30 min and incubated with primary antibodies (PSMB8, (1:200), Proteintech Group, IL, USA; cIAP2, (1:200), Proteintech Group, IL, USA) overnight at 4° C. The next day, sections were incubated with secondary antibodies (Maxin Biotechnologies, China) and treated with diaminobenzidine hydrochloride to visualize immunoreactivity. The immunohistochemical scoring was performed independently by two experienced pathologists, who had no knowledge of the clinicopathological information.

Nomograms are user-friendly clinical tools used to predict disease prognosis. The risk score and clinical parameters were subjected to univariate Cox regression analysis, and features with P values < 0.05 were subjected to multivariable COX regression analysis. Features with p values < 0.05 after multivariate analysis were incorporated into nomograms that were constructed to predict the 3- and 5-year OS rates. The nomogram was based on three independent prognostic factors: age, tumor stage, and the risk score. Factors corresponded to a specific point by drawing a line straight up to the point axis. The sum of the three factor points indicated the total points. By drawing a perpendicular line from the total point axis to the two-outcome axes, estimated 3- and 5-year OS probabilities were obtained. Observed 3- and 5-year OS rates were compared with predicted rates to further verify predictive performance. We assessed nomogram goodness-of-fit using calibration plots.

The estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) method was applied to calculate the immune score, stromal score, and ESTIMATE score of the patients, via the R package “estimate”. Tumor immune estimation resource (TIMER) analysis was conducted to evaluate the abundance of six types of immune cells (neutrophils, CD4 T cells, macrophages, CD8 T cells, dendritic cells (DCs), and B cells). The MCPcounter (microenvironment cell populations-counter) algorithm was also used to assess T cell, CD8 T cell, cytotoxic lymphocyte, B cell lineage, natural killer (NK) cell, monocytic cell lineage, myeloid DC, neutrophil, endothelial cell, and fibroblast abundance.

TIMER was used to analyze correlations between hub gene expression and the degree of lymphocyte infiltration. TISIDB was also used to analyze correlations between hub gene expression and immune molecule expression in BC. We used the GSCA Lite (A Web Server for Gene Set Cancer Analysis: http://bioinfo.life.hust.edu.cn/web/GSCALite) online tool to analyze the correlation between hub genes expression and sensitivity to current chemotherapeutic or targeted drugs for BC.

Statistical analyses were completed using R (version 3.6.3). Discontinuous data were presented as number (percentage), and continuous data were displayed as mean± standard deviation. The Wilcoxon rank sum test was utilized to compare two groups and the Kruskal-Wallis test to compare multiple groups. In addition, the survfit function of “survival” package in R was used to analyze the prognostic differences between the two groups, and the log-rank test was used to further evaluate the significance of prognostic differences between the two groups. Statistical significance was defined as p < 0.05.

We downloaded the gene expression profiles and clinical information of 1,069 BC patients from The Cancer Genome Atlas (TCGA). Based on gene expression, BC can be mainly classified into Luminal A, Luminal B, HER2-enriched, Basal-like, and Normal-like (35, 36). The ESTIMATE algorithm showed that the highest mean immune score of Normal-like subtype was highest among all five subtypes, followed by Basal-like subtype, HER2-enriched subtype, and Luminal A subtype. The Luminal B subtype cases had the lowest immune scores ( Supplementary Figure S2A , p < 0.0001). However, stromal scores, from high to low, were Normal-like, Luminal A, HER2-enriched, Luminal B, and Basal-like ( Supplementary Figure S2B , p < 0.0001). The mammary gland is a hormone-responsive organ- the endocrine system is closely related to its development and disease occurrence, therefore we performed correlation analyses between immune and stromal scores and hormone receptor status. As shown in Supplementary Figure S2C , patients with progesterone receptor positive (PR+) had lower immune scores when compared with progesterone receptor negative (PR-) patients (p < 0.01), while estrogen receptor positive (ER+) patients had lower scores when compared with estrogen receptor negative (ER-) patients (p < 0.0001). In contrast, PR+/ER+ patients had higher scores when compared with PR-/ER- patients, and ER+ patients had higher when compared with ER- patients in the stromal scores ( Supplementary Figure S2D , p < 0.0001). Thus, stromal and immune scores were significantly associated with BC subtypes and hormone receptor status.

To identify potential OS correlations with immune scores and/or stromal scores, we divided our cohort into top and bottom halves (high vs. low score groups) based on their scores. Kaplan-Meier survival curves showed that median OS in the low score group was longer when compared with the high score group when based on immune scores ( Supplementary Figure S2E , p = 0.01). Consistently, patients with lower stromal scores had longer median OS when compared with patients with higher stromal scores ( Supplementary Figure S2F , p = 0.85), although statistical significance was not observed.

To determine global gene expression profile correlations with immune scores and/or stromal scores, we compared Affymetrix microarray data in 1,069 BC patients. Heatmaps in Figure 1 showed distinct gene expression profiles of cases belong to immune scores/stromal scores groups. Based on immune scores, 943 genes were upregulated, and 71 genes downregulated in the high score group than the low score group ( Figure 1A , fold change > 1.5, p < 0.05). Similarly, 1,011 genes were upregulated, and 50 genes were downregulated in the high score group ( Figure 1B , fold change > 1.5, p < 0.05). Moreover, Venn diagrams ( Figures 1C, D ) showed that 498 genes were upregulated in the high-score group, while two genes were downregulated. We performed subsequent analyses by focusing on all DEGs obtained based on comparisons of immune and stromal scores. To determine potential DEGs functions, we performed functional enrichment analysis on 1,574 DEGs. Top Gene Ontology (GO) terms included immune system process, immune response, extracellular matrix, signalling receptor binding, and integrin binding ( Figures 1E–G ).

To explore individual DEG correlations with OS, we performed Kaplan-Meier survival curve analysis. In total, 421 DEGs out of 1,574 significantly predicted OS in the log-rank test (p < 0.05, selected genes are shown in Supplementary Figure S3 ).

To better understand interactions between prognostic value DEGs, we examined protein-protein interaction (PPI) networks in STRING. The network consisted of eight modules, which included 218 nodes and 704 edges. We selected the top three important modules for further analysis ( Supplementary Figure S4 ). For descriptive convenience, we termed these modules MCODE1, MCODE2, and MCODE3 modules, respectively. In MCODE1 ( Supplementary Figure S4A ), ACKR3, CXCR3, and CCR5 had higher degree values. In MCODE2 ( Supplementary Figure S4B ), several immune response key genes occupied the module center and included HLA-DRB5, HLA-DRB1, CD247, and LCK. In MCODE3 ( Supplementary Figure S4C ), IL2RG, CD8B, and CD8A were significant nodes, as they had the most connections with other module members.

Consistent with PPI network analysis, functional enrichment analysis of these genes also identified strong associations with immune responses. Top GO terms included extracellular region and extracellular space ( Supplementary Figure S5A ), immune response ( Supplementary Figure S5B ), and antigen binding and signalling receptor binding ( Supplementary Figure S5C ). Additionally, all pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis ( Supplementary Figure S5D ) were associated with immune responses.

To determine if genes identified by TCGA had prognostic significance in other BC patients, we downloaded and analyzed gene expression data from 435 BC patients (GSE42568, GSE88770, GSE48390, and GSE162228) from the GEO database. Interestingly, 15 genes were significantly and prognostically related to the validation set ( Supplementary Figure S6 , p < 0.05); NPY1R, CELSR2, STC2, SCUBE2, GIMAP2, HLA-DPB1, TFF1, CXCL14, KLRB1, BIRC3, IL18, PSMB8, CD1C, TNFAIP8, and IRF1.

Subsequently, we performed least absolute shrinkage and selection operator (LASSO) Cox regression analysis to select highly relevant genes from these 15 genes. Finally, 14 were identified as related to the TME in BC, and optimal values of the penalty parameter were determined by 10-fold cross-validation ( Figures 2A, B ). We then constructed a prognostic model based on these genes, with the risk score of each sample from the training cohort calculated according to this model. Based on median risk score, BC samples from the training cohort were divided into high- and low-risk groups. To assess the OS in these groups, Kaplan-Meier curves were generated and showed that OS in the high-risk group was worse than that in the low-risk group, indicating the validity of the risk score prediction ( Figure 2C , p < 0.0001). Additionally, the expression of the TME-related genes, survival status, and survival time distribution for patients according to risk scores are shown in Figure 2D . In terms of model diagnosis, the AUC of the time-dependent receiver operating characteristic (ROC) curves were 0.69 for 1-year survival, 0.74 for 3-year survival, and 0.74 for 5-year survival, respectively, suggesting acceptable stability of the risk model ( Figure 2E ). In addition, to explore if BC subtypes affect survival, we grouped patients according to subtypes and subsequently performed survival analyses. Clearly, no differences in survival due to subtypes were observed, suggesting that the BC subtype did not affect survival ( Supplementary Figure S7 , p = 0.26). Together, our risk model, constructed from TME-related genes, appeared to accurately predict prognosis in BC patients.

Next, to identify hub genes, we identified interactions between genes in the TME model by constructing a PPI network in STRING. The network included 13 nodes and six edges. PSMB8 and BIRC3 had the maximum neighboring genes and were identified as hub genes. The Kaplan-Meier analysis showed both were the prognostic indicators, and its high expression favored the prognosis ( Supplementary Figure S2 , p < 0.05). To verify this phenomenon still exists in the human body, we used immunohistochemistry to compare hub protein expression and identified high PSMB8 and BIRC3 expression trends in BC epithelial cells when compared with paracancerous cells ( Figure 3 ).

Univariate Cox regression analysis showed that risk score could predict the prognosis of BC patients ( Figure 4A , p < 0.0001). In the multivariable Cox regression analysis, risk score remained statistically significant ( Figure 4B , p < 0.0001), indicating our risk model was an independent prognostic factor for BC. Additionally, BC patients in the training cohort were regrouped into subgroups based on age (< 50 and ≥ 50 years old), and TNM stage (stage I, stage II, stage III, and stage IV). Regardless of subgroups, low-risk group patients still showed significantly longer survival ( Figures 4C, D , p < 0.05), which indicated excellent risk model independence.

To create a quantitative method to predict OS, we integrated the risk score and independent clinicopathological prognostic factors, including age and TNM stage, to construct a nomogram ( Figure 5A ).

To evaluate its prognostic value, we compared the concordance index (C-index) of the nomogram with TNM stage, and as shown in Table 1 , the nomogram improved the prediction accuracy for BC. We compared predicted 3- and 5-year survival probabilities with actual probabilities and observed the calibration curve showed good concordance between these probabilities, thereby reflecting high nomogram accuracy and dependability ( Figure 5B ). Taken together, the nomogram, which integrated risk score, showed good performance and applicability, and has potential as a clinical tool to predict prognosis in BC patients.

Relationships between prognostic risk score and clinical characteristics were further investigated in the training cohort. Age, T category, M category and TNM stage were significantly related to risk score, whereas gender and N category were not ( Figure 6A , p < 0.05). As observed Figure 6B , patients with HER2-enriched had the highest risk score, followed by Basal-like, Luminal B, and Luminal A subtypes, while Normal-like patients had the lowest scores (p < 0.0001). Association analysis with hormone receptor status showed that patients with PR+/ER+ had lower risk score when compared with PR-/ER- patients, and ER+ patients had lower risk score when compared with ER- patients ( Figure 6C , p < 0.0001).

To better visualize the clinicopathological features in individual patients and assess correlations with survival, we used an alluvial diagram which showed that risk categories in the prediction model accurately predicted patient survival ( Figure 6D ).

Association between the risk model and immune cell infiltration was assessed using several immune infiltration approaches. ESTIMATE algorithm data showed that immune, stromal, and ESTIMATE scores in the high-risk BC patient group were lower when compared with BC patients in the low-risk group ( Figure 7A , p < 0.0001). The TIMER algorithm showed that B cell, neutrophil, CD4 T cell, dendritic cell (DC), and CD8 T cell abundance, but not macrophage, was statistically higher in the low-risk group when compared with the high-risk group ( Figure 7B , p < 0.0001). Moreover, MCPcounter algorithm results showed that T cells, CD8 T cells, cytotoxic lymphocytes, B lineage cells, natural killer (NK) cells, monocytic lineage cells, myeloid DCs, neutrophils, endothelial cells, and fibroblasts were highly infiltrated in the low-risk group ( Figure 7C , p < 0.01). Thus, our risk model correlated well with different immune microenvironment components.

Given the significant correlation of our risk model with the BC immune microenvironment, we next examined relationships between the risk model and immune cell subtype infiltration using Pearson’s algorithm. As shown in Figure 7D , correlation values for B cells, CD4 T cells, CD8 T cells, DCs, and neutrophils with risk scores were −0.35, −0.48, −0.49, −0.43, and −0.39, respectively. As expected, immune cell infiltration levels were significantly and positively correlated with prognosis ( Figure 7D , p < 0.0001).

To further confirm model practicability and reliability, it was verified using a validation cohort. Risk scores, survival status, and gene expression are shown in Figure 8A . As expected, significant differences in OS were identified between groups, with longer OS in the low-risk group ( Figure 8B , p < 0.0001). Furthermore, relationships between risk score and the BC immune microenvironment were confirmed in the validation cohort. From ESTIMATE, TIMER, and MCPcounter analysis, the low-risk group was significantly associated with high immune cell infiltration levels From ESTIMATE analysis, the low-risk group was significantly associated with high stromal, immune, and ESTIMATE scores ( Figure 8C , p < 0.0001). In TIMER analysis, the abundance of the five aforementioned immune cell types, except macrophages, was statistically different between groups ( Figure 8D , p < 0.0001), and immune cell abundance (all types) was significantly higher in the low-risk group than the high-risk group. The MCPcounter algorithm showed that T cells, cytotoxic lymphocytes, B lineage, monocytic lineage cells, myeloid DCs, endothelial cells, neutrophils, and fibroblasts were in a high infiltration state in the low-risk group ( Figure 8E , p < 0.05). Therefore, our TME-related gene risk model was associated with BC prognosis and the immune microenvironment.

As neoadjuvant and adjuvant chemotherapies are reportedly related to immune infiltration (37), we evaluated if chemotherapy influenced BC prognosis. According to the NCCN Guidelines in Oncology, anthracycline + cyclophosphamide (AC), AC followed by taxane (AC-T), and taxane + cyclophosphamide (TC) are major chemotherapy regimens. The OS advantage was observed in the low-risk group, regardless of whether they received chemotherapy or not. And whether in high-risk group or low-risk group, patients who received chemotherapy had a better prognosis ( Figure 9A , p < 0.0001). In the low-risk group, the OS advantage was evident in patients who received TC and AC-T chemotherapy regimens when compared with those who received no chemotherapy ( Figure 9B , p < 0.05). In contrast, the chemotherapy benefits in the high-risk group were observed for AC, TC, and AC-T chemotherapy regimens ( Figure 9C , p < 0.05). More importantly, subgroup interaction evaluations suggested that better chemotherapy outcomes were achieved in low-risk patients regardless of the chemotherapy regimen ( Figure 9D , p < 0.05).

We also explored if the A-regimen was an indispensable chemotherapy agent in the low-risk group. As shown in Figure 9E , no significant differences in prognosis outcomes for patients treated with the A-regimen were identified, regardless of low- or high-risk (p > 0.05). Further subgroup analysis showed no significant differences in prognosis outcomes in low-risk patients who received the A-regimen when compared with those who did not ( Figure 9F , p > 0.05). These observations suggested that the low-risk group selected by this prediction model has the opportunity to exempt the A-containing adjuvant chemotherapy regimen.

Immune checkpoint blockade with immunotherapies, including CTLA-4, CD28, and CD274 are promising treatment approaches for several malignancies (38). However, the bottleneck problem of immune checkpoint inhibitors (ICI) in the treatment of eBC is the lack of precise biomarkers identifying populations who may benefit from these therapeutics. In our study, we determined the expression levels of several key immune checkpoint regulators and inflammatory mediators to provide reference biomarker candidates for precision immunotherapy in early drug-resistant patients. As presented in Figure 10A , CD274, CD28, and CTLA-4 expression levels were significantly higher in the low-risk group (p < 0.0001). The Pearson algorithm was used to analyze correlations between immune checkpoints and our risk model. Correlation values of CTLA-4, CD28, CD274 and risk score were -0.37, -0.43 and -0.33, respectively ( Figure 10B , p < 0.0001). Additionally, other immunomodulators or inflammatory mediators were increased in the low-risk group ( Figure 10C , p < 0.0001). A previous study reported that HLA affected ICI efficacy (39), therefore we analyzed correlations between HLA family expression and our model, and showed this expression was significantly higher in the low-risk group when compared with the high-risk group ( Figure 10D , p < 0.0001). We next explored ICI therapy responses, represented by the CTLA-4/PD1 inhibitors, by using the immunophenotype score (IPS), and showed that the IPS was slightly higher than that of the low-risk group in the patients treated with CTLA-4 and PD1 inhibitors ( Figure 10E , p < 0.05). Overall, these results suggested that our model predicted the immunotherapy benefits for patients and may be a more effective biomarker to predict the efficacy of immunotherapy.

We further analyzed DEGs between low- and high-risk groups in TCGA. In total, 396 DEGs (7 upregulated and 389 downregulated genes, FDR p-value < 0.05) were identified in the high-risk group when compared with the low-risk group. Of these, SLC7A5, PRAME, CRABP1, CBX2, CA9, CALML5, and CD24 were significantly overexpressed in the high-risk group ( Supplementary Figures S8A, B , FDR p-value < 0.05, fold change > 1.5). Furthermore, KEGG analysis showed that genes in the high-risk group were mainly involved in environmental information processing, human diseases, and organismal systems ( Supplementary Figure S8C ). From GO enrichment analysis, these genes in the high-risk group were mainly involved in extracellular matrix, vesicle, immune response, and antigen binding ( Supplementary Figures S8D–F ).

As shown in Figure 11A , BC patients in the high-risk group had a higher TMB than those in the low-risk group (p < 0.05). As suggested from previous studies, a high TMB leads to a poor prognosis in many cancers (40), consistent with our data. In correlation analysis between risk score and TMB, we identified a significant positive correlation ( Figure 11B , p < 0.05). Further survival analysis indicated that the low-TMB group showed a significant survival benefit ( Figure 11C , p < 0.05). Given the synergistic effect of TMB and the risk score, their effect on prognostic stratification was evaluated. As indicated from the results, TMB status did not interference the predictive ability of the risk score. Survival difference of the risk score subtypes was significant in both high- and low-TMB groups, and the subgroup with low risk-score and low TMB showed a better survival benefit, while the high-risk score and high TMB subgroup had a lower survival probability ( Figure 11D , p < 0.001). Combined, risk score may act as a prognostic BC indicator, which is independent of TMB and can effectively predict TMB and treatment sensitivity.

We used the TIMER database to explore the relationships between hub genes expression (PSMB8 and BIRC3) and the level of infiltrating lymphocytes. Upregulated PSMB8/BIRC3 expression was associated with increased B cell, CD8+ T cell, macrophage, neutrophil, DC, and other infiltrating lymphocyte infiltration ( Figures 12A, B , p < 0.05). Next, using the TISIDB database, we found that upregulated PSMB8/BIRC3 ( Figures 12C, D ) expression was associated with increased expression of immunostimulatory molecules, immunosuppressive molecules, MHC molecules, chemokines, and chemokine receptors, which provides important information for predicting potential therapeutic targets. Finally, we used GSCA Lite online tool to analyze the relationship between the expression of the hub genes and sensitivity to current immune or targeted therapies for BC ( Figure 12E ). PSMB8 expression levels were negatively correlated with sensitivity to many BC-targeting or chemotherapeutic drugs, including clofarabine and gemcitabine, and were positively correlated with abiraterone. BIRC3 expression levels were positively correlated with axitinib sensitivity and negatively correlated with dasatinib sensitivity. Thus, hub genes could function as new targets for predicting drug sensitivity and developing multi-targeted combined therapy for BC.