
==== Front
J Pharm Anal
J Pharm Anal
Journal of Pharmaceutical Analysis
2095-1779
2214-0883
Xi'an Jiaotong University

S2095-1779(24)00072-8
10.1016/j.jpha.2024.100975
100975
Original Article
Dissection of triple-negative breast cancer microenvironment and identification of potential therapeutic drugs using single-cell RNA sequencing analysis
Cheng Weilun a1
Mi Wanqi b1
Wang Shiyuan c1
Wang Xinran b
Jiang Hui b
Chen Jing b
Yang Kaiyue b
Jiang Wenqi b
Ye Jun yelinghao@imm.ac.cn
c⁎⁎⁎
Guo Baoliang baoliangguo2013@163.com
a⁎⁎
Zhang Yunpeng zhangyp@hrbmu.edu.cn
b⁎
a Department of General Surgery, The Second Affiliated Hospital of Harbin Medical University, Harbin, 150001, China
b College of Bioinformatics Science and Technology, Harbin Medical University, Harbin, 150081, China
c State Key Laboratory of Bioactive Substance and Function of Natural Medicines, Institute of Materia Medica, Chinese Academy of Medical Sciences & Peking Union Medical College, Beijing, 100050, China
⁎ Corresponding author. zhangyp@hrbmu.edu.cn
⁎⁎ Corresponding author. baoliangguo2013@163.com
⁎⁎⁎ Corresponding author. yelinghao@imm.ac.cn
1 These authors contributed equally to this work.

02 4 2024
8 2024
02 4 2024
14 8 1009752 2 2024
28 3 2024
30 3 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Breast cancer remains a leading cause of mortality in women worldwide. Triple-negative breast cancer (TNBC) is a particularly aggressive subtype characterized by rapid progression, poor prognosis, and lack of clear therapeutic targets. In the clinic, delineation of tumor heterogeneity and development of effective drugs continue to pose considerable challenges. Within the scope of our study, high heterogeneity inherent to breast cancer was uncovered based on the landscape constructed from both tumor and healthy breast tissue samples. Notably, TNBC exhibited significant specificity regarding cell proliferation, differentiation, and disease progression. Significant associations between tumor grade, prognosis, and TNBC oncogenes were established via pseudotime trajectory analysis. Consequently, we further performed comprehensive characterization of the TNBC microenvironment. A crucial epithelial subcluster, E8, was identified as highly malignant and strongly associated with tumor cell proliferation in TNBC. Additionally, epithelial-mesenchymal transition (EMT)-associated fibroblast and M2 macrophage subclusters exerted an influence on E8 through cellular interactions, contributing to tumor growth. Characteristic genes in these three cluster cells could therefore serve as potential therapeutic targets for TNBC. The collective findings provided valuable insights that assisted in the screening of a series of therapeutic drugs, such as pelitinib. We further confirmed the anti-cancer effect of pelitinib in an orthotopic 4T1 tumor-bearing mouse model. Overall, our study sheds light on the unique characteristics of TNBC at single-cell resolution and the crucial cell types associated with tumor cell proliferation that may serve as potent tools in the development of effective anti-cancer drugs.

Graphical abstract

Image 1

Highlights

• A common cluster of proliferating malignant epitheliums was identified in TNBC.

• TPX2, CDCA8, PLK1, UBE2S, RRM2, and TK1 are TNBC independent risk factors.

• EMT fibroblast and M2 macrophage promote aforementioned proliferating malignant epitheliums for tumor growth in TNBC.

• Drugs like pelitinib inhibit TNBC by suppressing multiple cell clusters in tumor microenvironment.

Keywords

Triple-negative breast cancer
Single-cell RNA sequencing
Epithelial-mesenchymal transition
Macrophage polarization
Therapeutic drugs
Bioinformatics
==== Body
pmc1 Introduction

Breast cancer is the most common malignant tumor type in women worldwide. According to the latest global statistical report, over 2.26 million new cases of breast cancer are diagnosed annually, accounting for 11.7% of all tumor cases. Breast cancer is responsible for more than 680,000 deaths on a yearly basis and poses a serious threat to the lives and health of patients [1].

Perou et al. [2] were the first research group to analyze microarray data of tumor tissues and classify breast cancer at the molecular level. In the expert consensus of St. Gallen International Expert Consensus on the Primary Therapy of Early Breast Cancer in 2011 [3], clinical pathological indicators were introduced as the typing criteria given their affordability and accessibility. Breast cancer is mainly divided into four subtypes (i.e., luminal A, luminal B, human epidermal growth factor receptor 2 (HER2) positive, and triple-negative breast cancer (TNBC)) based on evaluation of a set of molecular markers, including estrogen receptor (ER), progesterone receptor (PR), HER2, and Ki-67. Treatment strategies are selected according to molecular subtype, highlighting the significant variability of breast cancer.

Among the breast cancer subtypes identified to date, TNBC exhibits greater complexity at the molecular and cellular levels and thus warrants further in-depth investigation. Typically, TNBCs are extremely invasive with significant metastatic potential, high recurrence rates, and poor clinical outcomes [4,5]. Since TNBC does not respond to hormone receptor- or HER2-targeted drugs, chemotherapy currently remains the frontline accepted standard treatment [6]. In recent years, prognosis of TNBC patients has improved to some extent, owing to the implementation of personalized and precision therapy with small-molecule drugs in conjunction with novel immunotherapeutic agents [7]. Nevertheless, a proportion of patients continue to display insensitivity or resistance to existing treatments, highlighting the necessity for stable, reliable, and effective measures that encompass a broader patient population.

In this study, we extensively explored the mechanisms underlying breast cancer development using single-cell RNA sequencing (scRNA-seq) technology [8,9], which offers precise insights via analysis of gene expression profiles at single-cell resolution. Breast cancer is a heterogeneous disease co-constituted by tumor cells and their surrounding microenvironment, and therefore response to treatment varies from patient to patient. This tumor heterogeneity is often difficult to evaluate by classic clinical parameters such as histopathological features, tumor grading, and lymph node metastases. Bulk sequencing results only provide an average expression profile of tumor tissue, and true gene expression signals from rare cell clusters or cell types that drive tumorigenesis or influence drug resistance may be masked by the average gene expression profile. However, scRNA-seq solves this problem perfectly. The ability of single-cell technology to characterize the transcriptional profile of each cell has led to an unprecedented level of understanding of cellular heterogeneity. The scRNA-seq technique can be utilized to evaluate both tumor cells and the microenvironment [10,11] with the purpose of providing key insights into cell characteristics, tumor development, and drug response.

TNBC displays the most unique traits among the different subtypes of breast cancer. Early in 2018, Karaayvaz et al. [12] conducted an initial exploration of TNBC from a single-cell perspective, finding a subcluster of tumor cells with the largest scale of copy number variations (CNVs). Sebastian et al. [13] used scRNA-seq to identify the molecular characteristics of six cancer-associated fibroblast (CAF) subclusters in TNBC, and compared the heterogeneity of CAF in different tumor sites and healthy tissues.

Previous scRNA-seq studies of breast cancer tended to identify a particular group of cells and to look for specific clinical features associated with that cell type. In this study, we innovatively considered both the malignant epithelial cell cluster and other cells related to the malignant epithelium within TNBC microenvironment at single-cell resolution. All of these factors associated with poor prognosis were considered together as a whole and we calculated to produce several potential therapeutic drugs capable of acting simultaneously on multiple risk factors in the microenvironment, which could provide some new therapeutic regimens for the treatment of TNBC. We further identified pelitinib and other potential therapeutic drugs that could selectively target genes within TNBC microenvironment, offering novel strategies for the management of TNBC.

2 Materials and methods

2.1 Cell line and animals

The 4T1 murine TNBC line was a kind gift from Prof. Zhonggao Gao (Institute of Materia Medica, Peking Union Medical College, Beijing, China) and maintained in Dulbecco's modified Eagle medium (DMEM) containing 10% (V/V) fetal bovine serum (FBS; Gibco Life Technologies, Carlsbad, CA, USA) at 37 °C in a 5% CO2 humidified air incubator (Thermo Fisher Scientific Inc., Waltham, MA, USA).

Seven-week-old female BALB/c mice (body weight 18–20 g) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China). All animal procedures were carried out in strict accordance with the National Institutes of Health, USA and approved by the Animal Care and Welfare Committee of Institute of Materia Medica, Chinese Academy of Medical Sciences and Peking Union Medical College, China (Approval No.: 00004864).

2.2 Drug and antibodies

