Classification of breast tumors by triple hormone receptor protein expression

The triple-hormone receptor (THR) categorization is based on protein expression of ER, AR, and VDR, assessed by immunohistochemical (IHC) staining of formalin-fixed paraffin-embedded (FFPE) BrCa tissue microarrays (TMA) as described before [3, 4], producing four subgroups: THR-0, THR-1, THR-2, and THR-3, representing 7%, 11%, 28%, and 54% of BrCa, respectively (Fig. 1A). We previously showed in KM survival analysis that BrCa with fewer hormone receptors is associated with shorter OS [3] with a statistically significant hazard ratio (HR) in multivariate analysis: THR-0 (HR = 6.9, CI: 3.3–14.3, n = 104); THR-1 (HR = 5.3, CI: 2.7–9.9, n = 185); THR-2 (HR = 2.9, CI: 1.6–5.2, n = 429); and THR-3 (HR = 1.0, n = 998) in the Nurses’ Health Study dataset (n = 1716) [3].

Fig. 1

Breast cancer classification based on triple-hormone receptor (THR) expression. A Immunohistochemical (IHC) staining with ER, AR, and VDR of tissue microarrays (TMAs) from breast cancer patients (top) identifies four distinct subtypes (bottom). Hormone receptor positive tumors were identified as those with > 1% protein expression. B Heatmaps showing the expression of the top 300 (left) and top 50 (right) differentially expressed genes between THR-[0/1] and THR-[2/3] cell lines in the Cancer Cell Line Encyclopedia (CCLE) dataset. C Heatmap showing the expression of the THR-50 genes in human samples from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort. The heatmap shows the samples annotation including ER status measured by IHC, clinical 3-gene classifier groups (based on ER, HER2, and MIB-1), histological grade, and PAM-50 groups. Red = high expression, blue = low expression. D–E Kaplan–Meier survival plots show the difference in recurrence-free survival (RFS) (D) and overall survival (OS) (E) between METABRIC samples predicted as 0 (low-risk) and 1 (high-risk) by THR-50. High-risk samples have significantly worse RFS (HR = 1.5, 95%CI: 1.3–1.7, p < 0.0001) and OS (HR = 1.7, 95%CI: 1.5–1.9, p < 0.0001) compared to low-risk samples. Survival time is in months. The hazard ratios and 95% confidence intervals are shown. THR: triple-hormone receptors; HR: hazard ratio; CI: confidence interval

Derivation of a triple hormone receptor mRNA signature: THR-50

Utilizing data from the CCLE dataset, we developed an mRNA signature that can distinguish between THR-0/1 and THR-2/3 BrCa [27]. We selected the top 50 most significant genes (lowest p-value) from the 600 differentially expressed mRNAs to investigate further, referred to as THR-50 signature hereafter, which allowed us to examine the THR-IHC index in publicly available BrCa gene expression datasets (Fig. 1B and Additional Files 2 and 4).

Validation of THR-50 in human breast cancer

Analysis of the METABRIC cohort (n = 1904) [28] demonstrates that the median expression of THR-50 divides BrCa into two major clusters (Fig. 1C). The genes associated with THR-[2/3] CCLE (Fig. 1B) are highly expressed in ER-positive BrCa within the METABRIC dataset (Fig. 1C), as expected. THR-50 is also significantly associated with RFS (HR = 1.5) and OS (HR = 1.7) (Fig. 1D–E). Additionally, even after adjusting for important variables such as age, tumor stage, and grade using a multivariate Cox proportional hazards model, THR-50 remains significantly associated with RFS and OS (Table 1), underscoring the independent prognostic value of this signature.

Table 1 Multivariate survival analysis of THR-50Analysis of THR-50 across breast cancer subtypes

