The gene expression data were obtained from the Gene Expression Omnibus (GEO) database. The dataset GSE194040, derived from the I-SPY2 Trial (“I-SPY2-990” subproject), includes pretreatment transcriptomic and proteomic/phosphoproteomic profiles from patients with early-stage invasive breast cancer undergoing neoadjuvant therapy. A total of 988 pretreatment breast cancer samples with comprehensive clinical annotations, including treatment arm, pathological complete response (pCR) status, molecular subtype, and MammaPrint (MP) classification, were available. To eliminate the confounding effects of anti-HER2 targeted therapy, we included only the 743 HER2-negative patients who received anthracycline- and taxane-based neoadjuvant chemotherapy, comprising 537 non-pCR and 206 pCR cases. For external validation, nine additional GEO datasets containing pretreatment transcriptomic profiles of HER2-negative breast cancer patients treated with anthracycline- and taxane-based neoadjuvant regimens (GSE20194, GSE20271, GSE22093, GSE23988, GSE25066, GSE32646, GSE34138, GSE41998, and GSE163882) were retrieved, encompassing a total of 1,642 patients (1,281 non-pCR and 361 pCR). To reduce platform related bias when integrating the external validation cohorts, expression matrices from the nine GEO validation datasets were merged using the intersection of shared genes. Batch effects across datasets were then corrected using the ComBat function in the sva package, with the dataset identifier specified as the batch variable. Principal component analysis was performed before and after batch correction to visually assess whether inter cohort technical variation was mitigated. The ANN model was trained on the independent training cohort, whereas the merged and batch corrected validation cohort served for external evaluation, aiming to assess generalizability under cross platform conditions after mitigation of batch effects. Moreover, a comprehensive list of immune-related genes was downloaded from the ImmPort database (https://www.immport.org) for subsequent analyses.

To identify the transcriptomic differences associated with chemotherapy response across hormone receptor (HR) subtypes, the 743 HER2-negative breast cancer patients from the GSE194040 dataset were stratified into 364 HR-negative and 379 HR-positive cases. Within each HR subgroup, differential expression analysis between pCR and non-pCR samples was performed using the limma R package. Gene expression matrices were normalized with the normalizeBetweenArrays function, and differentially expressed genes (DEGs) were identified based on |log₂FC| > 1 and adjusted p < 0.05 (Benjamini–Hochberg correction) (19). The thresholds of |log₂ fold change| greater than 1 and adjusted p value less than 0.05 were selected to ensure that identified differentially expressed genes represented biologically meaningful expression changes while maintaining adequate statistical stringency. Applying the same criteria across hormone receptor positive and hormone receptor negative subgroups allowed for consistent comparison of transcriptional differences associated with treatment response, despite intrinsic biological heterogeneity between subtypes. Heatmaps and volcano plots were generated using pheatmap and ggplot2 to visualize the expression patterns and statistical significance of DEGs (20). Subsequently, the DEGs from the HR-positive and HR-negative cohorts were compared, and genes that were significantly dysregulated in both groups with concordant regulation direction were extracted as overlapping DEGs for downstream analyses. All computations and visualizations were conducted in R (version 4.5.0).

To further identify gene co-expression modules associated with chemotherapy response, we performed a weighted gene co-expression network analysis (WGCNA) based on the normalized transcriptomic profiles of 743 HER2-negative breast cancer patients (pCR vs. non-pCR) from the GSE194040 dataset (21). Genes with low variance (standard deviation < 0.3) were excluded to improve network robustness. Outlier samples were detected and removed through hierarchical clustering. The optimal soft-thresholding power (β) was determined using the pickSoftThreshold function to ensure a scale-free network topology (R² > 0.9). An adjacency matrix was constructed and transformed into a topological overlap matrix (TOM) to quantify gene interconnectedness, followed by hierarchical clustering to identify gene modules using the dynamicTreeCut algorithm (minimum module size = 50). Highly similar modules were merged at a cut height of 0.3 to obtain the final module structure. The relationships between module eigengenes and clinical traits (pCR and non-pCR status) were calculated using Pearson correlation, and a module–trait heatmap was generated to visualize the strength and significance of associations. Module significance and gene significance were computed to identify modules most strongly correlated with therapeutic response (22). Modules were prioritized for downstream analyses based on the strength and direction of their correlation with pCR status, statistical significance, and overall module significance, with preference given to modules showing the most robust and biologically plausible associations with treatment response.

To obtain robust candidate genes associated with neoadjuvant chemotherapy response, the genes from the module most significantly correlated with pCR in the WGCNA (23), the directionally concordant DEGs shared by HR-positive and HR-negative subgroups, and the immune-related genes from the ImmPort database were intersected to generate a key gene set for feature selection. Three complementary machine learning algorithms were then applied. First, least absolute shrinkage and selection operator (LASSO) regression was performed using the glmnet package in R (24), with tenfold cross-validation to identify the optimal penalty parameter (λ) that minimized deviance, and genes with nonzero coefficients at λ_min were retained. Second, a random forest (RF) classifier was constructed using the randomForest package with 500 decision trees to evaluate the importance of each gene based on the mean decrease in Gini index, and genes with importance scores above the optimal cutoff were selected. Third, a support vector machine–recursive feature elimination (SVM-RFE) approach was implemented via the e1071 package and an in-house script to iteratively eliminate less informative features and determine the optimal subset of genes with the lowest cross-validation error. Genes consistently identified by LASSO, RF, and SVM-RFE were considered as key feature genes for subsequent modeling and validation.

To develop and validate a predictive model for chemotherapy response, the intersecting genes jointly identified by the LASSO, RF, and SVM-RFE algorithms were used to construct an artificial neural network (ANN). The normalized expression profiles of 743 HER2-negative breast cancer patients from the GSE194040 dataset (training cohort) served as the training data, while an independent validation cohort comprising 1,642 HER2-negative breast cancer patients from nine GEO datasets (GSE20194, GSE20271, GSE22093, GSE23988, GSE25066, GSE32646, GSE34138, GSE41998, and GSE163882) was used for external validation. A feed-forward ANN model was implemented using the neuralnet package in R, with two output nodes (pCR and non-pCR) and a single hidden layer containing three neurons. Model training was performed with a maximum of 1,000,000 iterations to ensure convergence. The network structure and connection weights were visualized using the NeuralNetTools package. The trained model was then applied to the external validation dataset to predict pCR probabilities, and model performance was assessed using receiver operating characteristic (ROC) analysis via the pROC package. The area under the curve (AUC) and corresponding 95% confidence interval (calculated using the bootstrap method) were used to evaluate the discriminative ability and generalizability of the ANN model in predicting pCR among HER2-negative breast cancer patients (25).

The architecture of the artificial neural network was intentionally kept simple to reduce the risk of overfitting and to improve model interpretability. The number of hidden neurons was selected empirically through preliminary testing, in which multiple network configurations with different hidden layer sizes were evaluated for stability and predictive performance. A single hidden layer with three neurons provided a favorable balance between model complexity and generalization ability and did not yield substantial performance improvement when additional neurons were introduced. The maximum number of training iterations was set to 1,000,000 to ensure convergence of the network weights during optimization. Model training typically converged well before reaching the maximum iteration limit, which was used as an upper bound rather than a fixed stopping point. This strategy ensured stable model training without forcing unnecessary complexity. These hyperparameter choices were guided by a combination of empirical optimization within the training cohort and prior experience from published studies applying neural network models to transcriptomic data with limited feature dimensions.

To investigate the biological relevance of the identified feature genes, we first validated their expression differences between pCR and non-pCR groups in both the training and external validation cohorts. Normalized expression matrices were subset to the intersected feature genes, and two-sided Wilcoxon rank-sum tests were applied to compare gene expression levels, with results visualized as boxplots using the ggpubr package. Subsequently, immune infiltration analysis was performed only in the training cohort. The proportions of 22 immune cell types were estimated using the CIBERSORT algorithm (1,000 permutations) based on the LM22 leukocyte gene signature matrix (26). Differences in immune-cell fractions between pCR and non-pCR samples were evaluated using Wilcoxon tests, and boxplots were generated to display significantly altered immune populations. Finally, Spearman correlation analysis was conducted to examine associations between model gene expression and immune-cell infiltration levels, and the correlation network was visualized using the linkET package. To account for multiple hypothesis testing when comparing the relative proportions of 22 immune cell types, p values were further adjusted using the Benjamini–Hochberg false discovery rate method. Both nominal and FDR adjusted p values were considered when interpreting immune infiltration differences between groups.

GSEA was performed in the training cohort to interrogate pathways associated with the five key genes. For each gene, samples were dichotomized into high- and low-expression groups by the cohort median; because all five genes were overexpressed in pCR, we focused on the high-expression groups. A ranked gene list was generated by the difference in mean expression between high and low groups, and enrichment was tested with clusterProfiler (using MSigDB KEGG collection c2.cp.kegg.Hs.symbols.gmt, https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) (27). Results were summarized by normalized enrichment score (NES) and significance (nominal p < 0.05; adjusted q values were reported when available), and the top enriched pathways (up to five per gene) were visualized with enrichplot (gseaplot2). Analyses used R with limma for preprocessing and org.Hs.eg.db for gene annotation. For enrichment analyses, adjusted p values calculated using the Benjamini–Hochberg method were used to define statistical significance unless otherwise specified. Nominal p values were reported only for descriptive purposes and were explicitly labeled as such in the corresponding Results.

To further explore the regulatory mechanisms and therapeutic relevance of the key genes, drug enrichment, ceRNA, and transcription factor (TF) network analyses were performed. Drug enrichment was conducted using the DSigDB database (https://dsigdb.tanlab.org/) with the clusterProfiler package, applying the enricher function based on the DSigDB_All_detailed dataset. Significantly enriched drugs were defined as those with p < 0.05 and adjusted p < 0.05, and visualized using enrichplot and ggplot2. For post-transcriptional regulation, a ceRNA network was constructed by integrating miRNA–mRNA interactions predicted by miRanda, miRDB, TargetScan, and miRWalk, retaining miRNAs identified in at least four databases. Corresponding lncRNA–miRNA interactions were obtained from spongeScan, and the merged network was visualized in Cytoscape (v3.10.1). Upstream transcriptional regulation was investigated using the TRRUST v2 database (https://www.grnpedia.org/trrust/), where curated TF–target relationships were retrieved, and upstream TFs regulating the key mRNAs were identified to construct the TF–gene regulatory network for visualization and interpretation.

To experimentally validate the expression of key predictive genes, pretreatment tumor tissue samples were collected from 60 HER2-negative breast cancer patients, including 30 patients who achieved a pathological complete response (pCR) and 30 patients who did not achieve pCR (non-pCR) following anthracycline- and taxane-based neoadjuvant chemotherapy (28). Total RNA was extracted from tumor tissues and reverse-transcribed into complementary DNA. Quantitative real-time PCR (qRT-PCR) was performed to measure the mRNA expression levels of five immune-related genes (CCL2, CXCL10, CXCL13, HLA-E, and IGKV1D-8) using gene-specific primers. The primers used were as follows:

CCL2, forward 5′-CATCTCCTACACCCCACGAAG-3′ and reverse 5′-GGGTTGGCACAGAAACGTC-3′;

CXCL10, forward 5′-GCTCTAGAATTGCTGCCTTATCTTTCT-3′ and reverse 5′-CTCAACACGTGGGCAGGATTA-3′;

CXCL13, forward 5′-CTCAACGTCCCATCTACTTGC-3′ and reverse 5′-TCTTCAGGGTGTGAGCTTTCC-3′;

HLA-E, forward 5′-CCTGGAGGAGATACGGGAGAC-3′ and reverse 5′-GGACAGCCCTCCATGGT-3′;

IGKV1D-8, forward 5′-TCCTGGGGCTCTCATGCT-3′ and reverse 5′-GGATGGTGGGAAGATGGA-3′.

ACTB (β-actin) was amplified with forward primer 5′-AAGGAGCCCCACGAGAAAAAT-3′ and reverse primer 5′-ACCGAACTTGCATTGATTCCAG-3′.

All reactions were run in triplicate. Expression levels were normalized to ACTB and calculated using the 2^−ΔΔCt method. Differences in relative mRNA expression between the pCR (n = 30) and non-pCR (n = 30) groups were compared with two-sided Wilcoxon rank-sum tests.

All statistical analyses were conducted using R software (version 4.5.0). Differences in gene expression and immune-cell infiltration between groups were evaluated using two-sided Wilcoxon rank-sum tests. The relationships between module eigengenes and clinical traits in the WGCNA were assessed using Pearson correlation analysis. Gene set enrichment analysis (GSEA) and drug enrichment analysis were performed using the clusterProfiler package, with significance defined as nominal p < 0.05 and adjusted p < 0.05. Correlations between key gene expression levels and immune-cell fractions were examined using Spearman’s rank correlation. All tests were two-tailed, and a p < 0.05 was considered statistically significant.

To evaluate and mitigate potential batch effects among the nine external GEO validation cohorts, we merged these datasets based on the common gene set and applied ComBat batch correction using the cohort identifier as the batch factor. As shown in Figures 1A, B, samples exhibited clear cohort driven separation before correction, whereas the inter cohort clustering pattern was substantially attenuated after correction, indicating reduced platform related bias. To explore transcriptional differences associated with neoadjuvant chemotherapy response, differential expression analysis was performed separately in HR-negative and HR-positive HER2-negative breast cancer subgroups from the GSE194040 dataset. In HR- positive tumors, a total of 669 DEGs were identified (|log₂FC| > 1, adjusted p < 0.05), including 362 upregulated and 307 downregulated genes, clearly separating pCR from non-pCR samples (Figures 2A, B). Similarly, in HR-negative tumors, 89 DEGs were identified (86 upregulated and 3 downregulated), demonstrating distinct expression patterns between response groups (Figures 2C, D). Comparison of the two DEG sets revealed 77 overlapping genes with concordant regulation direction (Figure 2E), indicating a shared transcriptional signature related to chemotherapy sensitivity across HR subtypes. The difference in the number of differentially expressed genes between hormone receptor positive and hormone receptor negative tumors likely reflects inherent biological differences in transcriptional variability and immune activation between these subtypes rather than analytical bias introduced by the applied thresholds.

WGCNA was conducted using expression data from 743 HER2-negative patients to identify gene co-expression modules associated with pCR. The scale-free topology criterion (R² > 0.9) was achieved at β = 4 (Figures 3A, B). Hierarchical clustering identified four major co-expression modules (blue, brown, yellow, and turquoise; Figure 3C). Among them, the blue module showed the strongest positive correlation with pCR (r = 0.31, p < 0.01) and negative correlation with non-pCR (r = −0.31; Figure 2D). Among all identified modules, the blue module showed the strongest positive correlation with pathologic complete response, with a correlation coefficient of 0.31 and a statistically significant p value. This module contained a total of 1004 genes and exhibited a higher module significance and gene significance for pCR compared with other modules, supporting its prioritization for downstream analyses.

To identify robust immune biomarkers associated with chemotherapy response, genes from the WGCNA blue module, the ImmPort immune gene list, and the common DEGs with consistent regulation direction between HR subtypes were intersected, yielding 43 overlapping genes (Figure 4). These genes were subsequently subjected to multi-algorithmic feature selection. LASSO regression identified 17 genes with nonzero coefficients at the optimal λ (Figures 5A, B), while SVM-RFE achieved the highest classification accuracy (0.735) and lowest error (0.265) with 25 features (Figures 5C, D). Random forest analysis ranked gene importance and selected the top ten feature genes (Figures 5E, F), which were subsequently used for model construction and validation.

The intersection of LASSO, RF, and SVM-RFE analyses identified five shared key genes—CCL2, CXCL10, CXCL13, HLA-E, and IGKV1D-8—which were subsequently used to construct an artificial neural network (ANN) model for predicting chemotherapy response (Figure 6A). The ANN was developed using the GSE194040 cohort (n = 743) as the training set and externally validated in nine independent GEO datasets (n = 1,642). The model comprised five input nodes corresponding to the key genes, a single hidden layer with three neurons, and two output nodes representing pCR and non-pCR outcomes (Figures 6B, C). The ANN exhibited strong discriminative performance, achieving an AUC of 0.858 (95% CI: 0.829-0.888) in the training cohort and 0.773 (95% CI: 0.735–0.808) in the validation cohort (Figures 6D, E), indicating its robust and generalizable predictive capacity across independent datasets.

In both the training and external validation cohorts, all five key genes were significantly upregulated in the pCR group compared with the non-pCR group (Figures 7A, B). CIBERSORT-based immune infiltration analysis revealed distinct immune cell distributions between response groups (Figure 7C). The pCR group showed higher proportions of T cells CD8, T cells CD4 memory activated, T cells follicular helper, M1 macrophages, and activated dendritic cells, reflecting an inflammatory and immune-stimulating microenvironment. In contrast, non-pCR tumors exhibited increased infiltration of naïve B cells, Plasma cell, T cells CD4 memory resting, T cells regulatory (Tregs), M0 and M2 macrophages, and Mast cells resting, suggesting an immunosuppressive state. Correlation network analysis (Figure 7D) further revealed that CCL2, CXCL10, CXCL13, HLA-E, and IGKV1D-8 were positively correlated with multiple effector and antigen-presenting immune subsets, including T cells CD8, T cells CD4 memory activated, T cells follicular helper, Macrophages M1, and Dendritic cells activated, indicating their close association with immune activation and cytotoxic responses. Conversely, these genes were negatively correlated with B cells naïve, Plasma cells, T cells CD4 memory resting, T cells regulatory (Tregs), Macrophages M0, Macrophages M2, and Mast cells resting, which are characteristic of an immunosuppressive microenvironment. Collectively, these findings suggest that high expression of the five key genes is linked to a pro-inflammatory, antitumor immune milieu favorable for chemotherapy responsiveness.

GSEA revealed that the high‐expression groups of the five key genes were significantly enriched in immune‐related pathways (Figures 8A–E). Enriched pathways were considered statistically significant based on adjusted p values, whereas pathways reported using nominal p values are explicitly indicated and interpreted with caution. Specifically, the high CCL2 group was enriched in allograft rejection, cell adhesion molecules (CAMs), chemokine signaling pathway, cytokine–cytokine receptor interaction, and type I diabetes mellitus. The high CXCL10 group showed enrichment in antigen processing and presentation, cytokine–cytokine receptor interaction, hematopoietic cell lineage, intestinal immune network for IgA production, and Leishmania infection. The high CXCL13 group was enriched in allograft rejection, antigen processing and presentation, cytokine–cytokine receptor interaction, intestinal immune network for IgA production, and Leishmania infection. The high HLA-E group exhibited enrichment in antigen processing and presentation, cytokine–cytokine receptor interaction, Leishmania infection, natural killer cell–mediated cytotoxicity, and systemic lupus erythematosus. Finally, the high IGKV1D-8 group was enriched in allograft rejection, antigen processing and presentation, asthma, primary immunodeficiency, and type I diabetes mellitus.

A comprehensive ceRNA network was constructed to explore the post-transcriptional regulatory relationships of the key genes (Figure 9A). The network included 3 mRNAs (CCL2, CXCL10, CXCL13), 26 miRNAs, and 42 lncRNAs, forming a multilayered regulatory architecture. Among them, CCL2 was potentially regulated by multiple miRNAs such as hsa-miR-125a-5p, hsa-miR-30b-3p, and hsa-miR-195-5p, which were further modulated by lncRNAs including LINC00969, AC011284.3, and HCP5. CXCL10 was predicted to be targeted by hsa-miR-29b-3p, hsa-miR-34a-5p, and hsa-miR-340-5p, while CXCL13 was associated with hsa-miR-101-3p, hsa-miR-155-5p, and hsa-miR-450a-5p, all of which were interconnected with lncRNAs such as MIR155HG, LINC00511, and AC006449.2. These results delineate a complex ceRNA regulatory landscape potentially influencing gene expression and immune activity. Transcriptional regulation analysis based on the TRRUST v2 database identified several upstream transcription factors (TFs) strongly associated with the key genes (Figure 9B). CCL2 and CXCL10 were mainly regulated by NF-κB1, RELA, STAT1, STAT3, and IRF1, whereas CXCL13 was primarily controlled by IRF4 and SPI1. These TFs are well known for their roles in inflammatory and interferon-mediated signaling, suggesting transcriptional activation of these chemokine genes under immune-stimulating conditions. The drug enrichment analysis performed using the DSigDB database revealed a series of compounds potentially targeting the key genes (Figure 9C). The top significantly enriched molecules included pioglitazone, ibuprofen, roflumilast, rutin, methylprednisolone, beclomethasone, allopurinol, and acteoside. Several of these agents, such as pioglitazone (a PPAR-γ agonist) and roflumilast (a phosphodiesterase-4 inhibitor), have recognized immunomodulatory or anti-inflammatory effects, whereas ibuprofen and rutin were associated with cytokine pathway regulation. The enrichment of these compounds suggests that modulation of inflammatory and chemokine signaling might enhance chemotherapy sensitivity in HER2-negative breast cancer.

To experimentally validate the five identified immune biomarkers, we measured their mRNA expression levels in pretreatment breast tumor tissues from patients with and without pCR. qRT-PCR analysis demonstrated that CCL2, CXCL10, CXCL13, HLA-E, and IGKV1D-8 were all significantly upregulated in pCR cases compared to non-pCR cases (p < 0.0001; Figures 10A–E). These results were consistent with prior transcriptomic findings from the GSE194040 cohort and support the predictive value of these immune genes in determining neoadjuvant chemotherapy response.