Pelitinib (Cat. No.: HY-32718) was obtained from MedChemExpress (Monmouth Junction, NJ, USA). Rabbit monoclonal antibodies for β-catenin (Cat. No.: #8480), vimentin (Cat. No.: #5741), snail (Cat. No.: #3879), E-cadherin (Cat. No.: #3195), PI3 kinase (Cat. No.: #4257), phospho-protein kinase B (Akt) (Ser473) (Cat. No.: #4060), pan-Akt (Cat. No.: #4685), β-actin (Cat. No.: #4970), and goat anti-rabbit secondary antibody (Cat. No.: #7074) were acquired from Cell Signaling Technology Inc. (Danvers, MA, USA). Rabbit polyclonal antibodies for zonula occludens-1 (ZO-1) (Cat. No.: 21773-1-AP), F4/80 (Cat. No.: 28463-1-AP), CD206 (Cat. No.: 18704-1-AP), and CD86 (Cat. No.: 13395-1-AP) were acquired from Proteintech (Wuhan, China).

2.3 Data processing

Breast single-cell 10x Genomics transcriptome sequencing data were obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database (GSE176078, GSE161529, and GSE180878). The 16 samples used for first-step analysis included ER+ (GSM4909298, GSM4909301, GSM4909303, and GSM4909304), HER2+ (GSM4909289, GSM4909291, GSM4909292, and GSM4909293), TNBC (GSM4909281, GSM4909282, GSM4909283, and GSM4909284), and healthy control (GSM4909263, GSM4909265, GSM4909266, and GSM4909268) groups. The 13 TNBC samples selected for further analysis included GSM5354517, GSM5354525, GSM5354528, GSM5354530, GSM5354531, GSM5354532, GSM5354534, GSM4909281, GSM4909283, GSM4909284, GSM4909285, GSM4909286, and GSM4909287. Fibroblasts obtained from three healthy volunteers with no history of breast cancer were used as control samples (GSM5474104, GSM5474105, and GSM5474106).

The Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset [14] containing bulk sequencing data and paired clinical detailed information on 1980 primary breast cancer samples was sourced from the cBioPortal database (https://www.cbioportal.org/) [15].

2.4 scRNA-seq data quality control and analysis

Research and analysis were based on R (ver. 4.2.2) conducted using RStudio (ver. 2023.03.0 + 386). Low-quality cells with fewer than 200 genes were filtered with Seurat package (ver. 4.2.0) [16]. Batch effects were eliminated using Harmony package (ver. 0.1.0) [17]. Uniform manifold approximation and projection (UMAP) using the Seurat package was used to reduce multi-dimensional information to two dimensions.

2.5 Single-sample gene set enrichment analysis

The Gene Set Variation Analysis (GSVA) package (ver. 1.40.1) was used for single-sample gene set enrichment analysis (ssGSEA) [18] of hallmark gene sets from Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). The ClusterProfiler (ver. 4.7.1.002) package [19] was employed for functional enrichment of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [20], with the aim of extracting biologically significant patterns and functions from a vast amount of gene data.

2.6 CNV analysis

Infer CNV (ver. 1.8.1) [21] and Copy number KAryotyping of Tumors (CopyKAT) (ver. 1.1.0) [22] algorithms were employed to calculate the distribution of genomic copy numbers for a single cell and detect high-confidence diploid cells. Most stromal and immune cells have diploid copy number profiles. By identifying aneuploid cells in single-cell transcriptome data, tumor cells were identified and subclone structures within the tumor were revealed.

2.7 Pseudotime trajectory analysis

Pseudotime analysis was performed using Monocle2 (ver. 2.22.0) and Monocle3 (ver. 1.0.0) packages [23] with the aim of constructing dramatic cell transformations among different clusters over time. Dynamic changes in gene expression were captured based on pseudotime trajectory mapping.

2.8 Survival analysis

The Survival (ver. 3.2–11) package was employed for analyzing overall survival of patients in the METABRIC database. Patients were classified into high-expression and low-expression groups based on specific genes. Kaplan-Meier survival curves were generated for both groups.

2.9 Analysis of tumor purity

The Estimate of STromal and Immune cells in MAlignant Tumor tissues from Expression data (ESTIMATE) (ver. 1.0.13) algorithm [24] was used to identify tumor cells using TNBC transcriptional profiles. Following filtering of stromal and immune signatures from the public database, these two scores were applied to generate an ESTIMATE score, which was utilized to calculate tumor purity.

2.10 Gene co-expression analysis

High dimensional weighted gene co-expression network analysis (hdWGCNA) (ver. 0.2.01) [25] was conducted to identify gene modules co-expressed in all TNBC epithelial cells and explore core genes in the modules.

2.11 Cell-cell communication

The CellPhoneDB (ver. 4.0.0) program [26] was utilized to predict interactions between cell clusters in TNBC by analyzing the expression patterns of proteins involved in communication between different cell types.

2.12 Construction of protein-protein interaction networks

Information on protein interactions was obtained from the STRING database (https://string-db.org/). The characteristic genes of crucial subclusters of TNBC and genes involved in ligand-receptor interactions were mapped into an interaction network, which was visualized using Cytoscape (ver. 3.9.0) [27].

2.13 Drug repurposing

A Single-cell Guided Pipeline to Aid Repurposing of Drugs (ASGARD) (ver. 1.0.0) package [28] was used to score each drug for repositioning by simultaneously considering the situation of multiple cell clusters in scRNA-seq. Key cell subclusters in TNBC were analyzed using this technique, with the aim of identifying potential therapeutic drugs.

2.14 Tumor model and treatment

The 5 × 105 4T1 cells in 100 μL of phosphate-buffered saline (PBS) were injected under the right third pair of fat pads. The mice were equally and randomly divided into four groups (n = 8) when the tumor volume was 30–50 mm3: pelitinib low-dose group (5 mg/kg/day), pelitinib medium-dose group (10 mg/kg/day), pelitinib high-dose group (20 mg/kg/day) with oral administration, and control group. All the animals were treated for a period of 14 days. Tumor volume was calculated by the formula of (L × W2)/2 with a digital caliper every two days, where L is the longest tumor dimension (mm) and W is the shortest tumor dimension (mm). The toxicity was monitored by observing behavior and measuring body weight every day. At the termination of the experiment, all animals were euthanized by cervical dislocation. Tumor issues, as well as hearts, livers, spleens, lungs, and kidneys, were excised and fixed in 4% paraformaldehyde overnight or stored at −80 °C until further analysis.

2.15 Histopathology and immunohistochemistry (IHC) staining

Tumors, hearts, livers, spleens, lungs, and kidneys were harvested, fixed in 4% paraformaldehyde, embedded in paraffin, sectioned, and stained with hematoxylin and eosin (H&E). Sections for IHC were deparaffinized, rehydrated, and subjected to antigen retrieval. Sections were incubated overnight at 4 °C with primary antibodies F4/80, CD206, and CD86. In the next day, tissue sections were washed and incubated with secondary antibody at room temperature away from light for 1 h. After 3,3-diaminobenzidine tetrahydrochloride (DAB) staining, the IHC H-score [29] was calculated based on the intensity and area of immunohistochemical staining by ImageJ (ver.1.54h) software.

2.16 Western blot analysis

The tumor tissue samples were homogenized in ice-cold lysis buffer containing complete tablets, mini EASYpack (Roche, Basel, Switzerland), and 1 mmol/L phenylmethanesulfonyl fluoride (Solarbio, Beijing, China) for 3 min. The homogenate was centrifuged at 4 °C with a speed of 12,000 r/min for 15 min and the supernatant was collected. The protein concentration was determined using the bicinchoninic acid (BCA) assay (Thermo Fisher Scientific Inc.). Next, the proteins were separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to a polyvinylidene difluoride (PVDF) membrane. The membranes were blocked with 5% bovine serum albumin (BSA) for 1.5 h at room temperature, then incubated with the following primary antibodies: β-catenin, vimentin, Snail, E-cadherin, ZO-1, PI3 kinase, phospho-Akt (Ser473), pan-Akt, and β-actin overnight at 4 °C. In the next day, the membranes were rinsed using tris-buffered saline with Tween 20 (TBST) five times and incubated with horseradish peroxidase (HRP)-conjugated anti-rabbit secondary antibody (1:5000) for 1 h at room temperature. The protein expression bands were detected using enhanced chemiluminescence (Easybio, Beijing, China) and visualized using a Tanon-4600SF chemiluminescence imager (Tanon Science & Technology Co., Ltd., Shanghai, China).

2.17 Statistical analysis

Data analysis was calculated with GraphPad Prism (version 10.1.2(324)). Significance differences between different samples was revealed with Student's t-test, or one-way analysis of variance (ANOVA). P < 0.05 was considered to be statistical significance.

3 Results

3.1 Determination of cellular heterogeneity of the different breast cancer subtypes

To obtain an overview of global cell characteristics within breast tissue samples, we acquired scRNA-seq data (GSE161529) encompassing 16 evenly distributed samples in terms of molecular subtype (ER+ (luminal A and luminal B subtypes), HER2+, TNBC, and healthy breast tissues derived from individuals with no family history of breast cancer who underwent reduction mammoplasty). A total of 72,501 cells were retained for subsequent analysis after filtering and quality control to exclude abnormal cells.

UMAP dimensionality reduction clustering was conducted to visualize the cellular landscape of breast cancer molecular subtypes and healthy controls. Canonical cell markers were employed to annotate cell clusters, including epithelial cells (EPCAM), immune cells (CD45), fibroblasts (PDPN), endothelial cells (PECAM1), and perithelial cells (RGS5) [[30], [31], [32], [33], [34], [35]] (Figs. 1A−D). Significant differences in cell distribution and proportion were observed among the four categories. Perithelial cells were predominantly detected in healthy tissues while only trace amounts were present in tumor tissues. Conversely, healthy samples contained a limited number of immune cells while breast cancer samples exhibited immune cell infiltration. Meanwhile, fibroblasts were significantly more abundant in healthy samples (Figs. 1E and F). Comparison of the cell composition across different samples not only disclosed significant differences among the subtypes but also provided evidence of cellular heterogeneity within samples of the same subtype. The most significant differentially expressed genes (DEGs) identified in all the clusters were consistent with prior findings, validating the accuracy of our clustering results (Figs. 1G and H).Fig. 1 Expression profiling of 72,501 single cells in breast cancer subtypes and healthy controls. (A–D) Uniform manifold approximation and projection (UMAP) plots of 72,501 cells in 16 samples analyzed via single-cell RNA sequencing (scRNA-seq) from breast cancer subtypes and healthy controls: estrogen receptor (ER)+ breast cancer (A), human epidermal growth factor receptor 2 (HER)+ breast cancer (B), triple-negative breast cancer (TNBC) (C), and healthy group (D). (E) Relative proportions of cell types highlighting clear differences across clinical subtypes. (F) Variable proportions of the same cell type in different groups. (G) Average gene expression levels of cellular markers corresponding to different cell types. (H) Representative markers of cell types among the subtypes. (I) Heatmap showing genes of different cell types among subtypes enriched in different functions. NES: normalized enrichment scores.

Fig. 1

The ssGSEA method, a tool applied to analyze the function of each cell type based on characteristic genes, revealed distinct patterns in epithelial cells among the breast cancer subtypes. Estrogen response-related pathways were highly enriched in ER+ breast cancer and oxidative phosphorylation and peroxisome-related pathways in HER2+ breast cancer [36,37]. In TNBC, we observed strong enrichment of pathways related to the cell cycle, particularly e2f_targets and g2m_checkpoint [38], indicative of an active cellular state. In contrast to tumor cells, those from healthy epithelium demonstrated functional links to androgen response and protein secretion [39]. Cells from tumor tissues showed apparent separation while those from healthy tissues were aggregated, suggesting significant heterogeneity between tumor and healthy cell groups (Fig. 1I).

3.2 TNBC specificity based on highly proliferative and differentiated epithelial cells

Analysis of epithelial cells of the different subtypes revealed a unique feature of TNBC. Initially, epithelial cells from each tumor subtype were isolated and subsequently integrated with healthy epithelial cells. After re-clustering of all epithelial cells, distinct combinatorial patterns across the different subtypes emerged. In the ER+ group, clusters 0, 1, and 3 were mainly derived from healthy samples whereas clusters 2 and 4 predominantly originated from tumor samples. In the HER2+ group, tumor samples were the principal source of cluster 4/5/6, with a minor contribution from cluster 0. In the TNBC group, clusters 0 and 2 were predominantly obtained from tumor samples (Figs. 2A−C). Next, we examined highly expressed genes specific for each subcluster along with the major enriched functions within each cancer subtype, including response to stress, unfolded protein, and nuclear division (Figs. 2D and E). Additionally, potential aneuploid variations in malignant cells were investigated. Since malignant cells are typically aneuploid variants, we evaluated each cell type by gene CNV. Our data showed that the malignant cells identified were located in subclusters dominated by tumor cells as expected. (Fig. S1).Fig. 2 Epithelial cell characteristics of breast cancer subtypes. (A–C) Uniform manifold approximation and projection (UMAP) plots of emergence of tumor and healthy epithelial cells: estrogen receptor (ER)+ and healthy (A), human epidermal growth factor receptor 2 (HER2)+ and healthy (B), and triple-negative breast cancer (TNBC) and healthy (C). (D) Bubble plots of marker genes expressed in the respective subclusters among ER+, HER2+, and TNBC subtypes. (E) Main pathways for enrichment of marker genes in each subcluster. (F−H) Cluster 2 in ER+ (F), clusters 4 and 6 in HER2+ (G), and cluster 2 in TNBC (H) with higher proliferation and differentiation scores. (I) Violin plots displaying marker genes of cluster 2 in ER+, clusters 4 and 6 in HER2+, and cluster 2 in TNBC. (J) Venn diagram showing overlap of marker genes in the three clusters. (K) Gene Ontology (GO) enrichment of genes in the above three proliferative and differentiated epithelial clusters. (L) Scores of the validation set with 33 TNBC genes and 25 overlapping genes in the ER+ and HER2+ groups. PI3K: phosphoinositide 3-kinase; Akt: protein kinase B; mTOR: mechanistic target of rapamycin; mTORC1: mechanistic target of rapamycin complex 1; K-Ras: Kirsten rats arcomaviral oncogene homolog; TNF-α: tumor necrosis factor-α; NF-κB: nuclear factor-kappaB; IL-2: interleukin-2; STAT5: signal transducers and activators of transduction 5; JAK: Janus kinase; IFN-α: interferon-α; TGF-β: transforming growth factor-β; NES: normalized enrichment scores.

Fig. 2

GSVA was applied to calculate the proliferation and differentiation scores for each cell type (Figs. 2F−H). Among the ER+ subtypes, cluster 2, mainly composed of cells from cancer patients, displayed the highest proliferation and differentiation scores (Fig. 2A), and was defined as a malignant cluster. Similarly, clusters 4 and 6 of HER2+ and cluster 2 of TNBC subtypes exhibited high proliferation and differentiation potential. These clusters also showed significant overlap with aneuploid cells identified via CNV, suggesting a strong correlation with the malignant phenotype.

Subsequently, malignant tumor subcluster-specific DEGs were examined (Fig. 2I). Interestingly, we observed a number of similarities in the characteristic genes of ER+ and HER2+ cells. Among the 32 and 39 characteristic genes identified in the ER+ and HER2+ subtypes, respectively, 25 showed overlap. However, the pattern of differential gene expression in TNBC was distinct from those of the other two subtypes (Fig. 2J).

GO analysis revealed that characteristic genes of tumor cells from all three subtypes were functionally enriched in the following terms: cellular process, biological process, and regulation of the cell cycle. Interestingly, the TNBC subtype further exhibited unique enrichment in assembly and arrangement of cellular components and cellular metabolic processes [40,41], while ER+ and HER2+ isoform-specific genes shared similar functional expression patterns in relation to the regulation of cell communication, cell differentiation, and cell death [42,43] (Fig. 2K).

To ensure the robustness of our data on differences between TNBC and the other two subtypes, an additional scRNA transcriptome dataset (GSE176078) was acquired as a validation set with a view to establishing the reliability of our conclusions. Specifically, a TNBC signature comprising the 33 DEGs within cluster 2 was defined and the gene set applied to calculate the enrichment score of the validation set. TNBC samples exhibited significantly higher enrichment scores in the validation set. In parallel, the validation set was scored using the 25 genes shared by the ER+ and HER2+ subtypes. Similar scores were obtained for ER+ and HER2+ samples, which were statistically different from TNBC (Fig. 2L), further validating the specificity of TNBC and its unique gene expression pattern.

3.3 Genes involved in TNBC epithelial cell differentiation play a crucial role in disease grade and prognosis

We further investigated the process of epithelial cell differentiation with the aim of uncovering crucial factors underlying the transition from healthy epithelial cells to cancerous cells. After integration of healthy epithelium, pseudotime trajectory analysis was conducted for each subtype (Figs. 3A−C and S2). The pseudotime differentiation trajectories indicated that epithelial cells from cancer and healthy groups merged into one or more states. In these mixed states, cells from benign and malignant origins presented shared characteristics, suggesting that healthy cells begin to exhibit early stages of malignant potential at this time. In view of this finding, we further focused on oncogenes within these cell clusters that could potentially underlie the differentiation and transition process from normal epithelial cells into a malignant state. Pseudotime trajectories including genes with significant variations in expression were visualized, which showed distinct expression patterns in the initial and final states of differentiation (Figs. 3D−F). Trajectories of development are illustrated with arrows in Fig. 3A. The key genes responsible for the development of malignancy were identified. State 2 in ER+ and states 2 and 3 in HER2+ subtypes were established as transitional cell states. For the TNBC subtype, trajectories were significantly more complex. Accordingly, we grouped states 3 and 4 and states 5 and 6 together for analysis.Fig. 3 Profiles of significant differentially expressed genes (DEGs) in epithelial cells. (A–C) Epithelial cell transition along the pseudotime trajectory, colored by cell state (upper) and cell source (below): estrogen receptor (ER)+ (A), human epidermal growth factor receptor 2 (HER2)+ (B), and triple-negative breast cancer (TNBC) (C). 1–3: branching points in pseudotime trajectory. (D–F) Highly variable gene expression changes in ER+ (D), HER2+ (E), and TNBC (F) total trajectories. (G, H) Venn diagram displaying overlap of upregulated (G) and downregulated (H) genes in ER+ (state 2), ER+ (states 2 and 3), and TNBC (states 3 and 4 and states 5 and 6). (I, J) Significantly enriched pathways of upregulated (I) and downregulated (J) genes in Gene Ontology (GO) terms. (K) Expression of TNBC upregulated genes is increased with tumor grade in the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) database. (L) Expression of TNBC downregulated genes shows no significant differences among tumor grades. (M, N) Poorer prognosis of TNBC patients in the upregulated gene high-expression group (M) and downregulated gene-low expression group (N). ∗∗P < 0.01 and ∗∗∗P < 0.001. NS: not significant.

Fig. 3

Comparison of oncogenes among the subtypes revealed more upregulated and downregulated genes in TNBC relative to the ER+ and HER2+ subtypes (Figs. 3G and H and Table S1). Although most of the upregulated genes were different for each subtype, we observed some overlap in the GO enrichment results. The HER2+ subtype had fewer upregulated genes and lacked significant enrichment passways. The overlapping functions encompassed processes such as cell differentiation and tissue development [44] (Fig. 3I). Conversely, downregulated genes showed higher overlap across subtypes. In TNBC, the downregulated genes were markedly enriched in multiple processes associated with stress or stimulus as well as immune responses [45] (Fig. 3J). Individual characteristics of the different subtypes affected various aspects of cell differentiation.

To confirm our results, bulk sequencing data from the METABRIC dataset were obtained for the evaluation of clinical information. Our analysis revealed an intriguing relationship between expression levels of these genes and grades assigned using the Nottingham scoring system for breast cancer, which took into account factors such as glandular duct formation, nuclear pleomorphism, and mitotic morphology [46]. The expression of upregulated genes in TNBC increased significantly with grade, with the highest expression detected in grade 3 tumors. In contrast, expression patterns of downregulated genes displayed no substantial variations across cancer grades (Figs. 3K and L). We further explored the impact of these genes on patient prognosis. As expected, patients in the high expression group of upregulated genes displayed poorer prognosis, which could be attributed to their function in promoting tumor transformation [47,48]. Conversely, longer overall survival was associated with the high expression group of downregulated genes, supporting the involvement of these genes in inhibiting or delaying the development of malignancy in cells (Figs. 3M and N).

3.4 Comprehensive evaluation of the TNBC tumor microenvironment (TME)

Considering the unique gene expression pattern compared to other subtypes, we comprehensively investigated the cellular characteristics of TNBC. To this end, 13 primary TNBC samples with no preoperative adjuvant treatment from scRNA-seq (GSE176078 and GSE161529) were employed for analysis. Within the TNBC landscape, 7 main clusters were identified, specifically, epithelial cells, T lymphocytes, B lymphocytes, plasma cells, myeloid cells, fibroblasts, and endothelial cells. Most cell types tended to cluster closely, apart from epithelial cells, which displayed a more disperse distribution, emphasizing the high heterogeneity within this cell population (Fig. 4A). Feature plots and violin plots were generated to visualize distinctive and characteristic genes, which were further utilized to annotate clusters. Our data revealed exclusive expression of KRT19 and EPCAM in epithelial cells and CD3E in T lymphocytes (Figs. 4B and C). Each cell type featured a distinct set of highly expressed genes, with variations in gene expression patterns across different cell types. The corresponding heat map effectively illustrated the utility of highly variable genes in cell type discrimination (Fig. 4D).Fig. 4 Landscape of 42,488 triple-negative breast cancer (TNBC) single cells. (A) Uniform manifold approximation and projection (UMAP) plots of 42,488 single cells analyzed via single-cell RNA sequencing (scRNA-seq) derived from 13 primary TNBC patients. (B, C) Feature plots (B) and violin plots (C) displaying expression of marker genes in epithelial cells (KRT19 and EPCAM), T lymphocytes (CD3E), B lymphocytes (MS4A1), plasma cells (IGHG1), myeloid cells (AIF1), fibroblasts (COL3A1), and endothelial cells (PLVAP). (D) Heatmap showing expression of top significant differentially expressed genes (DEGs) among the seven main cell types in TNBC. (E) Average enrichment of genes in each cell type by hallmark gene set. (F) Distribution of the seven cell types in each TNBC sample. (G, H) Circle plot (G) and heatmap (H) showing the intensity and quantity of ligand-receptor interactions between the seven cell types. IL-2: interleukin-2; STAT5: signal transducers and activators of transduction 5; K-Ras: Kirsten rats arcomaviral oncogene homolog; IFN-α: interferon-α; IL-6: interleukin-6; JAK: Janus kinase; TNF-α: tumor necrosis factor-α; NF-κB: nuclear factor-kappaB; PI3K: phosphoinositide 3-kinase; Akt: protein kinase B; mTOR: mechanistic target of rapamycin; mTORC1: mechanistic target of rapamycin complex 1; TGF-β: transforming growth factor-β; NES: normalized enrichment scores.

Fig. 4

In gene functional enrichment analysis using 50 hallmark gene sets, the seven cell types could be broadly categorized into three groups: epithelial, immune, and stromal cells. Epithelial cells were highly enriched in tumor-related functions, including Kirsten ratsarcoma viral oncogenehomolog (KRAS) signaling, glycolysis, and oxidative phosphorylation pathways [[49], [50], [51]], immune cells exhibited greater activity in interferon response pathways [52], and stromal cells were strongly linked to cell apoptosis and angiogenesis (Fig. 4E). Significant differences in cell composition were observed among all samples belonging to the TNBC subtype. A high proportion of epithelial cells was detected in patients 11 and 12, indicative of more pronounced tumor characteristics. Conversely, patients 9 and 10 contained higher numbers of immune cells, reflecting increased sensitivity to immunotherapy. These findings highlight the substantial heterogeneity within TNBC, which could explain the differences in treatment response among TNBC patients (Fig. 4F). We further investigated the ligand-receptor interactions between different cell types. Notably, despite accounting for a minority of all cells, stromal cells exhibited strongest interactions with other cell types (Figs. 4G and H).

3.5 Identification of crucial cells and genes of TNBC via holistic analysis

For identification of tumor cells, we evaluated CNVs across all cell types using the CopyKAT algorithm. Since changes in genomic DNA copy number constitute a pivotal difference between tumor and normal cells, aneuploid cells exhibiting substantial CNVs were categorized as tumor cells while diploid cells were regarded as benign (Fig. 5A). Upon classifying cells into tumor or non-tumor groups within the UMAP cluster diagram, approximately 98.75% of tumor cells were identified within the epithelial cell group. Non-tumor cells mainly included other cell types, with the appearance of a few epithelial cells, which formed distinct clusters away from tumor cells (Figs. 5B and C). Furthermore, when applying ESTIMATE to score all cells, the tumor purity score was the highest among epithelial cells, which validated the robustness of our dimensionality reduction clustering of TNBC cells and tumor cell recognition process (Fig. 5D).Fig. 5 Epithelial cells and their differentiation trajectories in triple-negative breast cancer (TNBC) patients. (A) Copy number variation (CNV) analysis distinguishing aneuploid and diploid cells. (B) Cells in uniform manifold approximation and projection (UMAP) diagram recolored according to CNV. (C) Bar plot and pie chart depicting cell counts in all cell types (approximately 98.75% of aneuploid cells are positioned in epithelial cell clusters). (D) Tumor purity score for each cell type. (E) Cell counts in epithelial subcluster E1−E8 recolored according to CNV. (F) UMAP plots of 12,591 separated epithelial cells. (G) UMAP plots of epithelial subcluster E1−E8 recolored according to CNV. (H, I) Bar charts showing relative proportions of all the epithelial cells (H) and aneuploid cells (I) from each patient in epithelial subclusters E1−E8. (J, K) Pseudotime trajectories of malignant epithelial cells, colored by pseudotime (J) and subclusters (K). (L) Genes, such as BIRC5 and CDC20, were mainly expressed at the end of the trajectory. (M) Highly variable gene expression changes in the total trajectory. CHR: chromosome; CopyKAT: Copy number KAryotyping of Tumors.

Fig. 5

In view of the finding that tumor cells mainly originate from epithelial clusters, we extracted all cells from the epithelial cluster for analysis to explore the heterogeneity of TNBC tumors. All epithelia were re-classified into eight subclusters. In conjunction with CNV prediction, some healthy epithelial cells were included, although all samples were sourced from tumor tissue. E1 and E2 subclusters consisted predominantly of healthy epithelium, with the E3 cluster featuring a few healthy epithelial cells. These three clusters were positioned at the periphery of the overall cluster diagram, indicating significant disparities between healthy and tumor epithelial samples (Figs. 5E−G). When assessing the distribution of the eight epithelial subclusters, E1−E3 encompassed cells from nearly all sample sources. However, when focusing solely on tumor cell distribution, only the E8 cluster was common to all samples, while the remaining subclusters were distributed in a few samples. The observation that healthy epithelium exhibits high cellular similarity across TNBC samples whereas tumor epithelium displays remarkable variability further highlights the significant heterogeneity between tumor cells (Figs. 5H and I).

Notably, pseudotime trajectory analysis revealed that the E8 cluster shared among all samples was located at the terminal stage of differentiation. Intriguingly, while cells exhibited distinctive characteristics during the initial stages of tumor formation, a degree of similarity was observed with disease progression (Figs. 5J and K). According to Moran's Index, genes displaying the most significant changes in expression over time were classified into six major categories via hierarchical clustering. The fourth category, including BIRC5, CDC20, CDK1, and CENPF, showed elevated expression exclusively at the trajectory endpoint. This high expression profile showed significant overlap with the characteristic genes identified in the epithelial E8 cluster (Figs. 5L and M).

From a genetic perspective, all genes were clustered based on co-expression functions in epithelial cells, leading to classification into 10 functional modules (Fig. 6A). Module 10 was determined as the most distinctive module based on correlation analysis between functional modules and epithelial subclusters. In terms of gene function, Module 10 displayed minimal overlap and low correlation with the other modules. Notably, the eigengenes attributed to Module 10 exclusively exhibited marked expression levels in the E8 cluster. Additionally, the UMAP dimensionality reduction graph illustrated a significant correlation between eigengenes in Module 10 and the E8 epithelial subcluster (Figs. 6B−E).Fig. 6 Gene modules and functions. (A) Dendrogram of genes in triple-negative breast cancer (TNBC) epithelial cells divided into 10 modules. (B) Correlation between epithelial subclusters E1−E8 and gene function modules (Modules 1-10). (C) Correlations between different gene functional modules. (D) Uniform manifold approximation and projection (UMAP) plots of gene expression of Modules 1-10 at single-cell resolution. (E) Eigengenes in Module 10 showing significantly higher expression than the other modules. (F) Venn diagram showing intersections between genes at the end of pseudotime trajectory and eigengenes in Module 10. (G) Gene Ontology (GO) enrichment of 40 intersection genes. (H) Identification of TPX2, CDCA8, PLK1, UBE2S, RRM2, and TKI as risk genes and independent risk factors for TNBC. (I) Expression patterns of the six risk genes in TNBC and non-TNBC patients. (J) Poorer prognosis of the group with high and low expressions of the risk gene. ∗∗∗P < 0.001. CI: confidence interval; mRNA: messenger RNA.

Fig. 6

In view of the distinct cellular and genetic properties of the E8 epithelial subcluster, we further focused on its unique characteristics. To provide insights into the profile of E8, a comparative analysis was conducted between genes within the fourth category of pseudotime trajectory and those belonging to Module 10. Our data showed a substantial overlap between the groups, with a total of 40 shared characteristic genes (Fig. 6F). These crucial genes were enriched in key biological processes, including mitosis, organelle fisson, and chromosome segregation, which are highly related to cell proliferation (Fig. 6G).

Using the above 40 genes in univariate analysis of the METABRIC database, TPX2, CDCA8, PLK1, UBE2S, RRM2, and TK1 were identified as independent risk factors affecting overall survival. Although these genes were not exclusively expressed in TNBC, their expression levels were notably increased compared to other subtypes. Patients were further stratified into high and low expression groups based on gene levels. Examination of the 10-year survival rates revealed significantly more favorable prognosis in the low expression group (Figs. 6H−J).

3.6 Screening of potential therapeutic drugs for TNBC

Disease development is not only related to tumor cells but also influenced by the TME. Accordingly, we explored other cellular components in TME that could impact the E8 subcluster. Examination of the ligand-receptor interaction network between E8 and other cell types revealed relatively frequent interactions with fibroblasts, myeloid cells, and endothelial cells, but fewer interactions with lymphocytes. Given that TNBC contains limited endothelial cells, we further focused on E8-related fibroblasts and myeloid cells (Figs. 7A and B).Fig. 7 Identification of malignant associated fibroblasts and myeloid cells as well as potential therapeutic drugs. (A, B) Circle plot (A) and bar plot (B) showing ligand-receptor interactions between the epithelial subcluster E8 and the other cell types. (C, D) Uniform manifold approximation and projection (UMAP) plot (C) and pseudotime trajectory (D) of triple-negative breast cancer (TNBC) fibroblasts. (E) Normal score and epithelial-mesenchymal transition (EMT) score of fibroblast subclusters Fib_1 to Fib_4. (F) Ligand-receptor interactions between E8 and fibroblast subclusters. (G, H) UMAP plot (G) and pseudotime trajectory (H) in TNBC myeloid cells. 6 and 7: branching points in pseudotime trajectory. (I) Violin plots of expression of marker genes in monocytes (FCN1), dendritic cells (CD1E), M1 macrophages (CD80), M2 macrophages (MRC1 and CCL18), T cell receptor (TCR)+ macrophages (CD3D), proliferating cells (TYMS), and mast cells (TPSAB1). (J) Ligand-receptor interactions between E8 and myeloid cell subclusters. (K, L) Circle plot showing costimulatory interactions between E8 and fibroblast subclusters (K) and E8 and myeloid cell subclusters (L). (M) Protein-protein interaction network constructed from marker genes and ligand-receptor genes of three subclusters. (N) Potentially effective drugs for TNBC determined based on therapeutic scores. ∗∗∗P < 0.001. FDR: false discovery rate.

Fig. 7

Fibroblasts were divided into four subclusters, designated Fib_1 to Fib_4. Subcluster Fib_4 mainly exhibited characteristics of myofibroblasts [53] and was markedly different on the dimensionality reduction plot, showing no link to the other three subclusters on a pseudotime trajectory (Figs. 7C−E). To explore the direction of cell differentiation, we utilized a transcriptional profile of fibroblasts from healthy breast tissue samples obtained from GSE180878 as the control [54]. Fib_1 displayed the highest normal score, indicating that the subcluster is inherently similar to healthy fibroblasts and could serve as the origin for differentiation. Correspondingly, Fib_3 was located at the end of the differentiation trajectory.

KEGG enrichment analysis revealed significant enrichment of subcluster Fib_3 in the epithelial-mesenchymal transition (EMT) pathway. Fib_3 displayed a higher score in hallmark_EMT gene set scoring, leading to definition of this subcluster as EMT fibroblasts (Fig. 7E). In the context of ligand-receptor interactions, EMT fibroblasts displayed the most pronounced interactions with epithelial subcluster E8. Notably, LGALS9 played a crucial role in the TME [55], functioning as a ligand that could regulate multiple genes in E8 (Fig. 7F).

We identified a total of 7 distinct cell subclusters within myeloid cells of TNBC, which were annotated as monocytes (FCN1), dendritic cells (CD1E), M1 macrophages (CD80), M2 macrophages (MRC1 and CCL18), T cell receptor (TCR)+ macrophages (CD3D), proliferating cells (TYMS), and mast cells (TPSAB1) based on specific cell markers (Figs. 7G−I). Mast cells were not associated with other cells in the pseudotime trajectory. However, the developmental process from monocytes to dendritic cells and macrophages was clearly reflected [56]. Macrophages were positioned at the end of the trajectory, signifying a terminal differentiation state of myeloid cells. While M1 and M2 macrophages showed some tendency of mutual transformation, M2 cells were predominant and occupied the terminal position. Examination of the interactions between epithelial E8 and all myeloid cells revealed that M2 macrophages had the highest number of ligands interacting with E8. For example, NRP2, a gene expressed in M2 macrophages and fibroblasts, is reported to be involved in promoting efferocytosis and facilitating tumor growth (Fig. 7J). In terms of the costimulatory factors, significant associations of E8 with EMT fibroblasts and M2 were confirmed (Figs. 7K and L).

We further constructed a protein-protein interaction network based on characteristic genes of epithelial E8, EMT fibroblasts, and M2 macrophages, along with those involved in ligand-receptor interactions (Table S2). A network diagram of highly correlated genes with connectivity degrees ≥ 30 was generated (Fig. 7M). Within these three crucial cell subclusters, we focused on genes with connectivity degrees ≥ 1 and subsequently applied the ASGARD algorithm to assess the effectiveness of all drugs in the LINCS L1000 database (Table S3). Through this analysis, several potential therapeutic drugs were identified, including trametinib (a mitogen-activated protein (MEK)1/2 inhibitor), pelitinib (an epidermal growth factor receptor (EGFR) inhibitor), geldanamycin (an antibacterial drug), and roscovitine (a cyclin dependent kinase (CDK) inhibitor) (Fig. 7N). While these drugs have been clearly shown to inhibit specific genes or functions that lead to tumor suppression, they also possess additional anticancer mechanisms. All the above agents could effectively inhibit tumor progression by suppressing the EMT process in tumor epithelial cells and impeding M2 polarization of macrophages.

3.7 Anti-tumor activity and pathways of pelitinib in vivo

Pelitinib, among the drugs we screened, exhibited a high therapeutic score and has been documented favorable therapeutic results in various cancers [[57], [58], [59]], but few studies have been conducted in TNBC. Therefore, we proceeded to investigate the anti-tumor efficacy of pelitinib in TNBC using in vivo experiments.

In this study, pelitinib demonstrated significant inhibitory effects on tumor volume and tumor weight (P < 0.001), showing a dose-dependent relationship (Figs. 8A−C). Histopathological analysis revealed that the control group exhibited normal tumor tissue structure without degeneration or inflammation. However, increasing doses of pelitinib led to structural abnormalities in the tumor, including disorganized arrangements, nuclear pleomorphism, and infiltration of deeply stained inflammatory cells. Immunohistochemical analysis showed no change in the overall abundance of macrophages after pelitinib treatment compared to the control group. Nonetheless, there was a notable decrease in M2-polarized macrophages with higher drug concentrations, while M1-polarized macrophages increased (Fig. 8D). Western blot results indicated upregulation of E-cadherin and ZO-1, indicating epithelial characteristics after pelitinib treatment, along with downregulation of β-catenin and vimentin, representing mesenchymal traits. Additionally, the key transcription factor snail exhibited a significant decrease after pelitinib administration (Figs. 8E and F). The findings are consistent with previous bioinformatics data, supporting the role of pelitinib in inhibiting macrophage M2 polarization and impeding the tumor EMT process, thus confirming the robustness of our conclusions.Fig. 8 Pelitinib demonstrated anti-triple-negative breast cancer (TNBC) effects in animal experiments. (A) Changes in tumor volume in each group of mice. (B) Mouse tumor weight at the time of sacrifice. (C) Gross photograph of mouse tumors in different groups. (D) Hematoxylin and eosin (H&E) staining shows more structural abnormalities in tumors as a dose of pelitinib rises. Immunohistochemical analysis of F4/80, CD206, and CD86 of tumor sections in the control group and treatment with different doses of pelitinib. (E) Expression of epithelial-mesenchymal transition (EMT)-related proteins E-cadherin, zonula occludens-1 (ZO-1), β-catenin, vimentin, and snail, a key transcription factor in the EMT process, in different treatment groups by Western blot. (F) Western blot quantitation of proteins about EMT and phosphoinositide-3-kinase (PI3K)-protein kinase B (Akt) passway normalized relative to β-actin or Akt. (G) Expression of PI3K, Akt, and Akt in different treatment groups. ∗P < 0.05, ∗∗P < 0.01, and ∗∗∗P < 0.001.

Fig. 8

To further explore the relevant pathways to the inhibition of TNBC by pelitinib, we searched the LINCS L1000 database for genes upregulated in TNBC cells after pelitinib treatment and compared them with previously obtained genes with high connectivity. We found 487 common TNBC risk genes inhibited by pelitinib. KEGG enrichment analysis identified 10 signaling pathways inhibited by pelitinib (Fig. S3). We further investigated the phosphoinositide-3-kinase (PI3K)-Akt pathway, which showed the highest and most significant gene suppression by pelitinib (Fig. S4) [20]. Western blot analysis showed that compared to the control group, the expression of p-Akt in the low, medium, and high-dose pelitinib groups was significantly decreased in a dose-dependent manner. There were no statistically significant differences in the expression of Akt and PI3K among the groups (Figs. 8F and G).

The safety of pelitinib was also assessed through blood tests and histology. Pelitinib showed no significant myelosuppression or hematologic toxicity on complete blood counts. Similarly, plasma aspartate aminotransferase (AST) and alanine aminotransferase (ALT) levels, which reflect hepatic function, and uric acid (UA), creatinine (Cr), and blood urea nitrogen (BUN) levels, which reflect renal function, also remained stable after treatment, demonstrating the hepatic and renal safety of pelitinib. H&E stained sections of hearts, livers, spleens, lungs, and kidneys showed that the morphology of each organ was free of significant lesions at different dose levels, indicating the positive biological safety of pelitinib (Figs. S5 and S6).

4 Discussion

In this study, we comprehensively explored the distinct characteristics of different breast cancer subtypes, with the aim of providing insights into the complex interactions within tumor cells and the TME in cancer patients. This investigation revealed remarkable heterogeneity between and within breast cancer subtypes. In particular, the significant proliferative activity of TNBC was highlighted. Additionally, our results showed that genes affected during transitioning from normal to cancer cells reflect tumor grading and prognosis in TNBC. The TNBC subtype has unique complexity and presents a significant challenge in terms of both clinical treatment and basic research [60], highlighting the necessity to develop novel innovative therapeutic approaches.

Wu et al. [61] developed a single-cell method (scSubtype) to reveal neoplastic cell heterogeneity. TNBC-related biological modules are enriched in cell cycle, interferon response, antigen presentation, and EMT functions. Besides, the authors revealed nine tumor clusters with similar ‘tumor ecotypes’, in which Ecotype-3 points towards TNBC and exhibits the worst prognosis. Pal et al. [62] provided an overview of the healthy and diseased state of the breast from multiple perspectives, elucidating changes in the breast microenvironment and accompanying hormonal status. Single-cell profiling of 34 treatment-naive primary tumors revealed comparable diversity among cancer cells and immune landscapes. In our study, instead of limiting to a specific cell type, we considered multiple cell clusters associated with TNBC progression in the TME as a whole, and searched for effective therapeutic regimens that could act on these cells simultaneously.

In-depth examination of TNBC at the single-cell level revealed a specific subcluster of tumor epithelial cells consistently present in all patients, which was strongly associated with increased tumor cell proliferation. The majority of genes within this group exerted a pro-tumor effect, with TPX2, CDCA8, PLK1, UBE2S, RRM2, and TK1 emerging as independent risk factors underlying poorer patient prognosis.

Oncotype DX breast recurrence score [63] is a 21-gene prognostic and predictive assay that includes five proliferation-related genes. Among these genes, MKI67, CCNB1, and BIRC5 belong to the risk factors identified. Other genes not incorporated in the 21-gene score, such as RRM2, promote cell invasion or migration through the PI3K/Akt signaling pathway, regulate the EMT process, and play a pivotal role in the transformation of healthy to malignant epithelium. These genes collectively represent a signature of TNBC characterized by high proliferative activity.

Extensive communication between tumor cells and their surrounding components regulates proliferative capacity. The most relevant fibroblast clusters were highly enriched in the EMT pathway, a biological process that significantly enhances migration, causing cells to lose their epithelial characteristics and transform into a mesenchymal phenotype [64,65]. Remodeling of the microtubule skeleton triggers diffusive behaviors of integrins, promoting escape of malignant cells from the primary tumor site, acceleration of tumor cell migration, and acquisition of stronger invasion ability [66].

M2 macrophages also display complex connections with highly proliferative tumor cells [67]. Macrophages are crucial immune cells in TME, categorized into classically activated M1 and alternately activated M2 types with distinct functions [68]. M1 macrophages primarily function in T helper 1 (Th1) cell recruitment, pathogen resistance, and tumor killing primarily through natural and adaptive immune responses [69], while M2 macrophages are involved in tumor progression through various pathways. M2 macrophages typically exhibit anti-inflammatory and immunomodulatory properties that inhibit the immune defense system against tumors, contributing to tumor evasion of immune detection and clearance [70]. In addition, M2 macrophages secrete growth factors, such as vascular endothelial growth factor (VEGF) and epidermal growth factor (EGF), which promote the formation of new blood vessels and supply oxygen and nutrients required by tumors [71,72]. M2 macrophages further secrete matrix metallopeptidases (MMPs) [73], which help to degrade the extracellular matrix (ECM) and support tumor cell spread.

Non-specific chemotherapy remains the standard treatment for the non-surgical phase of TNBC [7]. To target the three crucial cell subclusters in TNBC, trametinib, AS-605240, pelitinib, geldanamycin, roscovitine, and other therapeutic drugs were screened. Within the same category of drugs, variations in modification groups or spatial arrangements could induce markedly different effects. Trametinib is a MEK1/2 inhibitor with anti-tumor activity approved by the U.S. Food and Drug Administration for the treatment of melanoma, non-small cell lung cancer, thyroid cancer, low-grade glioma, and other solid tumors [74,75]. Mammary gland-related clinical trials have been conducted on trametinib, reflecting its research value and potential [76]. AS-605240, a selective inhibitor of PI3Kγ, is released by mitochondria in injured organelles and contained in autolysosomes, acting on PI3Kγ/Akt/mammalian target of rapamycin (mTOR)/unc-51 like autophagy activating kinase 1 (Ulk1) signaling, which provides maladaptive feedback inhibition of autophagy. AS-605240 interferes with tumor growth through selective inhibition of PI3Kγ [77,78]. Geldanamycin shows strong adenosine triphosphatase (ATPase) inhibitory activity and targets heat shock protein (Hsp90), leading to caspase-dependent apoptosis and necroptosis [79]. Roscovitine, which primarily acts on CDK2/5, was one of the first CDK inhibitors to enter clinical trials [80,81]. Following its poor performance in previous clinical trials, roscovitine derivatives N&N1 and LGR1406 showed improved anti-tumor properties and induced apoptosis in a diverse set of cell lines [82]. In addition, more than 20 other drugs possessing therapeutic potential were included in our calculations, emerging as potential novel options for future treatment strategies.

We focus on pelitinib, which has been identified as a potent irreversible EGFR inhibitor, and has been studied for clinical use in other tumors such as lung cancer, liver cancer, and cervical cancer. In addition, pelitinib may exert its tumor-suppressive effect through other mechanisms or pathways that have not been confirmed. Studies have shown that pelitinib can upregulate the ABCB1/ABCG2 efflux transporter, thereby eradicating cancer stem cells in lung cancer [57]. Pelitinib can also induce the degradation of the transcription factor Twist1 during the EMT process by inhibiting the MEK kinase (MAPK)/Akt signaling pathway, thereby inhibiting the migration and invasion of hepatocellular carcinoma cells [58]. Furthermore, cervical cancer patients with PIK3CA mutations have higher sensitivity to pelitinib [59].

According to our analysis in this study, TNBC tumor weight and volume were significantly reduced after pelitinib treatment on 4T1 tumor-bearing model. IHC and Western blot results showed that it had a significant inhibitory effect on the EMT process and the polarization of M2 macrophages in tumors. In addition, we demonstrated that the anti-tumor effect of pelitinib is related to the inhibition of the PI3K-Akt signaling pathway. Finally, the high biological safety of pelitinib was also confirmed in the assessment of adverse events. All of these results demonstrate the efficacy and safety of pelitinib in the treatment of TNBC and provide a basic medical reference for the development of pelitinib as a new option for the treatment of TNBC.

The analytical process of this study could provide some help to guide individualized medication in clinical practice. This approach provides a broader idea of drug repurposing, and other anti-tumor mechanisms of existing drugs can be explored through a combination of bioinformatics analysis and experimental methods. Based on our analysis and following experiments, the drugs selected did indeed exhibit significant inhibitory effects on TNBC, verifying the reliability of the drugs screened in this study. In clinical practice, a patient's tumor tissue could be collected for scRNA-seq during surgery, and the transcriptional profile could be analyzed to identify potential effective drugs for that particular patient. When clinical patients have tried multiple lines of therapy and still have not achieved satisfactory efficacy, the analysis method of this study could be used to provide a wider range of potential drug choices and bring new hope to patients with late-stage disease. In addition, the first-line treatment of TNBC is currently dominated by chemotherapy, and in the high-risk patient population, the combination of the drugs screened in this study with chemotherapeutics may also be used in the future as an individualized regimen for combining first-line drugs.

5 Conclusions

In conclusion, this article analyzed the cellular characteristics of various molecular subtypes of breast cancer, emphasizing the unique features of TNBC from a single-cell perspective. Furthermore, our study revealed a shared proliferating malignant epithelium cluster among all TNBC patients, which is closely associated with EMT fibroblasts and M2 macrophages in TME. We described the genetic characteristics of these three clusters and identified potential therapeutic drugs accordingly, thereby expanding the options for TNBC treatment.

CRediT authorship contribution statement

Weilun Cheng: Writing – review & editing, Writing – original draft, Visualization, Software, Formal analysis, Conceptualization. Wanqi Mi: Visualization, Software, Methodology. Shiyuan Wang: Writing – review & editing, Validation, Investigation, Data curation. Xinran Wang: Validation, Formal analysis. Hui Jiang: Writing – review & editing. Jing Chen: Resources. Kaiyue Yang: Investigation. Wenqi Jiang: Data curation. Jun Ye: Writing – review & editing, Writing – original draft. Baoliang Guo: Supervision, Funding acquisition. Yunpeng Zhang: Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that there are no conflicts of interest.

Appendix A Supplementary data

The following are the Supplementary data to this article:Multimedia component 1

Multimedia component 1

Multimedia component 2

Multimedia component 2

Acknowledgments

This research was funded by the 10.13039/501100001809 National Natural Science Foundation of China (Grant Nos: 62172131 and 81872135 ) and the Outstanding Youth Foundation of Heilongjiang Province of China (Grant No.: YQ2021C026 ).

Peer review under responsibility of Xi'an Jiaotong University.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.jpha.2024.100975.
==== Refs
References

1 Sung H. Ferlay J. Siegel R.L. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries CA Cancer J. Clin. 71 2021 209 249 33538338
2 Perou C.M. Sørlie T. Eisen M.B. Molecular portraits of human breast tumours Nature 406 2000 747 752 10963602
3 Goldhirsch A. Wood W.C. Coates A.S. Strategies for subtypes – Dealing with the diversity of breast cancer: Highlights of the St. Gallen International Expert Consensus on the Primary Therapy of Early Breast Cancer 2011 Ann. Oncol. 22 2011 1736 1747 21709140
4 Howard F.M. Olopade O.I. Epidemiology of triple-negative breast cancer: A review Cancer J. 27 2021 8 16 33475288
5 Harbeck N. Gnant M. Breast cancer Lancet 389 2017 1134 1150 27865536
6 Gradishar W.J. Moran M.S. Abraham J. NCCN guidelines® insights: Breast cancer, version 4.2023 J. Natl. Compr. Canc. Netw. 21 2023 594 608 37308117
7 Li Y. Zhang H. Merkher Y. Recent advances in therapeutic strategies for triple-negative breast cancer J. Hematol. Oncol. 15 2022 121
8 Ding S. Chen X. Shen K. Single-cell RNA sequencing in breast cancer: Understanding tumor heterogeneity and paving roads to individualized therapy Cancer Commun. (Lond.) 40 2020 329 344 32654419
9 Sklavenitis-Pistofidis R. Getz G. Ghobrial I. Single-cell RNA sequencing: One step closer to the clinic Nat. Med. 27 2021 375 376 33664491
10 Zhang Y. Chen H. Mo H. Single-cell analyses reveal key immune cell subsets associated with response to PD-L1 blockade in triple-negative breast cancer Cancer Cell 39 2021 1578 1593.e8 34653365
11 Liu T. Liu C. Yan M. Single cell profiling of primary and paired metastatic lymph node tumors in breast cancer patients Nat. Commun. 13 2022 6823
12 Karaayvaz M. Cristea S. Gillespie S.M. Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq Nat. Commun. 9 2018 3588
13 Sebastian A. Hum N.R. Martin K.A. Single-cell transcriptomic analysis of tumor-derived fibroblasts and normal tissue-resident fibroblasts reveals fibroblast heterogeneity in breast cancer Cancers 12 2020 1307
14 Curtis C. Shah S.P. Chin S.-F. The genomic and transcriptomic architecture of 2, 000 breast tumours reveals novel subgroups Nature 486 2012 346 352 22522925
15 Cerami E. Gao J. Dogrusoz U. The cBio cancer genomics portal: An open platform for exploring multidimensional cancer genomics data Cancer Discov. 2 2012 401 404 22588877
16 Hao Y. Stuart T. Kowalski M.H. Dictionary learning for integrative, multimodal and scalable single-cell analysis Nat. Biotechnol. 42 2024 293 304 37231261
17 Korsunsky I. Millard N. Fan J. Fast, sensitive and accurate integration of single-cell data with Harmony Nat. Meth. 16 2019 1289 1296
18 Subramanian A. Tamayo P. Mootha V.K. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles Proc. Natl. Acad. Sci. U. S. A. 102 2005 15545 15550 16199517
19 Wu T. Hu E. Xu S. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data Innovation (Camb) 2 2021 100141
20 Kanehisa M. Furumichi M. Sato Y. KEGG for taxonomy-based analysis of pathways and genomes Nucleic Acids Res. 51 2023 D587 D592 36300620
21 Patel A.P. Tirosh I. Trombetta J.J. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma Science 344 2014 1396 1401 24925914
22 Gao R. Bai S. Henderson Y.C. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes Nat. Biotechnol. 39 2021 599 608 33462507
23 Cao J. Spielmann M. Qiu X. The single-cell transcriptional landscape of mammalian organogenesis Nature 566 2019 496 502 30787437
24 Yoshihara K. Shahmoradgoli M. Martínez E. Inferring tumour purity and stromal and immune cell admixture from expression data Nat. Commun. 4 2013 2612
25 Morabito S. Reese F. Rahimzadeh N. hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data Cell Rep. Methods 3 2023 100498
26 Garcia-Alonso L. Lorenzi V. Mazzeo C.I. Single-cell roadmap of human gonadal development Nature 607 2022 540 547 35794482
27 Shannon P. Markiel A. Ozier O. Cytoscape: A software environment for integrated models of biomolecular interaction networks Genome Res. 13 2003 2498 2504 14597658
28 He B. Xiao Y. Liang H. ASGARD is A single-cell guided pipeline to aid repurposing of drugs Nat. Commun. 14 2023 993
29 Detre S. Saclani Jotti G. Dowsett M. A “quickscore” method for immunohistochemical semiquantitation: Validation for oestrogen receptor in breast carcinomas J. Clin. Pathol. 48 1995 876 878 7490328
30 Hu C. Li T. Xu Y. CellMarker 2.0: An updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data Nucleic Acids Res. 51 2023 D870 D876 36300619
31 Huang W.C. Yen J.H. Sung Y.W. Novel function of THEMIS2 in the enhancement of cancer stemness and chemoresistance by releasing PTP1B from MET Oncogene 41 2022 997 1010 34974522
32 LeBlanc V.G. Trinh D.L. Aslanpour S. Single-cell landscapes of primary glioblastomas and matched explants and cell lines show variable retention of inter- and intratumor heterogeneity Cancer Cell 40 2022 379 392.e9 35303420
33 Zheng S. Zou Y. Tang Y. Landscape of cancer-associated fibroblasts identifies the secreted biglycan as a protumor and immunosuppressive factor in triple-negative breast cancer Oncoimmunology 11 2022 2020984
34 Hu L. Su L. Cheng H. Single-cell RNA sequencing reveals the cellular origin and evolution of breast cancer in BRCA1 mutation carriers Cancer Res 81 2021 2600 2611 33727227
35 Yan X. Xie Y. Yang F. Comprehensive description of the current breast cancer microenvironment advancements via single-cell analysis J. Exp. Clin. Cancer Res. 40 2021 142
36 Pecoraro M. Marzocco S. Franceschelli S. Trastuzumab and doxorubicin sequential administration increases oxidative stress and phosphorylation of connexin 43 on Ser368 Int. J. Mol. Sci. 23 2022 6375
37 Heublein S. Mayr D. Meindl A. Vitamin D receptor, Retinoid X receptor and peroxisome proliferator-activated receptor γ are overexpressed in BRCA1 mutated breast cancer and predict prognosis J. Exp. Clin. Cancer Res. 36 2017 57
38 Anurag M. Jaehnig E.J. Krug K. Proteogenomic markers of chemotherapy resistance and response in triple-negative breast cancer Cancer Discov. 12 2022 2586 2605 36001024
39 Vegunta S. Kling J.M. Kapoor E. Androgen therapy in women J. Womens Health (Larchmt) 29 2020 57 64 31687883
40 Edwards D.N. Ngwa V.M. Raybuck A.L. Selective glutamine metabolism inhibition in tumor cells improves antitumor T lymphocyte activity in triple-negative breast cancer J. Clin. Invest. 131 2021 e140100
41 Deepak K.G.K. Vempati R. Nagaraju G.P. Tumor microenvironment: Challenges and opportunities in targeting metastasis of triple negative breast cancer Pharmacol. Res. 153 2020 104683
42 Whittle J.R. Vaillant F. Surgenor E. Dual targeting of CDK4/6 and BCL2 pathways augments tumor response in estrogen receptor-positive breast cancer Clin. Cancer Res. 26 2020 4120 4134 32245900
43 Periasamy V.S. Riyasdeen A. Rajendiran V. Induction of redox-mediated cell death in ER-positive and ER-negative breast cancer cells by a copper(II)-phenolate complex: An in vitro and in silico study Molecules 25 2020 4504
44 Taurin S. Alkhalifa H. Breast cancers, mammary stem cells, and cancer stem cells, characteristics, and hypotheses Neoplasia 22 2020 663 678 33142233
45 Vishnubalaji R. Alajez N.M. Epigenetic regulation of triple negative breast cancer (TNBC) by TGF-β signaling Sci. Rep. 11 2021 15410
46 Böcker W. WHO classification of breast tumors and tumors of the female genital organs: Pathology and genetics Verh. Dtsch Ges. Pathol. 86 2002 116 119 12647359
47 Liu N. Wang X. Zhu Z. Selected ideal natural ligand against TNBC by inhibiting CDC20, using bioinformatics and molecular biology Aging 13 2021 23702 23725 34686627
48 Liu Y. Teng L. Fu S. Highly heterogeneous-related genes of triple-negative breast cancer: Potential diagnostic and prognostic biomarkers BMC Cancer 21 2021 644
49 Li J. Gao X. Zhang Z. CircCD44 plays oncogenic roles in triple-negative breast cancer by modulating the miR-502-5p/KRAS and IGF2BP2/Myc axes Mol. Cancer 20 2021 138
50 Li W. Tanikawa T. Kryczek I. Aerobic glycolysis controls myeloid-derived suppressor cells and tumor immunity via a specific CEBPB isoform in triple-negative breast cancer Cell Metab. 28 2018 87 103.e6 29805099
51 Evans K.W. Yuca E. Scott S.S. Oxidative phosphorylation is a metabolic vulnerability in chemotherapy-resistant triple-negative breast cancer Cancer Res. 81 2021 5572 5581 34518211
52 Wu Q. Nie D.Y. Ba-Alawi W. PRMT inhibition induces a viral mimicry response in triple-negative breast cancer Nat. Chem. Biol. 18 2022 821 830 35578032
53 O’Flanagan C.H. Campbell K.R. Zhang A.W. Dissociation of solid tumor tissues with cold active protease for single-cell RNA-seq minimizes conserved collagenase-associated stress responses Genome Biol. 20 2019 210
54 Gray G.K. Li C.M. Rosenbluth J.M. A human breast atlas integrating single-cell proteomics and transcriptomics Dev. Cell 57 2022 1400 1420.e7 35617956
55 Yoon H.K. Kim T.H. Park S. Effect of anthracycline and taxane on the expression of programmed cell death ligand-1 and galectin-9 in triple-negative breast cancer Pathol. Res. Pract. 214 2018 1626 1631 30139555
56 Saeed S. Quintin J. Kerstens H.H.D. Epigenetic programming of monocyte-to-macrophage differentiation and trained innate immunity Science 345 2014 1251086
57 To K.K. Poon D.C. Wei Y. Pelitinib (EKB-569) targets the up-regulation of ABCB1 and ABCG2 induced by hyperthermia to eradicate lung cancer Br. J. Pharmacol. 172 2015 4089 4106 25988710
58 Lee S. Kang E. Lee U. Role of pelitinib in the regulation of migration and invasion of hepatocellular carcinoma cells via inhibition of Twist1 BMC Cancer 23 2023 703
59 Lv X. Jia Y. Li J. The construction of a prognostic model of cervical cancer based on four immune-related LncRNAs and an exploration of the correlations between the model and oxidative stress Front. Pharmacol. 14 2023 1234181
60 Derakhshan F. Reis-Filho J.S. Pathogenesis of triple-negative breast cancer Annu. Rev. Pathol. 17 2022 181 204 35073169
61 Wu S.Z. Al-Eryani G. Roden D.L. A single-cell and spatially resolved atlas of human breast cancers Nat. Genet. 53 2021 1334 1347 34493872
62 Pal B. Chen Y. Vaillant F. A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast EMBO J. 40 2021 e107333
63 Sparano J.A. Gray R.J. Makower D.F. Prospective validation of a 21-gene expression assay in breast cancer N. Engl. J. Med. 373 2015 2005 2014 26412349
64 Dongre A. Weinberg R.A. New insights into the mechanisms of epithelial-mesenchymal transition and implications for cancer Nat. Rev. Mol. Cell Biol. 20 2019 69 84 30459476
65 Bracken C.P. Goodall G.J. The many regulators of epithelial-mesenchymal transition Nat. Rev. Mol. Cell Biol. 23 2022 89 90 34887545
66 Yuan J. Zhang Y. Liu Y. Diffusion behaviors of integrins in single cells altered by epithelial to mesenchymal transition Small 18 2022 e2106498
67 Chen Y. Zhang S. Wang Q. Tumor-recruited M2 macrophages promote gastric and breast cancer metastasis via M2 macrophage-secreted CHI3L1 protein J. Hematol. Oncol. 10 2017 36
68 Mantovani A. Marchesi F. Malesci A. Tumour-associated macrophages as treatment targets in oncology Nat. Rev. Clin. Oncol. 14 2017 399 416 28117416
69 Zhao X. Di Q. Liu H. MEF2C promotes M1 macrophage polarization and Th1 responses Cell. Mol. Immunol. 19 2022 540 553 35194174
70 Zhu X. Liang R. Lan T. Tumor-associated macrophage-specific CD155 contributes to M2-phenotype transition, immunosuppression, and tumor progression in colorectal cancer J. Immunother. Cancer 10 2022 e004219
71 Wheeler K.C. Jena M.K. Pradhan B.S. VEGF may contribute to macrophage recruitment and M2 polarization in the decidua PLoS One 13 2018 e0191040
72 Chen X. Gao A. Zhang F. ILT4 inhibition prevents TAM- and dysfunctional T cell-mediated immunosuppression and enhances the efficacy of anti-PD-L1 therapy in NSCLC with EGFR activation Theranostics 11 2021 3392 3416 33537094
73 Maybee D.V. Ink N.L. Ali M.A.M. Novel roles of MT1-MMP and MMP-2: Beyond the extracellular milieu Int. J. Mol. Sci. 23 2022 9513
74 Flaherty K.T. Robert C. Hersey P. Improved survival with MEK inhibition in BRAF-mutated melanoma N Engl J. Med. 367 2012 107 114 22663011
75 Maio M. Carlino M.S. Joshua A.M. KEYNOTE-022: Pembrolizumab with trametinib in patients with BRAF wild-type melanoma or advanced solid tumours irrespective of BRAF mutation Eur. J. Cancer 160 2022 1 11 34801354
76 Seo T. Noguchi E. Yoshida M. Response to dabrafenib and trametinib of a patient with metaplastic breast carcinoma harboring a BRAF V600E mutation, Case Rep. Oncol Med. 2020 2020 2518383
77 Schmid M.C. Avraamides C.J. Dippold H.C. Receptor tyrosine kinases and TLR/IL1Rs unexpectedly activate myeloid cell PI3Kγ, a single convergent point promoting tumor inflammation and progression Cancer Cell 19 2011 715 727 21665146
78 Li M. Sala V. De Santis M.C. Phosphoinositide 3-kinase gamma inhibition protects from anthracycline cardiotoxicity and reduces tumor growth Circulation 138 2018 696 711 29348263
79 Zhang Z. Li H. Zhou C. Non-benzoquinone geldanamycin analogs trigger various forms of death in human breast cancer cells J. Exp. Clin. Cancer Res. 35 2016 149
80 Zhang M. Zhang L. Hei R. CDK inhibitors in cancer therapy, an overview of recent development Am. J. Cancer Res 11 2021 1913 1935 34094661
81 Cicenas J. Kalyan K. Sorokinas A. Roscovitine in cancer and other diseases Ann. Transl. Med. 3 2015 135
82 Houzé S. Hoang N.T. Lozach O. Several human cyclin-dependent kinase inhibitors, structurally related to roscovitine, are new anti-malarial agents Molecules 19 2014 15237 15257 25251193
