We extracted and analyzed the sc-RNAseq data involved in this study from the Gene Expression Omnibus (GEO) database under the accession number GSE161529, encompassing four samples with BRCA1-mutant (BRCA1-MT1-4) and four samples with BRCA1 wild (BRCA1-WT1-4). Seurat package was used to generate the Seurat objects containing scRNA-seq gene expression matrices for main cell types. Low-quality cells were removed from each sample according to nFeature, nUMI per cell and, mitochondria content ( Supplementary Figure S1 ). Further, the expression matrices underwent normalization and scaling via the NormalizeData and ScaleData functions in the Seurat package. According to the top 2000 highly variable genes and 30 principal components, reduce dimension was employed via the Uniform Manifold Approximation and Projection (UMAP) algorithm. Harmony algorithm (16) was conducted to counteract the batch effect. The marker genes for each cell cluster were identified through the FindAllmarkers function, with the cutoff criteria as log₂FC > 0.25 and false discovery rate (FDR) < 0.05. With the aid of well-known cell markers and singleR package (17), major cell populations were annotated. The clustering of each major cell type was subsequently re-clustered using the workflow described above. These new clusters were identified as “subclusters”, representing different phenotypes within each major cell type. Tissue enrichment was assessed using the ratio of observed to expected cell number (Ro/e) in each tissue type. Specifically, Ro/e > 1 indicates enrichment of the cell cluster in that tissue, whereas Ro/e < 1 suggests under-representation (18).

The filtered spatial transcriptome data was fetched from the website, https://zenodo.org/record/4739739/. Seurat package was used to conduct subsequent analysis. According to the top 3000 highly variable genes, the expression matrix of each slice was normalized and scaled utilizing the SCTransform function. Subsequently, PCA dimensionality reduction analysis was performed for building an SNN graph with the default parameter. Using the anchor-based integration pipeline within Seurat package, we established a mapping between the spatial RNA-seq and the scRNA-seq data, yielding the prediction scores of subcluster to every spot within each slice.

Gene ontology (GO) analysis was conducted via clusterProfiler package (19). We utilized gene set variation analysis (GSVA) to calculate the enrichment scores of hallmark pathways for each cell. The differential activating pathways between two cell subclusters were identified using the limma package (20). The significant terms and pathways were identified with the cutoff of adjusted p-value <0.05.

The R packages, monocle2, and monocle3 were applied to explore the trajectory of cell differentiation among the selected clusters (21). For monocle2, the DDRtree algorithm was employed to reduce the dimensions. The mutual nearest neighbor algorithm was used to eliminate the batch effect for monocle3.

We utilized the milo algorithm for identifying the differential abundance of cell neighborhoods in two distinct conditions (22). More specifically, a KNN graph based on scRNA-seq data was constructed, followed by the data split into cell neighborhoods and the differential abundance analysis of selected cell subclusters was performed. In this section, we configured the parameters with k=40 and d=30, following the author’s recommendation.

The copykat software was utilized to infer genome-wide copy number variations from the single-cell gene expression count matrix (23). Immune cells were employed as reference cells, where single aneuploid cells with copy number variations were deemed as tumor cells, and diploid cells were predicted as normal epithelial cells. The CytoTRACE (24) algorithm was employed to validate the tumor cell predictions derived from CopyKAT algorithm.

For quantifying the difference in transcription factors (TF) activities between distinct cell clusters, cluster-specific gene regulatory networks were built utilizing pySCENIC (25) package with default parameters. Count matrix of scRNA-seq data was extracted from Seurat object and converted to the Loom format file as input for downstream analysis. In particular, we first utilized GRNBoost2 to form a co-expression network, followed by discovering regulons for every TF via RcisTarget. The motif database of humans is accessible through the website https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/. The activity of each regulon was measured through the AUCell algorithm. Subsequently, binary values for each regulon were obtained using AUC score thresholds automatically determined in the process. For one specific subcluster, we regarded TFs that activated in over 20% of cells in at least one subcluster as significant.