Next, we explored the association of THR-50 with survival outcomes across BrCa subtypes. Using the THR-50-derived risk score in the METABRIC cohort, we stratified patients into four equal groups with quartile 1 (Q1) and quartile 4 (Q4) representing the lowest and highest risk, respectively. We observed that patients in Q4 have significantly worse RFS compared to Q1 across ER-positive, ER-negative and HER2 + BrCa (Figure S2A, Additional File 5), as well as in the Luminal A (HR = 1.5) and Luminal B (HR = 2.4) BrCa subtypes (Figure S2B, Additional File 5). The distribution of THR-50 risk scores across different BrCa subtypes, defined by the PAM-50 and clinical 3-gene classification schemes, is illustrated in Figure S3, Additional File 5.

These findings were validated in the KMP cohort (n = 2,032, 50 studies) [29], where patients were stratified based on the average expression of THR-50, into low- and high-expression groups. We found that patients with low expression of THR-50 exhibit worse OS compared to those with high expression across ER + , AR + , ER-, Lum-A, Lum-B, HER2 + , and lymph-node positive (LN +) BrCa subtypes (Figure S4, Additional File 5). These results underscore the utility of the THR cell-of-origin signature in stratifying the risk of BrCa patients across diverse molecular and clinical subtypes.

THR-50 demonstrates promising performance relative to existing prognostic biomarker tests

Multigene biomarker tests such as Oncotype DX, PAM-50, MammaPrint, and EndoPredict are recommended by the American Society of Clinical Oncology (ASCO) for ER-positive, HER2-negative, and lymph node-negative BrCa [48,49,50]. However, they are not generally recommended by ASCO for ER-negative, HER2-positive, lymph node metastatic (> N1), or treated BrCa [51, 52].

In the KMP cohort (n = 2032), we evaluated the performance of THR-50 alongside PAM-50, MammaPrint, and Oncotype DX by categorizing patients based on the average expression of signature genes using optimal cutoffs. Remarkably, high average expression of THR-50 is significantly associated with better RFS compared to low expression in overall BrCa (HR = 2.04). Similarly, PAM-50 shows a significant association with RFS in overall BrCa (HR = 1.4), while MammaPrint and Oncotype DX do not (p = 0.052 and 0.13, respectively) (Figure S5, Additional File 5).

Next, we investigated the prognostic performance of these signatures across various BrCa groups. Noticeably, THR-50 demonstrates significant associations with RFS in multiple BrCa subgroups, including lymph node-positive (HR = 2.4), AR-positive (HR = 2.9), grade 2 (HR = 2.4), and grade 3 (HR = 1.6) BrCa (Fig. 2A). Similarly, PAM-50 demonstrates significant associations with RFS in the same groups, although with a slightly lower prognostic power compared to THR-50, except for grade 3 BrCa (p = 0.29) (Fig. 2A).

Fig. 2

THR-50 is significantly associated with recurrence-free survival (RFS) across different breast cancer clinical groups, outperforming established tests in the KMP cohort. A Kaplan–Meier survival plots comparing the prognostic power of THR-50 with PAM-50 using RFS in lymph-node positive, androgen receptor (AR) positive (AR +), grade 2 and grad 3 breast cancer. The analysis uses an independent validation cohort (KMP) comprising 2,032 samples from 50 gene expression datasets. Low (black line) and high (red line) expression groups are defined based on optimum cutoffs of the average expression levels of all signature genes. The reported p-values are derived from the log-rank test. The hazard ratios (HR) along with their corresponding 95% confidence intervals (CI) are shown. THR-50 results: lymph-node positive (HR = 2.4, CI:1.9–3.1, p = 1.5e-13), AR + (HR = 4.5, CI:2.0–11.0, p = 0.0001), Grade 2 (HR = 2.0, CI:1.1– 4.0, p = 0.02), Grade 3 (HR = 1.6, CI:1.1–2.5, p =  p = 0.02). PAM-50 results: lymph-node positive (HR = 1.4, CI:1.1–1.9, p = 0.0026), AR + (HR = 2.0, CI:1.1–3.9, p = 0.025), Grade 2 (HR = 1.8, CI:1.1–3.0, p = 0.021), Grade 3 (p = 0.29). B Kaplan–Meier survival plots comparing the RFS between patients with low and high average expression of THR-50 genes across different PAM-50 groups: Lum-A (HR = 2.1, CI = 1.5–3.1, p = 1.4e-05), Lum-B (HR = 1.8, CI = 1.4–2.3, p = 1.2e-05), HER2-like (HR = 2.3, CI = 1.6–3.2, p = 5.2e-07), basal-like (HR = 2.5, CI = 1.8–3.6, p = 7e-08). The plot also shows Oncotype DX (ONC-21), and MammaPrint (MAM-70) HR, 95% CI, and p-values. CI: 95% confidence interval. HR: hazard ratio

