Single-cell and bulk RNA sequencing analysis of B cell marker genes in TNBC TME landscape and immunotherapy

Objective: This study amied to investigate the prognostic characteristics of triple negative breast cancer (TNBC) patients by analyzing B cell marker genes based on single-cell and bulk RNA sequencing.

Methods: Utilizing single-cell sequencing data from TNBC patients, we examined tumor-associated B cell marker genes. Transcriptomic data from The Cancer Genome Atlas (TCGA) database were used as the foundation for predictive modeling. Independent validation set was conducted using the GSE58812 dataset. Immune cell infiltration into the tumor was assessed through various, including XCELL, TIMER, QUANTISEQ, CIBERSORT, CIBERSORT-ABS, and ssGSEA. The TIDE score was utilized to predict immunotherapy outcomes. Additional investigations were conducted on the immune checkpoint blockade gene, tumor mutational load, and the GSEA enrichment analysis.

Results: Our analysis encompassed 22,106 cells and 20,556 genes in cancerous tissue samples from four TNBC patients, resulting in the identification of 116 B cell marker genes. A B cell marker gene score (BCMG score) involving nine B cell marker genes (ZBP1, SEL1L3, CCND2, TNFRSF13C, HSPA6, PLPP5, CXCR4, GZMB, and CCDC50) was developed using TCGA transcriptomic data, revealing statistically significant differences in survival analysis (P<0.05). Functional analysis demonstrated that marker genes were predominantly associated with immune-related pathways. Notably, substantial differences between the higher and lower- BCMG score groups were observed in terms of immune cell infiltration, immune cell activity, tumor mutational burden, TIDE score, and the expression of immune checkpoint blockade genes.

Conclusion: This study has established a robust model based on B-cell marker genes in TNBC, which holds significant potential for predicting prognosis and response to immunotherapy in TNBC patients.

Introduction

Triple-negative breast cancer (TNBC) is a distinct subtype of breast cancer characterized by the absence of estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2) expression in immunohistochemistry (1). Accounting for approximately 15–20% of all breast cancers, TNBC exhibits aggressive clinical symptoms (2). Although chemotherapy remains a major therapeutic method for metastatic TNBC; its efficacy is limited, yielding a median overall survival time of 12-18 months (2). Consequently, there is an urgent need for innovative therapeutic strategies.

In recent years, immunotherapy has emerged as a promising intervention, demonstrating prolonged survival in various solid tumors (3–5). Notably, immune checkpoint inhibitor (ICI) monotherapy in TNBC patients has shown response rates ranging from 5% to 23% (6), while combined with chemotherapy, for early-stage TNBC patients exhibit pathologically complete response rates between 22% and 60% (7). Despite these advances, only a subset of patients benefits from immunotherapy, highlighting the need for the development of predictive models and the identification of novel biomarkers to better anticipate treatment outcomes and prognosis.

Single-cell RNA sequencing (scRNA-seq) technology and related data analysis methods has facilitated the exploration of molecular characteristics of immune cells within the tumor microenvironment (TME) (8). This unparalleled opportunity has provided insights into cancer immunity, enabling the establishment of genetic markers based on immune cell molecular characteristics to predict immunotherapy outcomes in cancer patients (9, 10).

Recent research has demonstrated that tumor-infiltrating B lymphocytes (TIL-Bs) in breast cancer, responding to B cell receptor (BCR) activation and generating immunoglobulin (Ig) in vivo (11–13). Thus, B cells may significantly influence the prognosis of breast cancer, particularly in patients with immunogenic TNBC. Here, our study aimed to conduct an integrated analysis of scRNA-seq data from TNBC samples, identify B cell marker genes, and subsequently develop a predictive model for accessing the prognosis of TNBC patients.

Methods

Data collection

Transcriptional RNA sequencing data from single-cell profiling datasets of four TNBC tumor tissue (GSM4909281, GSM4909282, GSM4909283 and GSM4909284; accessible at https://www.ncbi.nlm.nih.gov/geo/) were acquired for the purpose of investigating B cell marker genes. Concurrently, clinical profiles and transcriptomic data were obtained from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases.

Processing single-cell sequencing data

The single-cell sequencing data obtained from patients were subjected to comprehensive analysis through the R software and related packages. Seurat objects were constructed following the parsing of single-cell sequencing data sourced from TNBC samples represented by GSM4909281, GSM4909282, GSM4909283, and GSM4909284. The DoubletFinder package was applied to eliminate the cell doublets within the sample. Exclusion criteria for cells of suboptimal quality are as follows: 1) the number of features exceeding 500 and falling below 6000, 2) the mitochondrial genes expression below 10%, 3) the expression of erythroid genes less than 5%. In order to integrate data from multiple samples, the software known as “harmony” was employed.

The Uniform Manifold Approximation and Projection (UMAP) method was adopted to reduce the number of dimensions, displaying clustered cells with only two dimensions on a map. The subpopulations of tumor-associated B cells were determined by using the “SingleR” software as previous publications indicated (14, 15). Marker genes for different cell types were ascertained through Wilcoxon-Mann-Whitney test and the “FindAllMarkers” function analysis. The criteria for filtering marker genes of various cell types included a | log2 (fold change) | > 1 and the adjusted P value less than 0.05.

Functional enrichment

The functional annotation of marker genes associated with tumor-associated B cells was augmented through enrichment analysis using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The c2.cp.v2023.2.Hs.symbols.gmt dataset from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) was used for Gene Set Enrichment Analysis (GSEA) to elucidate the biological processes involved in the B-cell marker gene.