To explore the distinct levels of tumor heterogeneity between the BRCA1-WT group and the BRCA1-MT group, we calculated the intratumoral and intertumoral heterogeneity scores as previously reported (26). Notably, malignant tumor cells based-CopyKAT were chosen and re-clustered using the method described above with default parameters. Extracting the top 50 principal components as features for calculating the heterogeneity score. To ensure consistency in data scaling, z-score standardization is used for all heterogeneity scores.

The hdWGCNA package (27), which was designed for analyzing high-dimensional scRNA-seq data, was deployed to construct a scale-free network at the single-cell level utilizing default parameters.

In our effort to identify the critical molecular marker in one cell subgroup, we developed a reliable analysis pipeline via employing multi-scale data, encompassing bulk and scRNA-seq data. Initially, we identified the marker genes of myCAF utilizing FindAllMarkers function, with the threshold as adjusted p-value<0.05 and log2FC>0.25. Following this, Differential module eigengene analysis was conducted to discover the myCAF and BRCA1-WT specific module. For each module, the top 200 hub genes were extracted. Upon the overlap genes of these myCAF marker genes and hub genes ( Supplementary Table S1 ), the differential expression analysis was employed between BRCA1-WT samples and BRCA1-MT samples within the METABRIC dataset. The substantially elevated genes in the BRCA1-WT group were regarded as the key genes ( Supplementary Table S2 ).

Thirteen ICI pre-therapy RNA-seq cohorts with response information (n = 829) were collected in this study ( Supplementary Table S3 ). The expression values of RNA-seq data for each cohort were converted to TPM values, followed by log2-transformed. To eliminate unnecessary interference, we utilized the Combat function in the sva package to integrate the four cohorts encompassing Braun 2020 (n = 172) (28), Mariathasan 2018 (n = 298) (29), Liu 2019 (n = 119) (30), and Pender 2021 (n = 72) (31). We randomly allocated the integrated cohort (n = 661) into an 80% training set and a 20% validation set. Model performance was streamlined via 5-fold cross-validation. Utilizing the tidymodels package, eight machine learning algorithms were implemented to derive the predictive models for ICI therapy response. These algorithms included Support Vector Machine (SVM), Extreme Gradient Boosting (XGBoost), Random Forest (RF), LightGBM, naive Bayes (NB), multilayer perceptron (MLP) neural network, K-Nearest Neighbour (KNN). In terms of the validation set, the model with the maximum receiver operating characteristic curve (ROC) was deemed as optimal.

After surgically proving, paired tumor (n=9) and para-carcinoma tissues (n=9) from breast cancer patients were collected from the First Affiliated Hospital of Zhengzhou University. Samples were obtained with informed consent. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (Ethics review number: 2024-KY-0549; Zhengzhou, China). All methods were conducted in accordance with the relevant guidelines and regulations. Subsequently, the protein expression was quantified via high-sensitivity label-free quantitative proteomics sequencing, as described in the Supplementary Methods .

Fisher’s exact test was implemented for categorical variables, whereas Student’s t-test, Wilcoxon rank-sum test, ANOVA test, and Kruskal-Wallis test were conducted for continuous variables. Spearman correlation was used for continuous versus continuous variables. Survival analysis was performed through Kaplan-Meier curves and log-rank tests. To correct for multiple tests, the p-values were adjusted to the FDR utilizing the Benjamini-Hochberg approach where appropriate. Unless otherwise indicated, all statistical tests were two-tailed. Significance levels are indicated by asterisks (*p <0.05; **p <0.01; ***p <0.001, ****p <0.0001). Statistical and bioinformatics analyses mentioned above were carried out with R software (version 4.2.1).