THR-50 identifies significant prognostic subgroups even within PAM-50 categories, revealing distinct RFS outcomes, in Luminal A (HR = 2.2), Luminal B (HR = 1.8), HER2-like (HR = 2.3), and basal-like BrCa (HR = 2.5) (Fig. 2B).

In contrast, MammaPrint demonstrates significance only in Luminal B (HR = 1.3) and HER2 + BrCa (HR = 1.5) subtypes. Oncotype DX, in comparison, shows significant associations with RFS across PAM-50 subtypes, but with reduced prognostic efficacy compared to THR-50, except for basal-like BrCa, where it did not reach significance (p = 0.15) (Fig. 2B).

These results indicate that THR-50 exhibits a significant prognostic value across diverse BrCa subtypes, unlike currently available tests.

Derivation and validation of THR-70

The results shown above using THR-50 suggest that the CCLE THR signature can be used to filter human tumor gene expression data. Therefore, fine-tuning of the THR signature was conducted by overlapping the cell line expression profiles with human tumor tissue cohort, comprising 855 BrCa cases (BC855) [21, 32,33,34, 53].

Analysis of BC855 cohort reveals that THR categories differ from PAM-50 subtypes (Fig. 3A); demonstrating that each THR group encompasses diverse proportions of all six PAM-50 subtypes. For instance, the THR-1 cohort comprises Luminal A (15.7%), Luminal B (18.8%), HER2-enriched (18.5%), Claudin-low (13.9%), Normal-like (8.7%), and Basal-like (24.4%) PAM-50 clusters, indicating a balanced distribution (Fig. 3A). This distribution suggests that the poor outcomes in THR-1 are not attributed to a single PAM-50 cluster.

Fig. 3

Development and validation of THR-70. A Heatmap showing the expression of the three hormone receptors ER, AR, and VDR across the different triple hormone receptor (THR) groups (top). The Pie charts (bottom) show the percentages of PAM-50 subtypes in the different THR groups in the BC855 cohort. B Venn diagram showing the genes in common between the top differentially expressed genes (DEGs) between the THR-0/1 and THR-2/3 groups in the CCLE and BC855 cohorts, using p < 0.05 as a cut-off. The THR-70 signature comprises the top 70 DEGs in common between both cohorts based on SAM-fold expression. C Heatmap of the expression of THR-70 genes in normal breast epithelial clusters reported in Bhat-Nakshatri et al. Expression levels are z-score transformed. D Violin plots comparing the enrichment of THR-70 across normal breast epithelial clusters identified by Kumar et al. Signature scores (normalized U statistics between 0 and 1), shown on the Y axis, were computed using UCell. E–F Kaplan–Meier survival plots in the METABRIC cohort comparing the overall survival (OS) between patients predicted as low (Q1) and high-risk (Q4) by THR-70. The high-risk samples have significantly worse OS compared to low-risk samples in the PAM-50 basal (HR = 2.2, 95%CI: 1.2–4.1, p = 0.01), Claudin-low (HR = 6.6, 95%CI: 2.5–17.2), Luminal A (HR = 2.9, 95%CI: 2.1–4.0, p < 0.0001), and Luminal B (HR = 4.2, 95%CI: 2.5–6.9, p < 0.0001) (E). Additionally, in clinical 3-gene classifier ER-/HER2- (HR = 2.7, 95%CI: 1.6–4.5, p < 0.0001), ER + /HER2- high proliferation (HR = 4.8, 95%CI: 2.9–8.0, p < 0.0001), and ER + /HER2- low proliferation (HR = 2.9, 95%CI: 2.1–4.1, p < 0.0001) (F). Survival time is in months. Hazard ratios (HR) and 95% confidence intervals (CI) are shown. Statistically significant HRs are highlighted in red

