Genes linked to programmed cell death influence immune system imbalance in triple-negative breast cancer

We downloaded the single-cell RNA sequencing dataset GSE176078↗ (14)from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo↗) and extracted nine single-cell samples of triple-negative breast cancer (TNBC). We calculated and filtered cells that expressed more than 300 genes using the PercentageFeatureSet function, with mitochondrial gene expression below 15%, and red blood cell gene proportion less than 1%. The merged ScRNA-seq data were normalized, and the top 2000 highly variable genes were identified using the FindVariableFeatures function. We then dimensionality reduction on the selected top 2000 highly variable genes using the ScaleData function and the RunPCA function. Batch correction was performed using the Harmony algorithm. Cells were clustered using the ‘FindNeighbors’ and ‘FindClusters’ functions (resolution=0.3) to identify cell clusters, and further dimensionality reduction was carried out using the UMAP method. Finally, the FindAllMarkers function was used to screen for marker genes in six subgroups, which were annotated and visualized using references and the CellMarker 2.0 database (15). After the inferCNV package was employed for tumor cell identification, we used the Monocle2 package (16) for pseudotime analysis of the subgroups.

Additionally, we obtained transcriptome data of 121 TNBC patients with survival information and normal patients from the University of California Santa Cruz (UCSC) Xena Browser (https://xenabrowser.net/↗). We also downloaded datasets GSE58812↗ and GSE21653↗ from the GEO database for subsequent transcriptome-level validation. All data required for this study are publicly accessible, and we provide all the codes in the attached materials (Supplementary Code Data Package). Differential gene expression between TNBC and normal groups was calculated using the limma package (17), with a selection criterion of p-value < 0.05.

The programmed death gene set encompasses 13 distinct modes of cell death (18), including apoptosis, pyroptosis, ferroptosis, autophagy, entotic cell death, cuproptosis, parthanatos, netotic cell death, alkaliptosis, lysosome-dependent cell death, and oxeiptosis, as well as disulfidptosis (Supplementary Table S1). We intersected this gene set with the differentially expressed genes (DEGs) identified in triple-negative breast cancer (TNBC). Subsequently, a univariate Cox regression analysis was performed on the intersected genes, resulting in the selection of 37 prognostic genes with a p-value < 0.05 for further analysis (Supplementary Table S2).

We employed non-negative matrix factorization (NMF) to decompose the gene expression matrix and extract biologically meaningful patterns (19). To explore the effects of selected programmed death genes on the tumor microenvironment (TME), we specifically extracted cells from tumor patients. We then applied dimensionality reduction using the NMF R package (version 0.20.6) to analyze the expression of 37 prognostic genes associated with programmed death within the TME. Following this analysis, distinct cell subtypes were identified based on the single-cell RNA (scRNA) expression matrix.

We utilized the FindAllMarkers function to identify marker genes associated with each non-negative matrix factorization (NMF) subtype within the TME. The NMF-derived cell subgroups were characterized based on the following criteria: (1) marker genes exhibited an absolute log2 fold change (Log2FC) greater than 1 and a p-value less than 0.05; (2) each subgroup was named after the gene with the highest Log2FC among the 37 programmed death genes; if no gene met the criteria within an NMF subgroup, it was designated as ‘none’. We then employed the Add ModuleScore function to compute the signature scores of differentially expressed genes (DEGs) for these NMF subgroups. The FeaturePlot function was applied to visualize the distribution of these DEGs. A list of marker genes for each cluster is provided in. 1

We utilized the SCENIC package, version R3.1.4, to dissect the regulatory network of transcription factors (TFs) within Triple-Negative Breast Cancer. For enrichment analysis, we leveraged gene-motif rankings from the RcisTarget database, focusing on the hg19-tss-centered-10kb track, to elucidate the regulatory network linked to transcription start sites (TSS) in TNBC. Visualization of the derived regulatory networks was facilitated by the pheatmap package.

CellChat enables the quantitative inference and analysis of intercellular communication networks from single-cell RNA sequencing (scRNA-seq) data (20). We utilized CellChat to assess the communication processes between all cells and tumor cells within the tumor microenvironment (TME), as well as the interactions between PCD-related cell clusters and malignant epithelial cell subpopulations. The strengths of cell-cell communication networks among all NMF clusters were visualized using the netVisualCircle function.

We utilized the cluster Profiler R package (21) to perform Gene Ontology (GO) annotation and Genes and Genomes (KEGG) pathway enrichment analysis on differentially expressed genes, considering a false discovery rate (FDR) threshold of p < 0.05 as statistically significant. Moreover, we employed the scMetabolism package (22) to assess the metabolic activity of macrophage NMF subtypes, encompassing 85 KEGG pathways and 82 REACTOME metabolic pathways. Additionally, we used the PROGENy R package (23)to calculate the activity of 13 cancer-related pathways (EGFR, MAPK, WNT, PI3K, VEGF, JAK-STAT, TGFβ, TNFα, NFκB, Hypoxia, Estrogen, p53, and Trail) across each NMF PCD-related subtype.

We classified the NMF PCD-related subtypes based on marker genes identified using the FindAllMarkers function in the TCGA database. Meanwhile, we conducted univariate Cox proportional hazards regression analysis to identify subtypes with prognostic significance. Subsequently, the survminer R package was employed to generate survival curves for the high- and low-risk groups stratified by these subtypes. Consistent with this approach, we also processed the GSE58812↗ and GSE21653↗ datasets.

We identified PCD-related clusters that exhibit high expression in tumors and possess prognostic significance, selecting one particular cluster for further analysis. We determined the prognostic genes by intersecting the markers of this cluster with the differentially expressed genes between TNBC and normal tissues, followed by univariate Cox regression analysis of genes with a p-value < 0.05, resulting in 14 identified genes (). 1

To develop a consensus Immune-Related Landscape Score (IRLS) with high precision and reliability, we integrated 10 machine learning algorithms to form 97 distinct algorithmic combinations. The algorithms included Random Survival Forest (RSF), Elastic Net (Enet), Lasso, Ridge, Stepwise Cox, CoxBoost, Partial Least Squares Regression for Cox (plsRcox), Supervised Principal Component (SuperPC), and Support Vector Machine for survival (survival-SVM). The signature generation process included the following steps: (a) Univariate Cox regression identified differential genes in the TCGA-TNBC cohort; (b) We then applied 101 algorithmic combinations to fit predictive models within leave-one-out cross-validation (LOOCV) framework specific to the TCGA-TNBC cohort; (c) Each model was independently validated using two external datasets (GSE58812↗, GSE21653↗); (d) For each model, Harrell’s Concordance Index (C-index) was calculated across all validation datasets, and the model demonstrating the highest average C-index was deemed optimal.

We utilized the ESTIMATE algorithm to assess the levels of stromal and immune cells in malignant tumor tissues based on expression data. The algorithm, obtained from the public website (https://sourceforge.net/projects/estimateproject/↗), estimates stromal and immune scores based on specific biomarkers associated with stromal and immune cell infiltration in tumor samples. These scores were analyzed separately to predict the levels of stromal and immune cells in tumor tissues. We calculated the StromalScore and ImmuneScore for each sample and then determined the correlation of key genes with these scores.

We employed the MCP-counter method (24), which robustly quantifies the absolute abundance of eight immune cell types and two stromal cell populations from transcriptome data (T cells, CD8T cells, Cytotoxic lymphocytes, B lineage, NK cells, Monocytic lineage, Myeloid dendritic cells, Neutrophils). Subsequently, we calculated the correlation of key genes with these populations.

Then, we used the ssGSEA method from the Gene Set Variation Analysis (GSVA) R package to identify genes associated with 28 immune cell types, as reported in the literature (25).

Additionally, we obtained the 50 HALLMARK pathways from the h.all.v7.5.symbols.gmt file at the GSEA website and used the ssGSEA method to calculate pathway scores for samples. We fianally determined the correlation of key genes with these pathways.

The MDA-MB-231 (SCSP-5043) and 4T1 (SCSP-5056) were purchased from the National Collection of Authenticated Cell Cultures. All cell lines were cultured in appropriate culture media supplemented with 1% Penicillin/Streptomycin and 10% fetal bovine serum (Gibco), and maintained at 37 °C and 5% CO2. The cell lines have been STR-authenticated and were routinely tested for mycoplasma-free.

Female BALB/c mice (6 weeks old) were purchased from Pengyue Experimental Animal Breeding Co, LTD (Jinan, China) and bred in animal facilities under specific pathogen-free conditions. All mouse procedures and experiments for this study were approved by the Animal Care and Utilization Committee of 960 Hospital. 4T1 cells in the exponential growth phase were digested with trypsin and washed three times with cold PBS. After centrifugation, cells were resuspended in cold PBS. 100 μl cold PBS containing 1×10⁶ 4T1 cells was orthotopically inoculated on the right flank of the fourth pair of mammary fat pads of mice. Tumor growth was monitored with a caliper and volume was calculated as (length × width²)/2. According to the experimental plan, tumor cells were sampled on the 7th, 14^th, and 21st days after injection, all tumor volumes did not exceed 1500 mm³.

Total RNA was extracted by using the PARKeasy Tissue/Cell RNA Rapid Extraction Kit (SPARKjade, AC0202, Shandong, China). The cDNA was synthesized using a cDNA synthesis kit (RR037A; Takara, Kyoto, Japan). The primers (Shanghai Bioengineering, Stable 1.), cDNA, TB Green^® Premix Ex Taq™ II (Takara, RR820A, Kyoto, Japan), and non-RNase dH2O were thoroughly mixed and then placed into a real-time PCR detection system, as instructed by the manufacturer. Real-time PCR mixture samples were detected by use of the SLAN-96P Real-Time PCR Detection System (SLAN-96P, Hongshi, Shanghai, China) the thermo-cycling conditions were initiated at 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Relative gene expression was calculated using the 2^−ΔΔct method.

Breast tumor tissues were frozen in liquid nitrogen and ground into a powder in a mortar. Powdered breast tumor tissues were lysed in RIPA buffer containing protease and phosphatase inhibitors (Selleck), and protein quantification was performed using the BCA Protein Assay kit (Solarbio). Next, proteins were isolated using 10% sodium dodecyl sulphate–polyacrylamide gel electrophoresis (SDS-PAGE gel, Solarbio) and transferred to a 0.45-μm polyvinylidene fluoride (PVDF) membrane (Merck Millipore, Billerica, MA, USA). After blocking in 5%(M/V) skim milk for 1 h, the samples were incubated with primary antibody at 4°Covernight. The membranes were then incubated with goat anti-rabbit/mouse IgG-HRP (Santa Cruz) at 37°C in 1:5000 dilution for 1h. Active bands were identified using an enhanced chemiluminescence (ECL) kit (Merck Millipore, Billerica, MA, USA).

Breast tumor tissues were collected and fixed in 4% paraformaldehyde, after which the tissues were then embedded in paraffin. 4μm thick formalin-fixed, paraffin-embedded tissue sections were used to perform immunohistochemical staining. After the deparaffinization and rehydration process, sections were heated in a pressure cooker in Tris-EDTA buffer for antigen retrieval and then were blocked by peroxidase endogenous blocking solution for 15 min and protein block for 1h. Sections were incubated with primary antibodies against CEBPB (proteintech, 23431-1-AP), and α-SMA (proteintech, 14395-1-AP) overnight at 4°C in a humid chamber. After washing with PBS, sections were incubated with HRP-conjugated Goat Anti-Rabbit/Mouse secondary antibody at room temperature for 1h, followed with DAB for 3–5 min under the microscope. Nuclei were counterstained with hematoxylin. Finally, sections were dehydrated in a series of alcohols and mounted in neutral gum. The specimens were observed using a confocal laser scanning microscope (CLSM; Olympus, Tokyo, Japan), and the area of staining was quantified using ImageJ software.

We performed single-cell RNA sequencing analysis of nine TNBC samples from the GEO TNBC dataset GSE176078↗. We identified T cells, B cells, myeloid cells, epithelial cells, fibroblasts, and endothelial cells as the main cell types in the TME (Figure 1A). Subsequently, the ClusterGVis package was used to visualize the marker genes for each cell type and the GO and KEGG pathway analyses were performed to identify the pathways enriched in each cell type. The results were consistent with the biological functions of the corresponding cell types (Figure 1B). Our data showed that all the epithelial cells in the TNBC samples exhibited copy number variations. Therefore, we concluded that the epithelial cells in the TNBC samples were all malignant. This finding was also consistent with their biological characteristics (Supplementary Figure S1). Next, to determine the differences in the immune microenvironment between the TNBC and normal tissues, we used the “ssGSEA” and “ImmuneCellA” packages to analyze and compare the immune cell populations in the TME between 121 TNBC patients with survival information and normal patients from the UCSC database (Supplementary Table S5). The results demonstrated that several immune cell types, including cytotoxic immune cells such as CD8+ T cells, CD4+ T cells, and NK cells, as well as B cells, macrophages, and fibroblasts were significantly increased in the TNBC tissues compared to the normal breast tissues, thereby indicating significantly altered immune microenvironment in TNBC (Figures 1C, D). In the subsequent analyses, we focused primarily on the CD8+T cells, CD4+T cells, NK cells, Tregs, B cells, macrophages, and fibroblasts. We also performed analysis to identify differentially expressed genes between TNBC and normal breast tissues. After intersecting these genes with a set of 13 PCD genes, we performed univariate Cox regression analysis to identify 37 PCD-related prognostic genes (Supplementary Table S2). We analyzed the expression levels of these prognostic genes at the single-cell level and found that the 37 PCD-related prognostic genes were highly expressed across various cell types (Figure 1E). This suggested that the 37 PCD-related prognostic genes played a key role in multiple biological functions of TNBC.

We extracted the stromal cells and identified subpopulations of macrophages through dimensionality reduction and annotation. Next, based on the PCD-related prognostic genes, non-NMF clustering was performed and four distinct cell clusters, namely, HMOX1-Mac-C1, NUPR1-Mac-C2, CEBPB-Mac-C3, and Non-Mac-C4, were identified. Pseudotime analysis demonstrated that different PCD-related prognostic genes were differentially expressed during various developmental phases of macrophages. Moreover, the HMOX1-Mac-C1 and Non-Mac-C4 clusters were located at the initial part of the developmental pseudotime trajectory and the NUPR1-Mac-C2 cluster was positioned at the terminal end of the trajectory (Figures 2A, B). Cell communication analysis suggested that NUPR1-Mac-C2 acted as a communication receptor and participated in the intercellular communication with the malignant epithelial cells through the MDK-(ITGA6+ITGB1), SPP1-CD44, and TNF-TNFRSF1A pathways (Figure 2C). SCENIC analysis identified differences in the transcriptional activity between various clusters. Specifically, the NUPR1-Mac-C2 cluster showed high expression of transcription factors such as STAT1, STAT2, and IRF7, whereas the CEBPB-Mac-C3 cluster showed high expression of transcription factors such as JUND, CEBPB, JUNB, and FOSB (Figure 2D). Based on the metabolic analysis, the NUPR1-Mac-C2 cluster exhibited high metabolic activity because of high expression of genes involved the TCA cycle and other pro-tumorigenic pathways (Figure 2E), thereby suggesting a close association with tumor growth. Progeny analysis demonstrated associations between the NUPR1-Mac-C2 cluster and pathways such as JAK-STAT, NF-κB, and TNF-α, thereby further substantiating its role in promoting oncogenesis (Figure 2F).

Based on secondary dimensionality reduction and clustering of the immune cells, we identified the following cell types: NK cells, CD8+ T cells, CD4+ T cells, regulatory T cells, naive T cells, and B cells (Figure 3A). Subsequently, we performed cell subtyping analysis to assess the functional profiles of each cluster and found that the VDAC3-CD8-C1, GZMB-CD8-C2, and IRF2-Treg-C1 clusters exhibited exhausted T cell phenotype, whereas the GZMB-NK-C3 and LAMTOR1-Treg-C2 clusters displayed cytotoxic phenotype. Furthermore, the GZMB-CD4-C1 cluster demonstrated both cytotoxic and exhausted T cell phenotypes, thereby suggesting functional differences among distinct PCD genes between various cell types (Figure 3B). To elucidate differences in the transcriptional activities between these cell clusters, we performed the SCENIC analysis and found that the clusters enriched for the cytotoxic cell or exhausted T cell functions exhibited robust transcriptional activity; among these, GZMB-CD4-C1 cluster showed the highest transcriptional activity (Figure 3C). Furthermore, there were differences in the transcription factor profiles between distinct T cell clusters with different functions (Figure 3C). Subsequent progeny analysis demonstrated that oncogenic pathways were significantly inhibited in subpopulations of cells with low transcriptional activity (Figure 3H). These findings were consistent with the SCENIC analysis. Among the B cells, the IRF2-B-C4 cluster exhibited highest level of transcriptional activity (Figure 3D). Cellular communication analysis demonstrated varying degrees of cellular communication between different T and B cell clusters and the epithelial cells (Figure 3E). Moreover, functional gene set analysis demonstrated functional differences in T cell activity among various subpopulations (Figure 3F). Pseudotime analysis demonstrated that the PCD-related prognostic genes played distinct roles at different time points across various cell types and influenced cluster differentiation (Supplementary Figure S2). Cellular communication analysis demonstrated that the immune cell subsets exhibited varied intensities of signaling output, primarily through the IL-6, CCL, and EGF pathways. When the immune cell subsets acted as signaling receptors, the pathways related to MIF and CXCL showed differing strengths, thereby indicating functional disparities (Figure 3G).

Pseudotime analysis results demonstrated that different PCD-related prognostic genes exerted their effects in the fibroblasts at distinct time points (Figure 4A). Subsequent NMF clustering categorized fibroblasts into the following six clusters: VDAC3-CAF-C1, NUPR1-CAF-C2, CEBPB-CAF-C3, LAMTOR1-CAF-C4, ITM2C-CAF-C5, and Non-CAF-C6. According to pseudotime analysis, VDAC3-CAF-C1 represented a stable subpopulation of fibroblasts throughout the developmental and differentiation process; Non-CAF-C6 cluster was positioned at the initiation of development; and CEBPB-CAF-C3 was located at the midpoint and represented a cluster that differentiated toward subsequent diverse subpopulations (Figure 4E). These results suggested a regulatory role for the PCD-related prognostic genes during fibroblast development and differentiation (Figure 4E).

To further elucidate the roles of different subpopulations, we assessed the biological functions of each subpopulation of fibroblasts (27). CEBPB-CAF-C3 exhibited biological functions associated with both pan-pCAF and pan-myCAF gene signatures that were characterized by genes related with smooth muscle generation and vascular wound healing, as well as, pan-dCAF and pan-iCAF-2 gene signatures that were related with genes related with the inflammation-related pathways and extracellular matrix (ECM) remodeling. The central location of CEBPB-CAF-C3 in the developmental branch also explains its diverse biological functions. And the CEBPB-CAF-C3 cluster was enriched in transcription factors such as JUND and ZC3H11A (Figure 4C). In contrast, VDAC3-CAF-C1 and ITM2C-CAF-C5 were significantly enriched in pathways associated with tumor-related inflammatory responses (Figure 4B). The SCENIC analysis demonstrated that the VDAC3-CAF-C1 cluster was enriched in transcription factors such as JUND, FOS, ATF3, and other transcription factors associated with the inflammatory pathways and was consistent with the classification. Further progeny analysis demonstrated that the CEBPB-CAF-C3 cluster was significantly associated with the PI3K-AKT pathway. This suggested that the CEBPB-CAF-C3 cluster may play a role in the occurrence and development of tumors through the PI3K-AKT pathway (Figure 4D). To further analyze the relationships between different clusters, we used the cell communication analysis and found that VDAC3-CAF-C1 was closely associated with malignant epithelial cells and various pro-oncogenic pathways (Figure 4F). Subsequently, based on the genes related to the tumor microenvironment (TME), we performed functional enrichment analysis of PCD-related fibroblasts. Genes associated with ECM, MMPs, TGF-β, and proinflammatory pathways showed significantly high expression in the CEBPB-CAF-C3 cluster, whereas proinflammatory-related genes were highly expressed in the VDAC3-CAF-C1 cluster (Figure 4G). These results aligned accurately with the classification of the CAF subtypes.

Transcriptomic analysis was performed to determine the prognostic significance of the PCD-related clusters. First, we validated the differential expression of PCD-related prognostic genes between the TNBC and normal breast tissues and found that the prognostic genes were significantly upregulated in the TNBC tissues (p < 0.001) (Figure 5A). We then performed Gene Set Variation Analysis (GSVA) using marker genes for each cluster to compare the transcriptomic differences and found that CEBPB-CAF-C3, IM2C-CAF-C5, VDAC3-CD8-C1, GZMB-CD8-C2, LAMTOR1-CD8-C4, ITM2C-CD8-C5, GZMB-CD4-C1, IRF2-Treg-C1, ITM2C-B-C1, IR2-B-C4, and NUPR1-Mac-C2 were all highly expressed in the TNBC tissues (Figure 5B). Prognostic analysis further demonstrated that the differential expression of CEBPB-CAF-C3 exhibited a consistent pattern. This suggested that CEBPB-CAF-C3 potentially represented a key oncogenic subgroup that impacted patient prognosis (Figures 5C–H). Further validation in the GSE58812↗ dataset confirmed that CEBPB-CAF-C3 was associated with poor prognosis and served as a risk factor subgroup (Figure 5I).

We identified a set of CEBPB-CAF-C3-specific marker genes and intersected them with the DEGs in the TNBC tissues compared with the normal breast tissues. Univariate Cox regression analysis was used to screen the PCD-related prognostic genes and 14 prognostic genes were identified (Supplementary Table S4). A prognostic model was constructed with these genes using an integrative approach based on machine learning. In the TCGA dataset, 97 predictive models were fitted using the Leave-One-Out Cross-Validation method and the concordance index (c-index) was calculated for each model across all the validation datasets (Figure 6A).

The model using a combination of Lasso and Stepwise Cox regression with backward elimination was the most optimal with the highest average concordance index (c-index) of 0.748. Among the 14 prognostic genes analyzed, we identified seven key genes, namely, CEBPB, TANK, MAT2B, TMEM165, CCDC167, CYFIP1, and COX17. Subsequently, the prognostic model was evaluated with the TCGA-TNBC dataset. The survival outcomes were poorer for the high-risk TNBC group compared to the low-risk TNBC group. Moreover, the prognostic model showed excellent prognostic prediction performance with area under the curve values of 0.92, 0.91, and 0.91 for 1-, 3-, and 5-year overall survival, respectively (Figures 6B, C). The prognostic model was further validated in the GSE58812↗ and GSE21653↗ datasets, which also showed poorer survival outcomes for patients in the high-risk group (Figures 6D, E). This demonstrated good predictive power for the prognostic model.

To further elucidate the differences in immune infiltration within the model, we used Cibersort, EPIC, MCPcounter, and Timer algorithms to analyze and assess the prognostic differences between high- and low-risk groups. Our data showed significant immune suppression and depletion of various immune cell types across all four algorithms in the high-risk group compared with the low-risk group (Figure 7A). We then assessed the differences in immune infiltration based on the seven hub genes and found that the risk genes (HR > 1) also demonstrated a degree of immune suppression (Supplementary Figure S3). Furthermore, we downloaded the HALLMARK pathway files (h.all.v7.5.symbols.gmt) comprising 50 pathways from the GSEA website and used ssGSEA to evaluate the scores of each pathway in the TCGA-TNBC dataset (26). We then calculated the Pearson correlation co-efficients to assess the relationship between the seven key genes and the pathway scores and visualized the results. Our results showed significant differences in the gene function enrichment and the risk genes were enriched in the STAT pathway and TNF-α-related pathways (Figure 7B).

We performed a subtype communication analysis to further investigate the differences in communication between immune cell types (T cells, B cells, and myeloid cells) and the CAF subtypes. The results demonstrated that CEBPB-CAF-C3 communicated with the T cells, B cells, and myeloid cells (Figure 8A). Furthermore, CEBPB-CAF-C3 as a signaling ligand communicated with the T cells through the MDK-(ITGA4+ITGB1) and MDK-NCL receptor-ligand pathways, and through MDK-SDC1 receptor-ligand interactions with the B cells, and through multiple pathways, including ANXA1-FPR1, MDK-(ITGA4+ITGB1), MDK-LRP1, MDK-NCL, and RARRES2-CMKLR1 with the myeloid cells (Figure 8B). However, when CEBPB-CAF-C3 was used as a signaling receptor, we did not observe any communication with the immune cells (Figure 8C).

To further characterize the immunotherapy outcomes between the high- and low-risk groups, we performed Tumor Immune Dysfunction and Exclusion (TIDE) analysis to compare the immunotherapy responses across different clusters. Patients in the high CEBPB-CAF-C3 group showed significantly poorer responses to immunotherapy compared to other groups (Figure 9A). This demonstrated functional disparities of various genes within distinct cell types. Subsequently, we used the IMvigor210 dataset to validate the efficacy of the prognostic model in predicting immunotherapy outcomes. The high-risk group in the IMvigor210 dataset was associated with a significantly lower response rate and poorer therapeutic outcomes than the low-risk group (Figures 9B, C). Furthermore, patients in the high-risk group from the IMvigor210 dataset showed significantly worse prognosis (p<0.001) based on the Kaplan-Meier survival analysis and the log-rank test (Figure 9D). This demonstrated the broad clinical applicability of the prognostic model. Finally, we performed a TIDE analysis of data from the GSE58812↗ dataset to validate the immunotherapy prognosis and showed that the CEBPB-CAF-C3 group was associated with a poorer prognosis than other groups (Figure 9E). This further substantiated the immunosuppressive role of CEBPB-CAF-C3 in TNBC.

To validate the above findings, we established three different stages of BC progression in a mouse model of TNBC, specifically on days 7, 14, and 21 post-tumor implantation, and assessed the expression levels of the fibroblast marker proteins, α-SMA and CEBPB. The results from immunohistochemistry (IHC) (Figure 10A), western blotting (Figure 10B), and multiple immunofluorescence (Figure 10C) experiments confirmed the presence of the CEBPB+CAF subtype in the TNBC mouse model. Subsequently, we performed quantitative reverse-transcription polymerase chain reaction (qRT-PCR) to examine the expression levels of the seven key prognostic genes—CEBPB, TANK, MAT2B, TMEM165, CCDC167, CYFIP1, and COX17—in the CEBPB+CAF prognostic model in the tumor tissues. The seven key prognostic genes were highly expressed in the murine BC model (Figure 10D). This suggested a positive correlation of the prognostic genes with the occurrence of TNBC. We further investigated the CEBPB+ CAF subtype in the TNBC patient tissues and included both cancerous and adjacent normal tissues. QRT-PCR (Figure 10E), Western blotting (Figure 10F), IHC (Figure 10G), and multiple immunofluorescence (Figure 10H) experiments confirmed the high expression of the CEBPB+ CAF subtype in the BC tissues.

These results demonstrated the presence of the CEBPB+ CAF subtype in TNBC and its association with tumor progression. Therefore, the CEBPB+CAF subtype may serve as a valuable prognostic predictor for TNBC.