To understand gene-expression perturbations and generate a comprehensive map of the TME landscape of TNBC, we collected and analyzed eight specimens, comprising four samples with BRCA1-mutant (BRCA1-MT1-4) and four samples with BRCA1 wild (BRCA1-WT1-4) ( Figure 1A ). After filtering low-quality cells through a rigorous quality control pipeline, 32,386 cells remained for downstream analysis ( Supplementary Figures S1A, B ). Of these, 18,182 cells were derived from BRCA1-MT samples, while 14,204 cells were obtained from BRCA1-WT samples. The mean number of genes examined per cell was 1,331. Afterward, the standard workflow in Seurat was utilized to identify the subpopulations within scRNA-seq data, yielding seven unique cell subsets ( Figure 1B ). Leveraging the expression of well-characterized markers, we identified cell lineages, including Epithelial cells, immune cells (NK/T cells, B cells, and Myeloid cells), stromal cells (Fibroblasts, pericytes, and endothelial cells) ( Figure 1C , Supplementary Figure S2 ).

To gain further insight, we assess the difference in TME landscape between the BRCA1-WT group and the BRCA1-MT group, and the results indicated that the significantly elevated percent of immune and stromal cells were observed in the BRCA1-MT group, comparing to BRCA1-WT group with more epithelial cells ( Figure 1B ). Next, differential expression analysis for each cluster was executed to decode the transcriptional perturbations caused by BRCA1 mutation, totally yielding 2,199 differentially expressed genes (DEGs) ( Figure 1D ). Regarding epithelial cells, BRCA1-WT samples harbored more up-regulated DEGs than BRCA1-MT samples. Of note, upregulated genes in BRCA1-WT were mainly enriched in immune cells and Fibroblasts. Although BRCA1-WT tumors showed a higher proportion of epithelial cells, the functional disturbances reflected by DEGs were predominantly associated with immune and stromal compartments. This may suggest a compensatory or stress-induced transcriptional response in these less abundant cell types, highlighting the importance of integrating both compositional and functional perspectives. Altogether, our findings dedicated the cellular heterogeneity within TME caused by BRCA1 mutation, which was meticulously analyzed in the subsequent section.

To investigate the relative heterogeneity of immune cells between BRCA1-WT and BRCA1-MT samples, we explored the features and functions of T/NK cells. T/NK cells were re-clustered into 12 subpopulations, encompassing CD4+ naïve T cells (CD4+ Tn), CD4+ exhausted T cells (CD4+ Tex), CD4+ regulatory T cells (CD4+ Treg), CD8+ resident memory T cells (CD8+ Trm), four clusters of CD8+ effector memory T cell (CD8+ Tem), two clusters of CD8+ effector T cells (CD8+ Teff), CD8+ exhausted T cells (CD8+ Tex), NK cells, which displayed characteristic marker gene expression ( Figures 2A, B , Supplementary Figure S3A ). Notably, CD4⁺ Treg, CD4⁺ Tex, and CD8⁺ Tex cells exhibited high expression of canonical exhaustion markers, including PDCD1, HAVCR2, LAG3, CTLA4, and TIGIT ( Figure 2B ). Ro/e analysis further revealed that these subpopulations were preferentially enriched in the BRCA1-WT group, indicative of a more exhausted phenotype ( Figure 2C ). Furthermore, we found that the proportion of CD4+ Tex cells was positively correlated with the proportion of CD4+ Treg cells ( Supplementary Figure S3B ). The pseudotime analysis was conducted to deeply understand the immune dynamics. The results revealed that there were potential developmental branches, from CD4+ Tn to CD4+ Treg (Path1) or CD4+ Tex (Path2) ( Figure 2D ). We also identified the differential expression levels of marker genes during the differentiation and showed that FOXP3 and ICOS were strikingly upregulated along the branch of CD4+ Treg, consistent with their established roles in defining and sustaining the Treg lineage ( Supplementary Figure S3C ). It was suggested via trajectory analysis that the developmental trajectories were unclear in CD8+ T cells ( Supplementary Figure S3D ). Moreover, we employed a trajectory analysis based on CD8+ Tem and Teff cells, and the findings illustrated an apparent tendency to transform into CD8+ Tex subpopulation ( Figures 2E, F ). Moreover, Slingshot-based trajectory analysis was applied to further corroborate these results ( Supplementary Figures S3E, F ). Hence, we hypothesized that CD8+ Tex cells might derive from CD8+ Tem cells in TNBC, simultaneously developing the property of elevated exhaustion.