A predictive model of B cell marker gene

A predictive model for B cell marker gene expression in TNBC was developed to predict TNBC prognosis. In the univariate COX regression analysis, genes identified as B-cell markers were incorporated. The least absolute shrinkage and selection operator (LASSO) regression was performed to construct a penalty function, effectively compressing the coefficients of variables and reducing overfitting of the model caused by prognosis-related genes. Subsequently, the outcomes from the multivariate COX regression analysis were used to formulate a predictive B cell marker gene score (BCMG score) model.

The formula is as follows:

$$BCMG score={\beta }_1* Exp{r}_1+{\beta }_2* Exp{r}_2+\dots +{\beta }_n*Exp{r}_n$$

Here, Expr denotes the mRNA expression of the crucial gene, n denotes the number of genes included in the model, and β represents the associated regression coefficient determined in the multivariate gene COX regression analysis. The data was then classified into higher and lower- BCMG score groups based on the median value of the BCMG score.

To validate the model’s robustness, an independent third-party validation set, GSE58812, was employed to assess and confirm the predictive performance of the model.

Immune infiltration analysis

Based on the gene expression dataset of patients, the CIBERSORT algorithm was conducted to determine the magnitude of immune cell proliferation in each sample. Spearman correlation analysis was performed to examine the variations in immune cell infiltration between the two groups. A single-sample GSEA was carried out to further assess the differences in immune cell activity between the higher and lower- BCMG score groups. Estimation techniques were implemented to determine the differences of immune scores, stromal scores, overall scores, and tumor purity between the two groups. Multiple tools, including XCELL, TIMER, MCP counter, CIBERSORT, and CIBERSORT-abs, were applied to establish a correlation between B cell marker genes and immune cells.

Forecasting the patient’s response to immunotherapy