We found that 190 THR-associated mRNAs in cell lines (CCLE) overlap with THR signature in human BrCa tumors (BC855) (Fig. 3B). THR-70 refers to the top 70 genes identified by SAM fold change and p-value ranking among these 190 genes (Additional Files 3 and 4).

Next, by integrating THR-70 with the recently described single nucleus transcriptome data of healthy breast tissues [35], we show that THR-70 gene expression is enriched in normal human breast tissue. Interestingly, different THR-70 genes are enriched in proliferating (LASP-AP and LASP-BL) versus hormone sensing (LHS-HSα and LHS-HSβ) luminal epithelial breast cells (Fig. 3C).

Previously, we reported that the biomarker profiles of most human breast cancers, including the TNBC and basal-like cancers, are similar to normal luminal breast epithelium with the majority of tumors enriched for signatures derived from LHS and LASP cells [35, 54]. Consistent with this, we found that THR-70 is not enriched in normal basal-myoepithelial breast cells (BM_BAα and BM_BAβ) (Fig. 3C). Similar results are observed in the study by Kumar et al. [38], in which THR-70 is more enriched in luminal hormone-responsive cells (LummHR-SCGB, LummHR-active, and LummHR-major) compared to the luminal secretory and basal cells (Fig. 3D).

Having observed that THR-70 contains a breast cell-of-origin signature, we examined it in the KMP dataset revealing prognostic subgroups within PAM-50 categories including in Luminal A (HR = 1.6), Luminal B (HR = 2.0), HER2-like (HR = 2.1), and basal-like BrCa (HR = 2.2), as well as in lymph node-positive (HR = 2.8), AR-positive (HR = 1.6), grade 2 (HR = 2.5), and grade 3 (HR = 1.6) BrCa (Figure S6).

Next, we stratified patients within the METABRIC cohort based on THR-70, utilizing calculated risk scores, and observed that patients in the highest-risk category (Q4) have worse RFS and OS compared to those in the lowest-risk category (Q1) across all BrCa subtypes except HER2 + (Figs. 3E-F and S7-S8, Additional File 5).

To validate these findings in another dataset, we examined the THR-70 signature in the Meta-10 cohort, which comprises samples from ten different gene expression datasets [30] (Figures S9 and S10, Additional File 5), where THR-70 demonstrates significant association with RFS (HR = 2.5), distant metastasis-free survival (DMFS) (HR = 3.8), and survival in lymph node-positive (HR = 9.7), lymph node-negative (HR = 3.4), treated (HR = 4.4), and untreated (HR = 2.3) patients (Fig. 4).

Fig. 4

THR-70 is prognostic in the SurvExpress Meta-10 cohort, comprising samples from 10 different datasets. THR-70 shows a significant association with recurrence-free survival (RFS, HR = 2.5 CI: 2.0—3.1, p = 3.7e-16) (A) and distant metastasis-free survival (DMFS, HR = 3.7 CI: 2.7—5.6, p = 6.7e-15) (B). It is prognostic in both lymph node positive (LN + , HR = 9.7 CI: 5.0—18.7, p = 9.6e-12) (C) and negative (LN-, HR = 3.3 CI: 2.4—4.5, p = 2.7e-15) (D) disease, as well as in patients treated with endocrine therapy post surgery (ETPS, HR = 4.3 CI: 3.0—6.3, p = 1.4e-14) (E), and those who did not receive neoadjuvant treatment (NTPS, HR = 2.3 CI: 1.8—2.8, p = 8.7e-17) (F). 95% confidence interval (CI); hazard ratio (HR)

Cumulatively, these results indicate that ASCO-recommended multigene biomarker tests exhibit varying prognostic efficacy across different BrCa subgroups. In contrast, THR signatures demonstrate consistent prognostic power across BrCa subtypes. These findings are notable as conventional prognostic signatures rarely correlate with such diverse aspects of tumor biology across multiple datasets (Figures S7-S10, Additional File 5), suggesting that a cell-of-origin signature may influence various facets of the tumor phenotype.