To assess the functional phenotype of T cells, we measured the T-cell-associated signature scores via ssGSEA. The results delineated that both CD4+ Tex and CD8+ Trm subpopulations had higher proliferation scores. The highest cytotoxicity scores were enriched in CD8+ Teff subpopulation, showing the key cell killing capability ( Figure 2F ). For CD8+ T cells, we found that the cytotoxicity score was negatively associated with naïveness score, in line with the poor cell killing capability of more naïve T cells ( Supplementary Figure S3G ). It was revealed via GO analysis of CD8+ T cells that pathways linked to immune negative regulation were preferentially congregated in Tex cells, metabolism-associated pathways were predominantly observed in Trm cells, and T cell activation pathways were enriched in other subpopulations ( Supplementary Figure S3H ). Signature scores were compared between the BRCA1-WT and BRCA1-MT groups. As depicted in Figure 2H , the naïve scores of T cells were substantially elevated in the BRCA1-MT group compared to the BRCA1-WT group, while higher exhaustion and proliferation scores of T cells were observed in the BRCA1-WT group. Intriguingly, we found that the cytotoxicity score of CD8+ T cells in the BRCA1-WT group was significantly higher than the BRCA1-MT group ( Figure 2H ). Consistent with prior research, dysfunctional T cells displayed higher proliferative potential ( Figures S3I, J ).

For Tem, Teff, and Tex subpopulations of CD8+ T cells, we dug into the dynamic landscape of exhaustion scores along developmental trajectory via pseudotime analysis, discovering a significantly more rapid ascent of exhaustion scores in BRCA1-WT group ( Figure 2I ). SCENIC analysis was conducted to further unravel the underlying mechanism of a more dysfunctional phenotype in the BRCA1-WT group. Specifically, genes within the interferon regulatory factor family (including IRF2 and IRF9) exhibited increased transcriptional activity in dysfunctional T cells which preferentially enriched in BRCA1-WT group ( Figure 2G ). Notably, this activation pattern was also confirmed at the gene expression level. We compared the expression profiles of IRF family genes in exhausted T cells between the two groups, and the results were consistent with our hypothesis ( Supplementary Figure S4 ), collectively contributing to the formation of the immune-suppress microenvironment (32, 33). These findings illustrated that BRCA1 potentially induced the immune dysfunction of T cells, leading to a better understanding for the immunotherapy of BRCA1-associated TNBC.

In addition, we delved into the underlying function of NK cell subpopulation between the two groups. It was revealed by GSVA that the BRCA1-WT group was characterized by more activated metabolism-related pathways, such as fatty acid metabolism, glycolysis, and oxidative phosphorylation ( Figure 2K ). This result was supported by the finding of elevated cytotoxicity scores of NK cells in the BRCA1-WT group ( Figure 2J ). Taken together, these evidences indicated the activation of NK cells in the BRCA1-WT group.

Accumulating evidences indicated that myeloid cells are the most abundant leukocytes in breast tumours (34), and we herein decoded the cellular complexity of myeloid cells in the patients with BRCA1-WT and BRCA1-MT. Myeloid cells were subclustered into monocyte, Tumor-associated macrophage (TAM) including SPP1+ TAM, C1Qs+TAM, HSPs+ TAM, ACP5+ TAM, and CCL5+ TAM, classical dendritic cell (cDC) encompassing cDC, CD1A+ DC, LAMP3+ DC and PCLAF+ DC, plasmacytoid DC (pDC) as well as unknown subpopulations ( Figure 3A ). Combined with the marker genes of each subpopulation ( Figure 3B , Supplementary Figure S5A ), the unknown cell subpopulation was regarded as double cells and excluded from the subsequent analysis.