Spearman test was used to analyze the correlation of expression between model and immune genes. The TIDE score was incorporated to evaluate the effectiveness of immunotherapy in individual patients. Patients with lower TIDE scores indicated a decreased risk of immunological escape, which suggested a higher likelihood of effective immunotherapy. An online tool TIDE (http://tide.dfci.harvard.edu/) was used to examine the immunotherapy scores of each patient, and investigate the differences in immunotherapy results between the higher and lower- BCMG score groups.

Analysis of the tumor’s mutational burden

The raw data of the patient’s tumor mutations were obtained and downloaded from the TCGA database. The “maftools” package in R was used to generate waterfall graphs in accordance with predefined specifications. This study also further assessed the differences in tumor mutation burden as well as predicted discrepancies between higher and lower- BCMG score groups.

Statistics analyses

Kaplan-Meier method and the log-rank test were utilized for survival analysis. Alternately, the Wilcoxon rank-sum test was employed to compare the differences existing between the two groups. R (4.1.2, available at https://www.r-project.org/) was conducted throughout the entire data analysis process. P values of less than 0.05 on both sides were defined as statistically significant.

Results

Examination of single cells

The integration of single-cell sequencing data from four TNBC patients was achieved using the “Harmony” package. After excluding cells of suboptimal quality, the subsequent study consisted of the examination of 22,106 cells and 20,556 genes (Figure 1). The identification of potential marker genes for TNBC-associated B cells was conducted through the application of the Wilcoxon-Mann-Whitney test, revealing 116 distinct genes with significant differences (Table S1). The clinical characteristics of the training group TCGA and the testing group GSE58812 was detailed in Table S2.

FIGURE 1

Figure 1 Single cell sequencing profile of 22,106 cells from human samples of TNBC. (A–C) UMAP plots of 22,106 tumor-associated B cells here, colored by individual samples (A), 22 cell clusters identified when resolution equals 2 (B) and cell types (C). (D) Heatmap displaying the top 10 marker genes in each cell cluster. (E) Proportions of each cell type in each sample colored by cell types. (F) Expression analysis of BCMGs in 22 cell clusters. The intensity of the color indicates the average expression of the gene. Dot size indicates the percentage of cells expressing the gene.

Tumor-associated B-cell marker gene enrichment analysis

The investigation of genes enriched in tumors associated B cells yielded notable results. The GO and KEGG analysis revealed that marker genes of B cells predominantly influenced functional pathways specific to B cells. Significant differences were observed in cell pathway activity scoring across different cell types. Particularly, the unfolded protein response signaling pathway demonstrated notable activity in B cells. The B cell marker genes mainly participate in biological processes such as pid CD8 tcr dowstream pathway, pid IL12 pathway, WP cancer imunotherapy by PD 1 blockade, and WP T cell recpto signaling pathway. (Figure 2).

FIGURE 2

Figure 2 Exploring the differences in activity in different cell types by employing different enrichment analysis methods. Histogram of GO (A), KEGG (B) and GSEA (C) enrichment analysis.

Development of a predictive model for B-Cell marker gene expression

The evaluation of the 116 genes involved a univariate COX regression to identify the 10 genes associated with prognosis. In order to prevent overfitting, a lasso regression analysis was subsequently performed. The results of the lasso regression were further subjected to a multifactor COX regression analysis to construct predictive models. The final model incorporated the following variables: ZBP1, SEL1L3, CCND2, TNFRSF13C, HSPA6, PLPP5, CXCR4, GZMB, and CCDC50 (Figures 3A, B).

FIGURE 3

Figure 3 Construction of BCMG score model in the TNBC cohort. (A) Log (λ) change curves of regression coefficients and cross-validation for optimizing the parameter in LASSO regression. (B) Tenfold cross-validation of adjusted parameter choices in lasso regression. (C, D) Kaplan-Meier survival analysis of TCGA and GSE58812 of TNBC patients stratified by higher and lower- BCMG score. (E, F) ROC curve analysis for predicting the risk of death in the cohort of TCGA (E) and GSE58812 (F).

The predictive score was computed based on the expression levels of the selected genes using the following formula:

$$\begin{array}{c}\text{BCMG score}=(-0.0390\times ZBP1)+(-0.3450\times SEL1L3)+(-0.0474\times CCND2)+(-\\ 0.0716\times TNFRSF13C)+(0.2560\times HSPA6)+(0.5725\times PLPP5)+(-0.1924\times \\ CXCR4)+(-0.1424\times GZMB)+(-0.1246\times CCDC50)\end{array}$$

Patients from TCGA and GEO were divided into higher and lower- BCMG score groups according to the median value of the scores obtained from TCGA.

Model validation

The survival curves revealed that patients classified in the higher BCMG score group had a worse prognosis than those in the lower BCMG score group (Figures 3C, D). The time-dependent receiver operating characteristic (ROC) curves indicated that the predictive model had an outstanding discriminatory capability (Figures 3E, F). Importantly, these results were not only internally validated within the study group but also externally validated in an independent validation group (Figures 3E, F). According to the findings of the independent predictive study (Figures 4A, B), the predictive model was identified as a significant independent predictive factor in patients with TNBC. The correlation analysis with clinical characteristics further highlighted the close association between the predictive model and lymph node metastases (Figure 4). In-depth exploration of the correlation between the model genes and patient prognosis (Figure 5) revealed that the expression levels of ZBP1, SEL1L3, CCND2, TNFRSF13C, CXCR4, GZMB, and CCDC50 were higher in the lower BCMG score group compared to the higher group. Conversely, the higher BCMG score group exhibited significantly higher levels of expression for HSPA6 and PLPP5 (Figure 6).

FIGURE 4

Figure 4 Clinical analysis in the cohort of TCGA. (A) Univariate Cox regression analysis revealed the association between patients’ survival and clinicopathological parameters along with BCMG score. (B) Multivariate Cox regression analysis uncovered that only the BCMG score (P< 0.001) was an independent prognostic factor for TNBC patients. (C) The prediction of 3-, 5-, and 7-year survival for TNBC patients based on the prognostic nomogram derived from the BCMG score and other clinicopathologic feature. (D) Clinical correlation with age, stage, T, N, M and BCMG score shown in heatmap. (E) DCA curve analysis. (F) Calibration Curve illustrated the consistency between predicted and observed 3-, 5-, and 7- year survival rates in TNBC patients depending on the prognostic nomogram.

FIGURE 5

Figure 5 Survival analysis in TNBC patients based on BCMG score. Kaplan-Meier survival analysis of TCGA TNBC patients stratified by high and low CCDC50, CCND2, CXCR4, GZMB, HSPA6, PLPP5, SEL1L3, TNFRSF 13C, and ZBP1 (A–I), which were based on the best cutoff values.

FIGURE 6

Figure 6 Prognostic analysis of TNBC patients on the basis of BCMG score in TCGA and GEO. BCMG score distribution in TCGA (A) and GSE58812 (B). Scatterplot of survival status and survival time of TNBC patients in TCGA (C) and GSE58812 (D). Heatmap generated on the basis of identified gene expression in TCGA (E) and GSE58812 (F).

Analysis of tumor mutation burden

The analysis of the tumor mutation burden revealed an inverse relationship with the prognostic model. The waterfall plot visually illustrated that the mutation types were not identical between the higher and lower- BCMG score groups. The survival analysis results indicated that patients with lower tumor mutation load exhibited a more favorable prognosis (Figure 7).

FIGURE 7

Figure 7 Tumor mutation burden analysis. (A, B) Differential counting of tumor mutation burden between higher and lower - BCMG score groups. (C) Tumor mutation load comparison. (D) Correlation analysis of BCMG score and mutation burden. (E, F) Prognostic analysis of tumor mutation load in higher and lower- BCMG score groups.

Analysis of TME

The CIBERSORT analysis demonstrated the differential presence and activity of various immune cell types in the TMB between the higher and lower- BCMG score groups. Specifically, CD8⁺ T cells, activated memory CD4⁺ T cells, M1 macrophages, M2 macrophages and eosinophils were found to serve different functions between the two groups (Figures 8A, B). The overall activity of the majority of immune cells exhibited notable changes in both groups, further emphasizing the dynamic nature of the TME (Figure 8C). The ESTIMATE algorithm indicated that the lower BCMG score group had significantly higher tumor scores, interstitial scores, total scores, and tumor purity scores than the higher BCMG score group (Figure 8D). Furthermore, correlation analysis between the predictive scores and the TME analyzed by the four methods CIBERSORT-ABS, CIBERSORT, QUANTISEQ, and XCELL, demonstrated the potential influence of B-cell marker gene model on the TME (Figure 9).

FIGURE 8

Figure 8 Correlation of the B cell prognostic model genes with immune cell infiltration and immune checkpoint in TCGA cohort. (A) Analysis of immune cell infiltration in higher and lower- BCMG score groups. (B) The comparison of 22 immune cells’ infiltration level in higher and lower- BCMG score groups. (C) The comparison of 28 immune cells’ infiltration level in higher and lower- BCMG score groups. (D) TME score differences compared based on the ESTIMATE algorithm. *P < 0.05; **P < 0.01; ***P < 0.001.

FIGURE 9

Figure 9 BCMG score and immune cell correlation analysis. (A) Spearman correlation analysis showed that BCMG scores strongly correlated with tumor-infiltrating immune cells. (B–G) Correlation of predictive models and immune cells based on XCELL, TIMER, QUANTISEQ, CIBERSORT, and CIBERSORT-abs.

Predictions regarding immunotherapy

Correlation analysis revealed that ZBP1, GZMB, CCND2, and TNFRSF13C exhibited strong correlations with the majority of the genes that inhibit immune response (Figure 10A). Simultaneously, a significant difference was observed between the higher and lower- BCMG score groups concerning the expression of several immune blockade checkpoint genes (Figure 10C). The TIDE ratings were used to make predictions about the outcomes of immunotherapy treatment for patients classified as either high or low BCMG expression. The TIDE scores of patients in the high expressed group were significantly higher than those of patients in the low expressed group, which indicated a higher likelihood of immune evasion in the high expressed group and immunotherapy may be less effective for patients at high expression. In addition, there were observed differences between the two groups in terms of CD8 scores, CD274 ratings, MDSC scores, T cell dysfunction scores, and T cell exclusion scores (Figure 10B).

FIGURE 10

Figure 10 Prediction of response of TNBC patients to immunotherapy. (A) Analysis of differences in TIDE, IFNG, MSI Expr, Merck18, CD274, CD8, Dysfunction, Exclusion, MDSC, CAF, and TAM M2 between higher and lower- BCMG score groups. (B) Analysis of immune checkpoint blocking gene expression between higher and lower- BCMG score groups. (C) Relevance of model genes and immune checkpoint blocking gene.

Discussion

In this study, we investigated the role of B-cell marker genes in TNBC using scRNA-seq technology. We successfully established a novel prediction score for TNBC patients, derived from the identified B-cell marker genes within the TCGA database. This prediction score was rigorously validated for its predictive capabilities in an independent cohort of the GEO dataset. Our analysis highlighted specific genes, including ZBP1, SEL1L3, CCND2, TNFRSF13C, CXCR4, GZMB, and CCDC50, which exhibited protective characteristics. In contrast, HSPA6 and PLPP5 were identified as having negative implications. Subsequently, the nine gene signatures collectively contribute to a robust predictive model that emerges as a promising predictor for clinical prognosis in TNBC patients.

Patients with TNBC exhibit a significantly higher level of immune infiltration compared to other breast cancer subtypes, despite their overall poor prognosis. The higher immune infiltration is characterized by elevated tumor-infiltrating lymphocytes (TILs) (16), increased levels of PD-L1 expression in tumor cells and immune cells (17, 18), and a greater number of non-synonymous mutations, which produce tumor-specific neoantigens, activate neoantigen-specific T cells, and trigger anti-tumor immune responses (19, 20). This unique immune landscape suggests a greater potential for benefit from immunotherapy in TNBC comparison to other subtypes. ICIs can block immunosuppressive receptors and increase the cellular toxicity and proliferation of TILs (3, 21, 22). According to the findings of a few studies, the effectiveness of ICI monotherapy ranges from 5% to 23% (6), whereas the effectiveness of combination therapy ranges from 22% to 60% (7). Nonetheless, some patients experienced resistance during or after treatment, leading to immune escape and tumor recurrence. Dual checkpoint blockade, which have synergistic anti-tumor effects in advanced malignancies, such as anti-CTLA-4 and anti-PD-1, could potentially serve as a more effective therapeutic strategy (23). Furthermore, recent studies have explored innovative strategies, such as bispecific antibodies targeting TGF-β and PD-L1 (BiTP), which have demonstrated potent antitumor activity in TNBC. BiTP in murine TNBC models showed higher anti-tumor activity compared with solo anti-PD-L1 or anti-TGF-β (24). The rejection and exhaustion of CD8⁺ T cells are two key factors that contribute to the reduced tolerance of ICI therapy. When enhanced CD8⁺ T-cell infiltration occurs, there is a shift towards heat immunity (25). BiTP has a high binding affinity to the dual targets, which reduces collagen deposition, enhances CD8⁺ T-cell penetration and increases TILs (26).

Our study provides valuable insights into the molecular landscape of TNBC-associated B cells and their functional implications in tumor biology. These findings underscore the intricate interplay between the identified prognostic model and the immune landscape within the TME, and we also illustrate potential immunomodulatory roles of the B-cell marker genes in the context of TNBC. Moreover, the differential mutation burden observed in distinct BCMG score groups highlights the importance of considering mutational characteristics in conjunction with the prognostic model for a comprehensive understanding of the disease prognosis in TNBC patients.

Extensive investigations in the breast cancer TME have primarily focused on CD8⁺ T lymphocytes and natural killer (NK) cells (27). In contrast, TIL-B lymphocytes (TIL-Bs) cells may exhibit antigen-induced phenotype (28), and it is hypothesized that autoantibodies play a role in initiating tumor cell clearance (29). TIL-Bs have also been identified as antigen-presenting cells, contributing to the stimulation of an immune response against tumors (30). Moreover, there is a potential for TIL-Bs to form stromal clusters with T cells, engaging in functional cross-talk in both directions. The formation of BCR-immune complexes may also be facilitated by TIL-Bs through the upregulation of BCR pathway components (1). Besides, CD8⁺ T cells represent pivotal defense cells in TME and establish a functional cycle with dendritic cells (DCs) and NKs (31). Positive crosstalk is observed between CD8⁺ T cells and DCs, with CD8⁺ T cells inducing NK cell activity (32). Therefore, it is possible that B cells exert a significant influence on the prognosis of breast cancer, particularly in patients of immunogenic TNBC.

The scRNA-seq technology has allowed for a more in-depth exploration of the molecular properties of immune cells within the TME, offering potential insights into novel biomarkers (33). In this study, the investigation of B-cell marker genes provides valuable insights into the potential role of nine specific genes in influencing the prognosis of individuals diagnosed with TNBC. ZBP1, overexpressed in necrotic breast tumors, has been linked to the induction of interferon. In the absence of ZBP1, both the process of tumor necrosis and the suppression of metastasis are prevented (34). SEL1L3, with decreased expression in cancer cells during sustained endoplasmic reticulum stress (35), may serve as a predictive tool for survival outcomes and immunotherapy response in various cancers (36–39). The promoter of CCND2 is hypermethylated (40–42) in solid tumors, which leads to CCND2 hypo-expression (40, 43, 44) and promotes cell proliferation (45). TNFRSF13C could encode B-cell activating factor receptor, which regulates B-cell proliferation, development and maturation (46, 47). High levels of CXCR4 levels were expressed in more than 40% of breast tumor tissues (48) and in 75% of TNBC patients (49), and overexpression of CXCR4 in cancer cells contributes to tumor growth, invasion, metastasis, and recurrence (50). GZMB, known for stimulating anti-tumor immune responses and inhibit tumor growth (51), is correlated with favorable prognosis in breast cancer tissues (52, 53). CCDC50 is known to mediate apoptosis via the NF-κB pathway (54) and had prognostic predictive value in lung adenocarcinoma (55). The two-way role of HSPA6 in tumor, acting as both a potential target for tumor inhibition and risk factor of tumor development and tumor progression, underscores the complexity of its functions in cancer biology (56). Our study suggested a negative correlation between high HSPA6 expression and prognosis, which contradicts the findings reported by Shen et al., where high HSPA6 expression was positively correlated with longer overall survival (57). This discrepancy may arise from the multifaceted functions of HSPA6, affecting different aspects of cancer biology. PLPP5 participates in the regulation of many cancer-associated transduction pathways such as the JAK/STAT (58), and its overexpression in pancreatic and small cell lung cancer cells may promote proliferation and survival (59). However, in breast cancer cell lines, down-regulation of PLPP5 inhibits tumor growth and increases apoptosis (60), which is consistent with our study.

Despite these intriguing findings, there are certain limitations in the current study. Additional research is needed to get a deeper understanding of molecular interactions and mechanisms underlying these gene functions. Furthermore, it is essential to verify these findings in cell lines, animal models, and tissue samples to ensure their relevance and applicability to clinical settings.

Conclusion

In conclusion, this study established a robust model utilizing B cell marker genes to predict prognosis and immunotherapy response in TNBC patients, contributing novel insights into the understanding of TNBC molecular landscape and its implications for therapeutic interventions. The discovery of these gene signatures adds to our understanding of TNBC and offers a promising avenue for improving patient outcomes and guiding therapeutic strategies in this challenging subtype of breast cancer. Further research and clinical validation may continue to refine the application of these gene signatures in the clinical management of TNBC patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

Conceptualization: CZ and HL; Data curation: TX; Methodology: XW; Formal analysis: YL; Writing-original draft: FZ; Review and editing of the manuscript: CZ and XL. All authors contributed to the article and approved the submitted version.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The study was financed by Natural Science Foundation of Hubei Province (2023AFB165) and Renmin Hospital of Wuhan University Cross-Innovation Talent Project (JCRCWL-2022-003).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1245514/full#supplementary-material

References

1. Harris RJ, Cheung A, Ng JCF, Laddach R, Chenoweth AM, Crescioli S, et al. Tumor-infiltrating B lymphocyte profiling identifies IgG-biased, clonally expanded prognostic phenotypes in triple-negative breast cancer. Cancer Res (2021) 81(16):4290–304. doi: 10.1158/0008-5472.CAN-20-3773

CrossRef Full Text | Google Scholar

2. Garrido-Castro AC, Lin NU, Polyak K. Insights into molecular classifications of triple-negative breast cancer: improving patient selection for treatment. Cancer Discovery (2019) 9(2):176–98. doi: 10.1158/2159-8290.CD-18-1177

CrossRef Full Text | Google Scholar

3. Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med (2018) 24(9):1449–58. doi: 10.1038/s41591-018-0101-z

CrossRef Full Text | Google Scholar

4. Motzer RJ, Rini BI, McDermott DF, Redman BG, Kuzel TM, Harrison MR, et al. Nivolumab for metastatic renal cell carcinoma: results of a randomized phase II trial. J Clin Oncol (2015) 33(13):1430–7. doi: 10.1200/JCO.2014.59.0703

CrossRef Full Text | Google Scholar

5. El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet (2017) 389(10088):2492–502. doi: 10.1016/S0140-6736(17)31046-2

CrossRef Full Text | Google Scholar

6. Keenan TE, Tolaney SM. Role of immunotherapy in triple-negative breast cancer. J Natl Compr Canc Netw (2020) 18(4):479–89. doi: 10.6004/jnccn.2020.7554

CrossRef Full Text | Google Scholar

7. Nanda R, Liu MC, Yau C, Shatsky R, Pusztai L, Wallace A, et al. Effect of pembrolizumab plus neoadjuvant chemotherapy on pathologic complete response in women with early-stage breast cancer: an analysis of the ongoing phase 2 adaptively randomized I-SPY2 trial. JAMA Oncol (2020) 6(5):676–84. doi: 10.1001/jamaoncol.2019.6650

CrossRef Full Text | Google Scholar

8. Chen H, Ye F, Guo G. Revolutionizing immunology with single-cell RNA sequencing. Cell Mol Immunol (2019) 16(3):242–9. doi: 10.1038/s41423-019-0214-4

CrossRef Full Text | Google Scholar

9. Zhong R, Chen D, Cao S, Li J, Han B, Zhong H. Immune cell infiltration features and related marker genes in lung cancer based on single-cell RNA-seq. Clin Trans Oncol (2021) 23(2):405–17. doi: 10.1007/s12094-020-02435-2

CrossRef Full Text | Google Scholar

10. Wu T, Wu X, Wang HY, Chen L. Immune contexture defined by single cell technology for prognosis prediction and immunotherapy guidance in cancer. Cancer Commun (Lond) (2019) 39(1):21. doi: 10.1186/s40880-019-0365-9

CrossRef Full Text | Google Scholar

11. Hu X, Zhang J, Wang J, Fu J, Li T, Zheng X, et al. Landscape of B cell immunity and related immune evasion in human cancers. Nat Genet (2019) 51(3):560–7. doi: 10.1038/s41588-018-0339-x

CrossRef Full Text | Google Scholar

12. Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature (2020) 577(7791):549–55. doi: 10.1038/s41586-019-1922-8

CrossRef Full Text | Google Scholar

13. Garaud S, Buisseret L, Solinas C, Gu-Trantien C, de Wind A, Van den Eynden G, et al. Tumor infiltrating B-cells signal functional humoral immune responses in breast cancer. JCI Insight (2019) 5(18):e129641. doi: 10.1172/jci.insight.129641

CrossRef Full Text | Google Scholar

14. Chung W, Eum HH, Lee HO, Lee KM, Lee HB, Kim KT, et al. Single-cell RNA-seq enables comprehensive tumour and immune cell profiling in primary breast cancer. Nat Commun (2017) 8:15081. doi: 10.1038/ncomms15081

CrossRef Full Text | Google Scholar

15. Hu Q, Hong Y, Qi P, Lu G, Mai X, Xu S, et al. Atlas of breast cancer infiltrated B-lymphocytes revealed by paired single-cell RNA-sequencing and antigen receptor profiling. Nat Commun (2021) 12(1):2186. doi: 10.1038/s41467-021-22300-2

CrossRef Full Text | Google Scholar

16. Denkert C, von Minckwitz G, Darb-Esfahani S, Lederer B, Heppner BI, Weber KE, et al. Tumour-infiltrating lymphocytes and prognosis in different subtypes of breast cancer: a pooled analysis of 3771 patients treated with neoadjuvant therapy. Lancet Oncol (2018) 19(1):40–50. doi: 10.1016/S1470-2045(17)30904-X

CrossRef Full Text | Google Scholar

17. Mittendorf EA, Philips AV, Meric-Bernstam F, Qiao N, Wu Y, Harrington S, et al. PD-L1 expression in triple-negative breast cancer. Cancer Immunol Res (2014) 2(4):361–70. doi: 10.1158/2326-6066.CIR-13-0127

CrossRef Full Text | Google Scholar

18. Gatalica Z, Snyder C, Maney T, Ghazalpour A, Holterman DA, Xiao N, et al. Programmed cell death 1 (PD-1) and its ligand (PD-L1) in common cancers and their correlation with molecular cancer type. Cancer Epidemiol Biomarkers Prev (2014) 23(12):2965–70. doi: 10.1158/1055-9965.EPI-14-0654

CrossRef Full Text | Google Scholar

19. Luen S, Virassamy B, Savas P, Salgado R, Loi S. The genomic landscape of breast cancer and its interaction with host immunity. Breast (2016) 29:241–50. doi: 10.1016/j.breast.2016.07.015

CrossRef Full Text | Google Scholar

20. Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science (2015) 348(6230):69–74. doi: 10.1126/science.aaa4971

CrossRef Full Text | Google Scholar

21. Ribas A, Wolchok JD. Cancer immunotherapy using checkpoint blockade. Science (2018) 359(6382):1350–5. doi: 10.1126/science.aar4060

CrossRef Full Text | Google Scholar

22. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med (2012) 366(26):2443–54. doi: 10.1056/NEJMoa1200690

CrossRef Full Text | Google Scholar

23. Nikoo M, Rabiee F, Mohebbi H, Eghbalifard N, Rajabi H, Yazdani Y, et al. Nivolumab plus ipilimumab combination therapy in cancer: Current evidence to date. Int Immunopharmacol (2023) 117:109881. doi: 10.1016/j.intimp.2023.109881

CrossRef Full Text | Google Scholar

24. Li T, Wang X, Niu M, Wang M, Zhou J, Wu K, et al. Bispecific antibody targeting TGF-beta and PD-L1 for synergistic cancer immunotherapy. Front Immunol (2023) 14:1196970. doi: 10.3389/fimmu.2023.1196970

CrossRef Full Text | Google Scholar

25. Mortezaee K, Majidpoor J. Mechanisms of CD8(+) T cell exclusion and dysfunction in cancer resistance to anti-PD-(L)1. BioMed Pharmacother (2023) 163:114824. doi: 10.1016/j.biopha.2023.114824

CrossRef Full Text | Google Scholar

26. Yi M, Wu Y, Niu M, Zhu S, Zhang J, Yan Y, et al. Anti-TGF-beta/PD-L1 bispecific antibody promotes T cell infiltration and exhibits enhanced antitumor activity in triple-negative breast cancer. J Immunother Cancer (2022) 10(12):e005543. doi: 10.1136/jitc-2022-005543

CrossRef Full Text | Google Scholar

27. Ramakrishnan R, AsSudani D, Nagaraj S, Hunter T, Cho HI, Antonia S, et al. Chemotherapy enhances tumor cell susceptibility to CTL-mediated killing during cancer immunotherapy in mice. J Clin Invest (2010) 120(4):1111–24. doi: 10.1172/JCI40269

CrossRef Full Text | Google Scholar

28. Singh M, Al-Eryani G, Carswell S, Ferguson JM, Blackburn J, Barton K, et al. High-throughput targeted long-read single cell sequencing reveals the clonal and transcriptional landscape of lymphocytes. Nat Commun (2019) 10(1):3120. doi: 10.1038/s41467-019-11049-4

CrossRef Full Text | Google Scholar

29. Scott AM, Wolchok JD, Old LJ. Antibody therapy of cancer. Nat Rev Cancer (2012) 12(4):278–87. doi: 10.1038/nrc3236

CrossRef Full Text | Google Scholar

30. Rossetti RAM, Lorenzi NPC, Yokochi K, Rosa M, Benevides L, Margarido PFR, et al. B lymphocytes can be activated to act as antigen presenting cells to promote anti-tumor responses. PLoS One (2018) 13(7):e0199034. doi: 10.1371/journal.pone.0199034

CrossRef Full Text | Google Scholar

31. Zheng X, Jiang K, Xiao W, Zeng D, Peng W, Bai J, et al. CD8(+) T cell/cancer-associated fibroblast ratio stratifies prognostic and predictive responses to immunotherapy across multiple cancer types. Front Immunol (2022) 13:974265. doi: 10.3389/fimmu.2022.974265

CrossRef Full Text | Google Scholar

32. Uzhachenko RV, Shanker A. CD8(+) T lymphocyte and NK cell network: circuitry in the cytotoxic domain of immunity. Front Immunol (2019) 10:1906. doi: 10.3389/fimmu.2019.01906

CrossRef Full Text | Google Scholar

33. Zhao J, Shi Y, Cao G. The application of single-cell RNA sequencing in the inflammatory tumor microenvironment. Biomolecules (2023) 13(2):344. doi: 10.3390/biom13020344

CrossRef Full Text | Google Scholar

34. Baik JY, Liu Z, Jiao D, Kwon HJ, Yan J, Kadigamuwa C, et al. ZBP1 not RIPK1 mediates tumor necroptosis in breast cancer. Nat Commun (2021) 12(1):2666. doi: 10.1038/s41467-021-23004-3

CrossRef Full Text | Google Scholar

35. Fridman WH, Pages F, Sautes-Fridman C, Galon J. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer (2012) 12(4):298–306. doi: 10.1038/nrc3245

CrossRef Full Text | Google Scholar

36. Zhang C, Ge C. A Simple Competing Endogenous RNA Network Identifies Novel mRNA, miRNA, and lncRNA Markers in Human Cholangiocarcinoma. BioMed Res Int (2019) 2019:3526407. doi: 10.1155/2019/3526407

CrossRef Full Text | Google Scholar

37. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. Science (2015) 347(6220):1260419. doi: 10.1126/science.1260419

CrossRef Full Text | Google Scholar

38. Mei Y, Chen MM, Liang H, Ma L. A four-gene signature predicts survival and anti-CTLA4 immunotherapeutic responses based on immune classification of melanoma. Commun Biol (2021) 4(1):383. doi: 10.1038/s42003-021-01911-x

CrossRef Full Text | Google Scholar

39. Zhong H, Wang J, Zhu Y, Shen Y. Comprehensive analysis of a nine-gene signature related to tumor microenvironment in lung adenocarcinoma. Front Cell Dev Biol (2021) 9:700607. doi: 10.3389/fcell.2021.700607

CrossRef Full Text | Google Scholar

40. Castro M, Grau L, Puerta P, Gimenez L, Venditti J, Quadrelli S, et al. Multiplexed methylation profiles of tumor suppressor genes and clinical outcome in lung cancer. J Transl Med (2010) 8:86. doi: 10.1186/1479-5876-8-86

CrossRef Full Text | Google Scholar

41. McCullough LE, Chen J, Cho YH, Khankari NK, Bradshaw PT, White AJ, et al. Modification of the association between recreational physical activity and survival after breast cancer by promoter methylation in breast cancer-related genes. Breast Cancer Res (2017) 19(1):19. doi: 10.1186/s13058-017-0811-z

CrossRef Full Text | Google Scholar

42. Lehmann U, Berg-Ribbe I, Wingen LU, Brakensiek K, Becker T, Klempnauer J, et al. Distinct methylation patterns of benign and Malignant liver tumors revealed by quantitative methylation profiling. Clin Cancer Res (2005) 11(10):3654–60. doi: 10.1158/1078-0432.CCR-04-2462

CrossRef Full Text | Google Scholar

43. Evron E, Umbricht CB, Korz D, Raman V, Loeb DM, Niranjan B, et al. Loss of cyclin D2 expression in the majority of breast cancers is associated with promoter hypermethylation. Cancer Res (2001) 61(6):2782–7.

Google Scholar

44. Evron E, Dooley WC, Umbricht CB, Rosenthal D, Sacchi N, Gabrielson E, et al. Detection of breast cancer cells in ductal lavage fluid by methylation-specific PCR. Lancet (2001) 357(9265):1335–6. doi: 10.1016/S0140-6736(00)04501-3

CrossRef Full Text | Google Scholar

45. Hung CS, Wang SC, Yen YT, Lee TH, Wen WC, Lin RK. Hypermethylation of CCND2 in lung and breast cancer is a potential biomarker and drug target. Int J Mol Sci (2018) 19(10):3096. doi: 10.3390/ijms19103096

CrossRef Full Text | Google Scholar

46. Gross JA, Dillon SR, Mudri S, Johnston J, Littau A, Roque R, et al. TACI-Ig neutralizes molecules critical for B cell development and autoimmune disease. impaired B cell maturation in mice lacking BLyS. Immunity (2001) 15(2):289–302. doi: 10.1016/S1074-7613(01)00183-2

CrossRef Full Text | Google Scholar

47. Smulski CR, Eibel H. BAFF and BAFF-receptor in B cell selection and survival. Front Immunol (2018) 9:2285. doi: 10.3389/fimmu.2018.02285

CrossRef Full Text | Google Scholar

48. Salvucci O, Bouchard A, Baccarelli A, Deschenes J, Sauter G, Simon R, et al. The role of CXCR4 receptor expression in breast cancer: a large tissue microarray study. Breast Cancer Res Treat (2006) 97(3):275–83. doi: 10.1007/s10549-005-9121-8

CrossRef Full Text | Google Scholar

49. Chu QD, Panu L, Holm NT, Li BD, Johnson LW, Zhang S. High chemokine receptor CXCR4 level in triple negative breast cancer specimens predicts poor clinical outcome. J Surg Res (2010) 159(2):689–95. doi: 10.1016/j.jss.2008.09.020

CrossRef Full Text | Google Scholar

50. Chatterjee S, Behnam Azad B, Nimmagadda S. The intricate role of CXCR4 in cancer. Adv Cancer Res (2014) 124:31–82. doi: 10.1016/B978-0-12-411638-2.00002-1

CrossRef Full Text | Google Scholar

51. Zhang Z, Zhang Y, Xia S, Kong Q, Li S, Liu X, et al. Gasdermin E suppresses tumour growth by activating anti-tumour immunity. Nature (2020) 579(7799):415–20. doi: 10.1038/s41586-020-2071-9

CrossRef Full Text | Google Scholar

52. Zhou Y, Zheng J, Bai M, Gao Y, Lin N. Effect of pyroptosis-related genes on the prognosis of breast cancer. Front Oncol (2022) 12:948169. doi: 10.3389/fonc.2022.948169

CrossRef Full Text | Google Scholar

53. Zhong Y, Peng F, Luo X, Wang X, Yang B, Tang X, et al. A pyroptosis-related gene signature for prognostic and immunological evaluation in breast cancer. Front Oncol (2022) 12:964508. doi: 10.3389/fonc.2022.964508

CrossRef Full Text | Google Scholar

54. Bohgaki M, Tsukiyama T, Nakajima A, Maruyama S, Watanabe M, Koike T, et al. Involvement of Ymer in suppression of NF-kappaB activation by regulated interaction with lysine-63-linked polyubiquitin chain. Biochim Biophys Acta (2008) 1783(5):826–37. doi: 10.1016/j.bbamcr.2007.09.006

CrossRef Full Text | Google Scholar

55. Xu Y, Wang Y, Liang L, Song N. Single-cell RNA sequencing analysis to explore immune cell heterogeneity and novel biomarkers for the prognosis of lung adenocarcinoma. Front Genet (2022) 13:975542. doi: 10.3389/fgene.2022.975542

CrossRef Full Text | Google Scholar

56. Song B, Shen S, Fu S, Fu J. HSPA6 and its role in cancers and other diseases. Mol Biol Rep (2022) 49(11):10565–77. doi: 10.1007/s11033-022-07641-5

CrossRef Full Text | Google Scholar

57. Shen S, Wei C, Fu J. RNA-sequencing reveals heat shock 70-kDa protein 6 (HSPA6) as a novel thymoquinone-upregulated gene that inhibits growth, migration, and invasion of triple-negative breast cancer cells. Front Oncol (2021) 11:667995. doi: 10.3389/fonc.2021.667995

CrossRef Full Text | Google Scholar

58. Kilbas PO, Akcay IM, Doganay GD, Arisan ED. Bag-1 silencing enhanced chemotherapeutic drug-induced apoptosis in MCF-7 breast cancer cells affecting PI3K/Akt/mTOR and MAPK signaling pathways. Mol Biol Rep (2019) 46(1):847–60. doi: 10.1007/s11033-018-4540-x

CrossRef Full Text | Google Scholar

59. Mahmood SF, Gruel N, Nicolle R, Chapeaublanc E, Delattre O, Radvanyi F, et al. PPAPDC1B and WHSC1L1 are common drivers of the 8p11-12 amplicon, not only in breast tumors but also in pancreatic adenocarcinomas and lung tumors. Am J Pathol (2013) 183(5):1634–44. doi: 10.1016/j.ajpath.2013.07.028

CrossRef Full Text | Google Scholar

60. Bernard-Pierrot I, Gruel N, Stransky N, Vincent-Salomon A, Reyal F, Raynal V, et al. Characterization of the recurrent 8p11-12 amplicon identifies PPAPDC1B, a phosphatase protein, as a new therapeutic target in breast cancer. Cancer Res (2008) 68(17):7165–75. doi: 10.1158/0008-5472.CAN-08-1360

CrossRef Full Text | Google Scholar