THR-50 and -70 are associated with hormone and immune gene set signatures

To gain deeper biological insights into the functional roles of THR-50 and THR-70, we conducted gene set variation analysis (GSVA). Consistent with their derivation based on hormone receptor protein expression, both THR signatures are enriched for AR and ER pathways, which are negatively correlated with PAM-50, MammaPrint, and a proliferation signature (PCNA-131) (Fig. 5A).

Fig. 5

The THR signatures are enriched in ER, AR, and immune pathways in gene set variation analysis (GSVA). A Heatmap showing the Spearman correlation coefficients between different breast cancer signatures (rows) and key cancer-associated pathways and biological processes (columns). Positive correlation (red), negative correlation (blue), p-value < 0.05 (*), false discovery rate (FDR) < 0.05 (#). B Immune cell type enrichment score heatmap (rows) in different breast cancer signatures (columns). Positive enrichment (red) and negative enrichment (blue)

Traditional tissue-based prognostic signatures often reflect differences in tumor proliferation due to marked outcome differences between high-grade/high-proliferation tumors with poor outcomes and better-outcome tumors with low proliferation. It has been reported that proliferation-related genes are overrepresented in 88% of the BrCa prognostic signatures examined [55], and their removal substantially diminishes the prognostic efficacy across a majority of the 47 published signatures [56]. This suggests that many established BrCa prognostic tests may primarily operate as surrogate markers for proliferation [57]. Accordingly, the most notable positive associations observed in PAM-50 and MammaPrint through GSVA are with cell cycle and apoptosis pathways (Fig. 5A). In contrast, we aimed to reduce this proliferation bias by filtering THR signature from human tumors (BC855) with cell lines (CCLE-600) that exhibit uniformly high proliferation rates. Combining tumor and cell line data ensured that THR-70 is less influenced by proliferation effects (Fig. 5A), cell line artifacts, and non-tumor signals from tissue samples.

Notably, analysis of immune enrichment scores reveals that THR signatures are considerably associated with tumors with higher levels of central memory, gamma delta (γδ), CD4 + T, Th17, T follicular helper (Tfh), NK, and MAIT cell infiltrates (Fig. 5B). In contrast, PAM-50, MammaPrint, and PCNA-131 signatures are associated with tumors that are infiltrated with myeloid lineages such as neutrophils, dendritic cells, and monocytes (Fig. 5B). Some of these immune cells are known to influence BrCa outcomes, with CD8 + T cells generally associated with better outcomes and Treg/Th17 cells linked to poorer outcomes. Interestingly, γδT and Tfh are involved in anti-tumor cytotoxicity and antibody generation, respectively [58], suggesting that THR signatures may interact with both anti-tumor and pro-tumor immune infiltrates. Further research is needed to fully elucidate this immunological landscape, underscoring the intricate relationship between breast cell-of-origin and immune response in BrCa. In summary, these results indicate that THR signatures capture aspects of BrCa tumor biology that are not captured by standard prognostic tests.

Unsupervised clustering of breast cancer samples using THR signatures reveals distinct subtypes

We next explored the utility of the THR signature for de novo classification of BrCa. Using unsupervised clustering, we grouped samples in the METABRIC cohort based on their expression of THR-70, and subsequently analyzed the specific survival rates associated with each identified group. Our findings reveal that THR-70 divides BrCa into five distinct clusters: E1, E2a, E2b, E3, and PQNBC (Fig. 6A).

Fig. 6

Unsupervised clustering based on THR-70 uncovers five distinct breast cancer groups in the METABRIC cohort. A Heatmap showing the expression of THR-70 genes in the METABRIC cohort. Five distinct groups were identified: E1, E2a, E2b, E3, and PQNBC. B Kaplan–Meier survival plots comparing the 20 years recurrence-free survival (RFS) rates between different breast cancer groups identified by the by the 3-gene (ER, HER2, Mib-1 IHC) classifier (clinical), THR-50 combined with i20 (THR-50i), and THR-70 combined with i20 (THR-70i) signatures. THR: triple-hormone receptor. Pi + : pentaplex-negative (ER, PR, AR, VDR, and HER2), immune-positive tumors. Pi-: pentaplex-negative, immune-negative tumors. Survival time is in months

These clusters encompass a spectrum of ER-positive (E1-E3) clusters (Figure S11A, Additional File 5). Based on their overlapping survival curves, we consolidated E2a and E2b into a unified prognostic group named E2 (Figure S11B-C, Additional File 5). The pentaplex-negative and quadruple-negative BrCa (PNBC and QNBC) clusters include breast tumors that are negative for (ER, PR, ±HER2, ±AR, and ±VDR) (Fig. 6A).

The THR heatmap clusters differ from existing categories of BrCa. For instance, while the THR-E clusters predominantly consist of Luminal A (Lum-A), Luminal B (Lum-B), and low-grade tumors, they also encompass other subtypes like HER2 +, claudin-low, and high-grade tumors, albeit less frequently (Figure S12, Additional File 5).

Likewise, the PQNBC cluster contains multiple PAM-50 subtypes, including basal-like (49.7%), claudin-low (34.5%), and HER2-like (14.7%) (Figure S12, Additional File 5). These results suggest that the THR groups represent distinct classifications rather than merely renaming existing categories.

THR-70 clusters can be further stratified utilizing an immune signature

An immune signature consisting of 20 genes (referred to as i20, Additional File 4) was used to further divide the PQNBC cluster into two subgroups (PQNBC.i + and PQNBC.i-), with the PQNBC.i + subgroup characterized by higher immune infiltration than the PQNBC.i- subgroup. The combined THR-immune signatures are referred to as THR-50i and THR-70i.

Notably, when compared to the clinical 3-gene classification scheme, both THR-50i and THR-70i KM charts show significantly fewer survival curve crossovers (Fig. 6B). In contrast, ER-/HER2- tumors, typically considered poor outcome subtypes, exhibit survival curve overlap with ER-positive clusters (ER + LP/HP) (Fig. 6B). Similar findings were observed when comparing THR-70i to PAM-50, with THR-70i demonstrating clearer separation between different BrCa groups compared to PAM-50 (Fig. 7). These results align with previous studies showing multiple KM curve crossovers among PAM-50 subtypes [24, 59,60,61,62,63,64], highlighting that THR-immune signatures provide improved separation of outcome groups with reduced overlap.

Fig. 7

THR-70 improves the identification of ER-negative (ER-) and ER-positive (ER+) subgroups with distinct survival rates compared to the clinical 3-gene classifier and PAM-50. A Kaplan–Meier (KM) survival plots comparing the 20-year recurrence-free survival (RFS) in ER-negative breast cancer groups identified by clinical 3-gene classifier: HER2 + HR = 1.5, 95%CI: 1.1–2.0, p = 0.001 vs. TNBC (left panel); PAM-50 classifier: basal HR = 1.6, 95%CI: 1.2–2.2, p = 0.01 vs. claudin-low (middle panel), and THR-70i: PQNBC.i- HR = 15.7, 95%CI: 8.5–29.0, p < 0.0001 vs. PQNBC.i + (right panel). B KM plots comparing the 20-year RFS in ER + groups identified by clinical 3-gene classifier: ER + HP HR = 1.7, 95%CI: 1.4–2.1, p < 0.0001 vs. ER + LP (left panel), PAM-50: luminal B HR = 1.8, 95%CI: 1.5–2.2, p < 0.0001 vs. Luminal A (middle panel), and THR-70i: E1 HR = 2.1, 95%CI: 1.6–2.7, p <  p < 0.0001; E2 HR = 1.6, 95%CI: 1.3–2.0, p = 0.0001, VS. E3 (right panel). Survival time is in months. The hazard ratios (HR) and 95% confidence intervals (CI) are shown. HP: high proliferation, LP: low proliferation

As a cell-of-origin signature, THR genes exhibit mutations in less than 1–2% of BrCa cases. In contrast, PAM-50 signature genes have mutations in 35% of BrCa cases (Figure S13-14, Additional File 5). It has been proposed that early mutations might become redundant or non-essential as tumors progress, potentially altering the prognostic relevance of mutation-based signatures over time. Conversely, cell-of-origin signatures like THR may maintain stability throughout the tumor's lifespan.

Interestingly, i20 was also able to divide the ER-/HER2- group, defined by the clinical 3-gene classification scheme, into two distinct subgroups with significantly different survival rates (Figure S15B-C, Additional File 5), suggesting that i20 immune signature can further enhance the granularity of existing classifiers.

Stratification of ER + , TNBC, and HER2 cancer with THR.i signature

In ER-negative BrCa, we observed a 1.5-fold difference in survival probability using the clinical 3-gene (ER/HER2/Mib-1) and PAM-50 classification methods. In contrast, there is a 15-fold difference in survival probability between predominantly ER-negative PQNBC.i + and PQNBC.i- cohorts (HR = 15.7, 95%CI: 8.5–29.0). Therefore, compared to the 3-gene classifier and PAM-50, the THR-70i signature demonstrates a ten-fold improvement in distinguishing ER-negative breast tumor subtypes with markedly poor and favorable outcomes (Fig. 7A).

For ER-positive BrCa, while the 3-gene and PAM-50 methods identify two ER-positive subtypes with hazard ratios differing by 1.5- to 1.8-fold, THR-70 delineates three distinct ER-positive clusters (E1, E2, and E3) with a survival range differing by 2.1-fold (Fig. 7B).

Combining cell-of-origin, immune, and genetic biomarkers

We previously demonstrated that HER2-amplified BrCa does not align with a specific cell subtype in normal human breast. Consistent with their pathogenesis, HER2 + cancers exhibit marker profiles spanning various normal luminal cell types [3]. Therefore, we did not anticipate that HER2 + tumors would form a distinct cluster based on cell-of-origin signatures, which indeed was observed in THR cluster heatmaps (Fig. 6). Thus, to further stratify BrCa we coupled THR-70i with HER2 + (THR-70Hi) that identified six BrCa groups with different survival estimates: PQNBC.i-, PQNBC.i + , E1, E2, E3, and HER2 + (Figs. 8A and S15A-C, Additional File 5). The PQNBC subtype includes triple-negative (ER, PR, HER2) breast cancers (TNBC) that may lack either AR or VDR (quadruple-negative) and both AR and VDR (pentaplex-negative).

Fig. 8

THR-70 coupled with immune signature (i20) and HER2 (THR-70Hi) captures more granular breast cancer groups compared to PAM-50. A Kaplan–Meier (KM) survival chart shows 20-year recurrence-free survival (RFS) of different patient subgroups identified by THR-70, i20, and HER2 classifier (THR70-Hi): E3 (purple), E2 (black), E1 (blue), HER2 + (yellow), PQNBC.i- (green) and PQNBC.i + (red). PNBC subtype includes breast cancers that are negative for ER, PR, HER2, AR and VDR. QNBC subtype includes breast cancers that are negative for ER, PR, HER2, ∓ AR or VDR. B KM survival chart shows 20-year RFS of different patient subgroups identified by PAM-50: Lum-A (purple), Lum-B (blue), HER2-like (yellow), basal-like (green), and claudin-low (red). C Multivariate analysis of Hazard ratios (HR) and 95% confidence interval (95% CI) for RFS by THR-70Hi and PAM-50 breast cancer groups using a Cox proportional hazards model

The survival curves of THR-70Hi groups generally do not overlap, demonstrating survival differences ranging up to 5.8-fold in both univariate (Fig. 8A) and multivariate (Fig. 8C) survival analyses. In comparison, the PAM-50 clusters exhibit up to a 3.6-fold survival difference range among its groups (Fig. 8B). However, the basal-like subtype's survival curve crosses over HER2, Luminal B, and Luminal A subtypes around 5-, 10-, and 20 years, respectively, complicating its assessment (Fig. 8B).

Next, we examined the relationship between THR and HER2 in more detail. First, we note that HER2 status has no effect on RFS in the THR-70 PQNBC cluster (p = 0.28; Figure S16A, Additional File 5), emphasizing the domina