In contrast to enhanced transcriptional perturbations of myeloid cells within the BRCA1-MT group ( Figure 1D ), the Ro/e analysis demonstrated that TAMs predominantly enriched in the BRCA1-WT group but not C1Qs+ TAM and CCL5+ TAM ( Figure 3C ). The trajectory analysis of monocytes and TAMs indicated that along seperate developmental paths, monocytes derived from three TAM subpopulations including SPP1+ TAM, C1Qs+ TAM, and ACP5+ TAM ( Figure 3D ). Subsequently, the signature scores in each TAM subpopulation were evaluated to illustrate the functional phenotype. As depicted in Figure 3E , the M1 and M2 scores showed no significance among these TAM subpopulations. However, we found a strikingly positive correlation between M1 and M2 scores, reflecting the cellular heterogeneity of TAMs within TME (r=0.31, p < 0.001) ( Figure 3F ). In response to the function of their respective marker genes, we observed that SPP1+ TAMs showed increased angiogenesis scores, whereas C1Qs+, HSPs+, and CCL5+ TAM subpopulations showed increased phagocytosis scores ( Figure 3E ). In contrast with the relationship between M1 and M2 scores, a significant negative correlation was observed in phagocytosis and angiogenesis scores, indicating their superiority in illustrating the features of TAMs (r = -0.29, p < 0.001) ( Figure 3G ). BRCA1 mutation might comprehensively enhance the activation of TAMs within the TME, as supported in Figure 3H . TAMs from BRCA1-MT tumors exhibited significantly higher M1 polarization scores, as well as elevated angiogenesis and phagocytosis scores, compared to those from BRCA1-WT tumors ( Figure 3H ). These results indicate that TAMs in the BRCA1-MT group are functionally more active, potentially contributing to the dynamic remodeling of the immune landscape.

We then interrogated the features of DC subpopulations and conducted a comparison between the two groups. Monocle3 discovered that PCLAF+ DC cells with high proliferation differentiated into LAMP3+ DC cells with comprehensively elevated function along pseudotime trajectory, which was validated by correlation analysis between signature scores and pseudotime ( Figures 3I-K , Supplementary Figures S5B - S4C ). Besides, there was a significant positive link between activation score and migration score as well as the tolerance score of DC cells ( Figure 3L , Supplementary Figure S5D ). PCLAF+ DC, as an immature DC subpopulation enriched in the BRCA1-WT group, induced the proliferation of Treg cells through interacting with Tumor cells ( Figure 3C ) (35). Awing to the immature property of PCLAF+ DC, we excluded this subpopulation and conducted the comparison of signature scores between the two groups. The result indicated that all three signature scores were strikingly higher in the BRCA1-WT group than in the BRCA1-MT group ( Figure 3M ). Taken together, we proposed the hypothesis that activated DCs in the BRCA1-WT group could mediate immune tolerance within the TME.

Next, we investigated the B cell subclusters in TNBC. The B cells were subdivided into six subpopulations based on canonical marker genes ( Figures 4A, B ). Utilizing Ro/e analysis, we found that plasma3 cells were preferentially distributed in the BRCA1-WT group ( Figure 4C ). FDCSP, one marker gene of this subcluster, was associated with the activation and differentiation of B cells ( Supplementary Figure S6A ). For delving into the complexity of function within B cells, GO analysis was employed. Plasma3 cells were enriched in ROS regulation pathways, such as response to oxidative stress, and detoxification, which served as the signal regulator of B cells (36) ( Figure 4D ).

Having evaluated immune cell heterogeneity, our attention subsequently shifted to epithelial cells. Tumor epithelial cells were distinguished from non-malignant epithelial cells via the CopyKAT package and were further divided into 19 cell subpopulations ( Supplementary Figures S6B, C ). Compared to the immune and stromal cells, significant differences were observed in epithelial cells across distinct specimens, revealing substantial intrasample heterogeneity ( Figures 4E , 1B ). The expression of signature for distinct tumor subsets was identified ( Supplementary Figure S6D ). Intriguingly, Epi15, characterized by high expression of ISG15 was found mainly in the BRCA1-WT4 specimen. ISG15 is regarded as a critical proto-oncoprotein that enhances the proliferation and metastasis of TNBC via inhibiting ubiquitin pathway (37, 38), and may also mediate tumor immunity via the JAK/STAT signaling pathway (39, 40).

