Single-Cell Transcriptome Analysis Revealed Heterogeneity and Identified Novel Therapeutic Targets for Breast Cancer Subtypes

Abstract

Breast cancer (BC) is a heterogeneous disease, which is primarily classified according to hormone receptors and HER2 expression. Despite the many advances in BC diagnosis and management, the identification of novel actionable therapeutic targets expressed by cancerous cells has always been a daunting task due to the large heterogeneity of the disease and the presence of non-cancerous cells (i.e., immune cells and stromal cells) within the tumor microenvironment. In the current study, we employed computational algorithms to decipher the cellular composition of estrogen receptor-positive (ER⁺), HER2⁺, ER⁺HER2⁺, and triple-negative BC (TNBC) subtypes from a total of 49,899 single cells' publicly available transcriptomic data derived from 26 BC patients. Restricting the analysis to EPCAM⁺Lin^- tumor epithelial cells, we identified the enriched gene sets in each BC molecular subtype. Integration of single-cell transcriptomic with CRISPR-Cas9 functional screen data identified 13 potential therapeutic targets for ER⁺, 44 potential therapeutic targets for HER2⁺, and 29 potential therapeutic targets for TNBC. Interestingly, several of the identified therapeutic targets outperformed the current standard of care for each BC subtype. Given the aggressive nature and lack of targeted therapies for TNBC, elevated expression of ENO1, FDPS, CCT6A, TUBB2A, and PGK1 predicted worse relapse-free survival (RFS) in basal BC (n = 442), while elevated expression of ENO1, FDPS, CCT6A, and PGK1 was observed in the most aggressive BLIS TNBC subtype. Mechanistically, targeted depletion of ENO1 and FDPS halted TNBC cell proliferation, colony formation, and organoid tumor growth under 3-dimensional conditions and increased cell death, suggesting their potential use as novel therapeutic targets for TNBC. Differential expression and gene set enrichment analysis in TNBC revealed enrichment in the cycle and mitosis functional categories in FDPS^high, while ENO1^high was associated with numerous functional categories, including cell cycle, glycolysis, and ATP metabolic processes. Taken together, our data are the first to unravel the unique gene signatures and to identify novel dependencies and therapeutic vulnerabilities for each BC molecular subtype, thus setting the foundation for the future development of more effective targeted therapies for BC.

Keywords: TNBC; breast cancer; molecular subtypes; single cell analysis; therapeutic targets.

Figure 1

BC heterogeneity delineated by single-cell transcriptome analysis. Single-cell analysis of a total of 16,350 single cells derived from 12 ER⁺ BC (a), 7824 single cells derived from 3 HER2⁺ BC (b), 11,487 cells derived from 2 ER⁺HER2⁺ BC (c), and 14,238 single cells derived from 9 TNBC (d) employing AltAnalyze algorithm depicted as heat map. The lower panels represent UMAP dimensionality reduction visualization of the identified cell clusters. The text on the left of each heatmap indicates enriched cell-type markers from the default gene-set enrichment analysis and corresponding “Z” score p value.

Figure 2

Identification of gene biomarkers enriched in each BC molecular subtype. (a) Bar chart depicting the number of differentially expressed genes (DEGs) in ER⁺ vs. HER2⁺, ER⁺ vs. ER⁺HER2⁺, HER2⁺ vs. ER⁺HER2⁺, TNBC vs. ER⁺, TNBC vs. HER2⁺, and TNBC vs. ER⁺HER2⁺. (b) PCA illustrating the segregation of each BC molecular subtype based on the identified gene markers. (c) Heatmap depicting the enriched gene markers associated with TNBC, ER⁺, HER2⁺, and ER⁺HER2⁺ employing the MarkerFinder algorithm. The text on the left of each heatmap indicates enriched Go functional categories in each BC subtype and the corresponding enrichment p values. (d) Expression of selected genes enriched in TNBC (upper panel) or ER⁺ (lower panel) based on RNA–Seq data. **** p < 0.00005.

Figure 3

Discriminative analysis based on the identified gene markers. Anova analysis for the identified gene markers enriched in TNBC (a), ER⁺ (b), HER2⁺ (c), and ER⁺HER2⁺ (d). Sum score from the identified gene markers for each BC subtype were used for Anova analysis. (e) OPLS-DA score plot for the different BC molecular subtypes based on the identified gene markers. (f) ROC analysis of the identified gene markers comparing the four BC molecular subtypes.

Figure 4

Identification of novel therapeutic targets for each BC molecular subtype employing CRISPR-Cas9 perturbational effects. Violin plots illustrating the perturbational gene effects of the 13 ER⁺ essential genes (a), 44 HER2⁺ essential genes (c), and 29 TNBC essential genes (e). PPI network enrichment analysis for the identified essential genes in (b) ER⁺ (number of nodes: 13; number of edges: 11; average node degree: 1.69; avg. local clustering coefficient: 0.623); (d) HER2⁺ (number of nodes: 44; number of edges: 225; average node degree: 10.2; avg. local clustering coefficient: 0.625); and (f) TNBC (number of nodes: 29; number of edges: 141; average node degree: 9.72; avg. local clustering coefficient: 0.638). CRISPR-Cas9 perturbational gene effects data were retrieved from the Dependency Map database. Y-axis represents the perturbational effect scores for the indicated genes.

Figure 5

Prognostic significance of the identified TNBC therapeutic targets. Kaplan-Myer RFS analysis for ENO1 (a), FDPS (b), CCT6A (c), TUBB2A (d), and PGK1 (e) based on median expression in a cohort of 442 basal BCs. Log–rank p value is indicated on each plot. Expression of ENO1 (f), FDPS (g), PGK1 (h), CCT6A (i), and TUBB2A (j) as a function of TNBC molecular subtypes in an independent cohort of TNBC (n = 360). Anova p value is indicated on each plot.

Figure 6

Colony formation and viability assessment of TNBC cells models in response to ENO1 and FDPS depletion. CFU of MDA-MB-231 and BT-549 in response to FDPS and ENO1 depletion (a). Quantification of CFU in MDA-MB-231 and BT-549 in response to FDPS and ENO1 depletion (b). Data are presented as means ± S.E., n = 3. *** p < 0.0005. Upper panels show representative live (green) and dead (red) staining of MDA-MB-231 (c) and BT-549 (d) in response to FDPS and ENO1 depletion. (Scale bar = 1000 μM). Quantification of the number of dead (red) cells under each treatment condition is shown in lower panels. Data are shown as mean number of dead cells in three different fields (10×) ± S.D, n = 3. *** p < 0.0005.

Figure 7

Targeted depletion of FDPS and ENO1 reduces migration and organotypic growth of TNBC cell models. Cell migration (scratch assay) for MDA-MB-231 (a) and BT-549 (b) TNBC cell models in response to FDPS and ENO1 depletion. Quantification of relative (%) wound area under each treatment condition is shown in the right panels. Inhibition of 3D organoid formation of MDA-MB-231 (c) and BT-549 (d) TNBC models in response to FDPS and ENO1 depletion. (Scale bar = 1000 μM). Quantification of number of organoids under each treatment condition relative to control is shown in the right panels. Data are shown as mean relative number from three different fields (10×) ± S.D, n = 3. *** p < 0.0005.

Figure 8

GO enrichment analysis associated with elevated ENO1 and FDPS expression in TNBC. GO enrichment tree in a cohort of 360 TNBC patients divided according to median ENO1 (a) and FDPS (b) expression. Enrichment p value is indicated for each functional category.