We assessed the tumor heterogeneity scores across TNBC samples, demonstrating a significant positive correlation between intra- and inter-tumoral heterogeneity (r = 0.97, p < 0.001) ( Figure 4F ). Furthermore, marginally elevated heterogeneity scores were found in the BRCA1-WT group when compared to the BRCA1-MT group ( Figure 4G ). Intriguingly, the expression levels linked to human leukocyte antigen class I (HLA-I) displayed inconsistency across TNBC samples, implying potential deficits in the antigen-presenting capabilities of tumor cells ( Figure 4I ). Specifically, the mean expression level of HLA-A correlated positively with the percent of immune cells ( Figure 4H ). Despite higher immune cell proportions in BRCA1-MT samples, there were significant variations among these samples ( Figure 4I ).

Stromal cells in TME exhibited obvious functional and phenotypical heterogeneity and comprised multiple subsets. In stromal cells, we identified five cell types with distinct functional features, including myofibroblastic CAFs (myCAFs), inflammatory-like CAFs (iCAFs), pericytes, endothelial cells, and unknown cells ( Figure 5A ). The unknown subset was regarded as double cells and excluded from the subsequence analysis. As shown in Figure 5B and Supplementary Figure S7 , myCAFs were characterized by the high expression of POSTN and ACAT2, which showed immune-suppressive and pro-invasive TME (41), while iCAFs exhibited elevated cytokine-related gene expression, such as CXCL12, and CXCL14. It was revealed by distribution analysis that myCAFs were primarily gathered in the BRCA1-WT group, while iCAFs were mainly enriched in the BRCA1-MT group ( Figure 5C ). Notably, the iCAF3 subset, highly expressing MEG3, showed a similar distribution in the two groups. As previously characterized, MEG3+ CAFs could modulate metalloprotease-associated gene expression, including MMP1, MMP3, and MMP9, which was the signature of pre-CAFs (42, 43). Thus, iCAF3 subset was regarded as the original of CAFs, as evidenced by CytoTRACE ( Figure 5E ). Monocle analysis further demonstrated that the iCAF3 subset differentiated into myCAFs and iCAFs along separate pathways, as corroborated by Slingshot ( Figures 5D, F ). These findings suggested that BRCA1 mutation status might be a potential key driver of CAF differentiation direction. Using SCENIC analysis, iCAFs exhibited the activation of some essential TFs, such as JUN, FOSB, and FOS ( Figure 5G ), reflecting high response probability to ICIs (44). Furthermore, GSVA analysis revealed that the TNFα pathway was preferentially activated in iCAFs ( Figure 5H ). Taken together, we deduced that BRCA1 mutation mediated homologous recombination deficiency (HRD), thus promoting the activation of the TNFα pathway via modulating the JUN and FOS expression, and ultimately induces the iCAF phenotype.

For myCAFs that were enriched in the BRCA1-WT group, the TGF-β signaling pathway was found to promote the myCAF transformation. GSVA analysis indicated that epithelial-mesenchymal transition and apical junction pathways were also observed in myCAFs, suggesting that myCAFs might facilitate cancer metastasis ( Figure 5H ). We employed ssGSEA to assess the clinical significance of signature genes from myCAF in TNBC patients, demonstrating a correlation between elevated levels of this subset and poor prognosis ( Figure 5I ). In summary, our findings suggested that BRCA1 mutations could induce CAF reprogramming within the TME, thereby mediating distinct clinical outcomes in TNBC.

To dissect the intercellular interaction network in TNBC, CellChat was performed. Unsurprisingly, the BRCA1-MT group had more interaction strengths compared to the BRCA1-WT group, reflecting that BRCA1 mutation drove complex communication networks ( Figure 6A ). As shown in Figure 6B , the iCAFs exhibited the elevated strength of outgoing signals in the BRCA1-MT, while myCAFs had the increased strength of outgoing signals in the BRCA1-WT group. Furthermore, we identified the differential molecular crosstalk between CAF and other cell subpopulations ( Figures 6C, D ). Notably, the iCAF-mediated CXCL12-CXCR4 receptor/ligand pair was significantly enhanced in the BRCA1-MT group, while myCAF-mediated CD99-CD99 was strikingly increased in the BRCA1-WT group. More specifically, the receivers of increasing CXCL12 signaling network were observed to involve nearly all immune cell types in the BRCA1-MT group, reflecting its central role in driving immune regulation (45). For the CD99 signaling network, we observed decreased crosstalk communications between myCAFs and other immune cell types in the BRCA1-WT group, indicating its role in inhibiting immune infiltration ( Figure 6E ).

We reasoned that finer-resolution changes in cell subpopulations may be masked when analyzing the relative frequencies of broad cluster-based cell type identifications. To obtain subtler differences of cell subsets within TME, we grouped the cells into “neighborhoods” through k-nearest neighbor clustering, as implemented via the milo analysis framework. Differential abundance analysis of cell neighborhoods revealed that BRCA1-WT patients exhibited significant enrichment for myCAFs and BRCA1-MT patients were characterized by the enrichment of iCAFs ( Figures 7A, B ). To further elucidate the interactions among cell types at the spatial scale, four BRCA1-WT TNBC samples were included in this study. For every spot within the spatial data, we calculated the enrichment scores of each cell type, and correlation analyses were performed in each spatial RNA-seq sample ( Figure 7C ). The enrichment scores of myCAFs and iCAFs exhibited a significant positive correlation. In contrast, myCAFs and iCAFs demonstrated a clear spatial exclusivity ( Figure 7D ). These results reflected the homogeneity and heterogeneity within CAF subgroups. In addition, we observed a high degree of colocalization among CD4+ Tex, Bm, LAMP3+ DC, and CD8+ Tem, suggesting the presence of tertiary lymphoid structures ( Supplementary Figure S8A ). The mutual exclusivity relationship was also observed between epithelial cells and plasma cells ( Supplementary Figure S8B ).

We examined the enrichment scores of CAF subclusters in specific CAF spots. It was evident that myCAFs had substantially increased enrichment scores, consistent with their primary role in the CAF phenotype atlas of BRCA1-WT samples (p < 0.001) ( Figures 7B , E ). To pinpoint functional gene modules associated with the myCAF subcluster in the BRCA1-WT group, we conducted a co-expression network analysis using hdWGCNA on CAF cells, obtaining nine gene modules ( Figure 7F ). Notably, both CAF-M3 and CAF-M6 were significantly linked to the BRCA1-WT phenotype and the myCAF subcluster (log2FC > 1, p < 0.05) ( Figure 7G ). Employing a multi-scale identification framework, ISG15 was regarded as a crucial gene within myCAFs from BRCA1-WT patients, as devised in Figure 7H . Next, the association between ISG15 expression and myCAFs was also validated through spatial RNA-seq data. In spots expressing ISG15, the enrichment scores of myCAFs substantially exceeded others (p < 0.0001) ( Figure 7I ).

The expression level of ISG15 was significantly upregulated in the BRCA1-WT group compared to the BRCA1-MT group (p< 0.05) ( Figure 8A ). Given this, we reasoned that a subset of patients within the BRCA1-WT group, who exhibit BRCA1-like characteristics, may influence ISG15 expression. Thus, the TNBC samples were split into BRCA1-like and non-BRCA1-like groups as depicted by Chen et al (46). The comparison of ISG15 expression showed that the elevation of ISG15 expression was observed in the BRCA1-like group compared to the non-BRCA1-like group ( Supplementary Figure S9A ). We next evaluated the genetic vulnerabilities in eight TNBC cell lines via the CRISPR gene essentiality data extracted from the DepMap database. The result suggested that ISG15 was the most essential gene for SUM229PE cells, one of the BRCA1-WT cell lines ( Supplementary Figure S9B ). Our proteomics cohort (in-house cohort) indicated a substantial upregulation of ISG15 expression in breast cancer tissues compared to normal tissues, a finding supported by paired comparisons in the GSE109169 cohort. ( Figure 8B , Supplementary Figure S9C ). In summary, our study indicated that ISG15 was significantly upregulated in breast cancer, however, its functional role in TNBC specimens with BRCA1 mutation remained elusive.

Next, we evaluated the relationship between immunotherapeutic response and BRCA1-related classification systems including BRCA1 mutation and BRCA1-like via the TIDE (Tumor Immune Dysfunction and Exclusion) algorithm. No significant association between the scores for dysfunction and exclusion and the BRCA1 status (p > 0.05) ( Figure 8C ). Exclusion score was significantly elevated in the BRCA1-like group relative to the non-BRCA1-like group, and a higher dysfunction score was observed in the non-BRCA1-like group ( Supplementary Figure S9D ). These results indicated that the classification system based on BRCA1 status showed poor performance for immunotherapeutic response, and the development of a novel predictive model is particularly necessary.

As previous research indicated, ISG15 played a key role in regulating the immune infiltration of the TME, and its upregulation could synergically augment the therapeutic efficacy of ICIs (47). Utilizing protein-protein interaction analysis, we identified 50 ISG15-associated genes ( Supplementary Figure S9E ). Logistic regression was conducted 1,000 times to select the robust predictive features for immunotherapeutic response, yielding 11 genes as the input features for developing the model. Based on the 12 cohorts with immunotherapeutic response information, we used eight machine learning to devise the predictive model. A 5-fold cross-validation approach was employed to streamline the model and the ROC value was utilized to evaluate the performance of models ( Figure 8D ). Comparison of the performance of eight predictive models in the validation set showed that the greatest AUC of 0.60 (95% CI: 0.51–0.69) occurred in the KNN model, which was utilized for subsequent analysis ( Figures 8E, F ). For the testing set, the KNN model had an AUC of 0.61 (95% CI: 0.51–0.71), as shown in Figure 8G . To further explore the clinical significance of the KNN model, survival analysis was performed on TNBC patients from the METABRIC cohort. Increased response probability was associated with favorable overall survival (p < 0.05) ( Figure 8H ). In the TNBC patients from the TCGA cohort, we observed a similar trend of overall survival (p = 0.182) ( Figure 8I ).

Unlike the classification systems of BRCA1 mutation and BRCA1-like, the high response probability group exhibited significantly lower exclusion scores (p < 0.05) and slightly higher dysfunction scores (p = 0.19), as depicted in Figure 8J . Utilizing ssGSEA, we found that elevated immune infiltration scores were observed in the high response probability group relative to the low response probability group ( Figure 8K ).

The global feature importance for the KNN model was assessed using the DALEX package, identifying PCNA, STAT1, IFI6, and OAS1 as the four most significant genes ( Supplementary Figure S9F ). PCNA, a key immunological checkpoint molecule, is suggestive of a favorable response to NK cell-based immunotherapy (48). The stability of PCNA was augmented by ISG15 modification (ISGylation), substantiating the model’s validity in predicting immune therapy responses (49). ISG15 and STAT1 collectively activated IFN-related pathways, which promoted the release of IFNγ within immune cells, ultimately modulating accumulated immune infiltration (50). Furthermore, the expression of IFI6 was markedly positively associated with the abundance of B and T cells within TME, indicating a key role in immune regulation (51). The RNA sensor OAS1 and its upregulation responded to the amplification of IFN-1 and immune suppression genes, as evidenced by the partial elevation of immune suppressive cells in the high response probability group (50). Taken together, our predictive model has the potential to accurately identify TNBC patients who benefit from ICI therapy compared to the prediction based on BRCA1 mutation status.