Association of pathology-assessed cellularity with clinico-pathological features and patient’s outcomeWhile tumour purity can significantly impact quantification and interpretation of trancriptomic data, it is important to appreciate that both intrinsic and extrinsic factors can influence it. To obtain insights into the relative impact of these two aspects on breast cancer sample tumour content, we investigated whether pathology-assessed tumour cellularity is associated with biologically and clinically relevant features in breast cancer. Associations were evaluated in the TCGA, Metzger-Filho, Park and NA-PHER2 datasets, where tumour purity quantification by expert pathologists was available (Table 1). We found multiple significant associations as reported in Fig. 1a. In the TCGA, the invasive lobular cancers (ILC) had lower purity than other subtypes (Fig. 1a and Supplementary Fig. 1A). At the same time, in the Metzger-Filho dataset, including only ILCs, cellularity was significantly different between distinct ILC subtypes, with the lobular Classic subtype showing the lowest cellularity (Supplementary Fig. S3B). The proliferation marker Ki67 had a weak correlation with cellularity in the Metzger-Filho dataset (ρ = 0.19, p = 0.047), but the same association was not confirmed in Park and NA-PHER2 datasets (Fig. 1a, Supplementary Figs. 2M, 3C and 4B). Notably, different subtypes were included in the different datasets. On the contrary, high grade was consistently associated with higher cellularity (Fig. 1a, Supplementary Figs S2A, S3A). Breast cancer subtypes, either defined by the combination of ER and HER2 status or according to PAM50 classification, showed a significant association with cellularity in both TCGA and Park datasets. In particular, triple-negative or basal-like tumours had, on average, the highest cellularity (Fig. 1a, Supplementary Figs. 1E, F and 2E, F). Finally, in NA-PHER2 samples, stromal tumour-infiltrating lymphocytes (sTILs), but not intraepithelial tumour-infiltrating lymphocytes (iTILs), were negatively correlated with tumour cellularity (ρ = −0.33, p = 0.018) (Fig. 1a and Supplementary Fig. 4C, D).Fig. 1: Association of pathologist’s estimated tumour purity with molecular and clinico-pathological variables in breast cancer.a Landscape of association of available molecular and clinico-pathological variables with cellularity in four breast cancer datasets: TCGA (n = 1073), Metzger-Filho (n = 117), Park (n = 112) and NA-PHER2 (n = 52) pre-treatment samples. (* association between purity and continuous variables was assessed by Spearman’s correlation, association with two categorical groups was assessed by Student’s t test, and association with multiple categorical groups was evaluated by one-way ANOVA). b Variance component analysis (VCA) for each dataset computed for samples with no missing information (TCGA, n = 690; Park, n = 72; Metzger-Filho, n = 111; NA-PHER2, n = 52). The analysis estimated the proportion of total variance explained by the provided variables. c Forest plot of Cox regression univariate analysis evaluating association of cellularity with overall survival in TCGA (n = 1073) and Metzger-Filho (n = 117) datasets. Samples were evaluated overall and stratified by subtype (TCGA: 426 Luminal, 162 HER2+, 113 TN; Metzger-Filho: 100 Luminal). d Cellularity changes in on-treatment biopsy (n = 86) compared to pre-treatment (n = 112) in the Park dataset. The impact of the timepoint on tumour purity was evaluated by Student’s t test and VCA. e Same analysis as in d for the NA-PHER2 dataset (n = 52, pre-treatment biopsy; n = 40, on-treatment biopsy).To quantify the contribution of each clinico-pathological variable to the overall tumour cellularity variance, we performed a variance component analysis (VCA) for each dataset (Fig. 1b). In line with the analysis above, the histological type and grade, together with the molecular subtypes, explained the highest percentage of total variance. However, 52–100% of cellularity variance in each dataset was not explained by the included variables.Next, we evaluated whether tumour cellularity was associated with distant metastasis free survival in the TCGA and Metzger-Filho datasets, overall and for each subtype defined by ER and HER2 status (Fig. 1c). Cellularity was significantly associated with prognosis only in lobular luminal cases of the Metzger-Filho dataset (p = 0.04).Finally, in the Park and NA-PHER2 datasets, transcriptomic profiles were obtained from pre- and on-treatment biopsies. Because of treatment-induced tumour cell death, an overall reduction in tumour cellularity could be expected and was observed (p = 0.009 and p = 0.003, respectively). The biopsy timepoint explained 6.12% of the variance in the Park dataset and 16.3% in the NA-PHER2 dataset (Fig. 1d, e).In summary, this analysis denoted that, in breast cancer, intrinsic tumour biology factors can affect tumour cellularity. However, over half of the variability observed in clinical specimens undergoing molecular characterization was not explained by the main clinico-pathological features and could be related to tumour sampling.Development of a Breast Cancer Purity Score (BCPS)Cellularity estimation by the pathologist is not always available and might not refer to the same tumour region undergoing molecular characterization. Consequently, we aimed at identifying a gene expression signature able to estimate tumour purity in a bulk transcriptomic analysis of clinical breast cancer samples. As detailed hereafter, to generate the BCPS we interrogated four distinct datasets: Bruna, NA-PHER2, NeoTRIP21 (surgical samples) and Metzger-Filho (Fig. 2 and Table 1).Fig. 2: BCPS generation.Workflow involving four distinct datasets leading to the definition of the BCPS. In the NA-PHER2 and Metzger-Filho datasets, the correlation between tumour purity and the expression values of each gene was computed. In the Bruna dataset, primary tumours were compared to matched patient-derived xenografts to identify candidate tumour-specific and stroma-specific genes, exploiting the loss of human stroma during engraftment. In the NeoTRIP dataset, the ROC curve AUC was estimated for each gene considering surgical samples with medium/high or low/no cellularity, as annotated by expert pathologists. By applying for each analysis the indicated thresholds, 5 tumour-associated genes and 4 TME-associated genes were selected to generate the BCPS.We correlated gene expression with tumour cellularity in the NA-PHER2 (n = 92) dataset, identifying 733 genes with a positive correlation with tumour cellularity (ρ > 0.35) and 1009 genes with a negative correlation (ρ < −0.35). The same analysis was performed in the Metzger-Filho (n = 117) dataset, identifying 49 genes with a positive (ρ > 0.35) and 284 genes with a negative (ρ < −0.35) correlation with cellularity. Hence, we considered surgical samples from the NeoTRIP dataset grouped into two categories based on expert pathologist evaluation: high/mid tumour content or low/no tumour content. We assessed the ability of each gene to distinguish between the two classes by ROC curve analysis, identifying 570 genes with AUC ≥ 0.7 and 330 genes with AUC < 0.3. Results from the analysis in Fig. 1 helped us defining the BCPS development strategy. Indeed, the association between cellularity and breast cancer subtypes, suggested performing the correlation with cellularity separately in datasets including only one subtype (Table 1).As a complementary strategy, in the Bruna dataset, transcriptomic profiles of patient-derived xenografts (PDXs) were compared with matched clinical samples from which the xenografts originated. Since the human stroma is completely lost during engraftment and replaced by mouse stroma24, genes expressed by the human tumour microenvironment are expected to be downregulated in such comparison. At the same time, tumour specific genes are expected to be similarly expressed or upregulated in PDXs. The fold changes for selected genesets, expected to be tumour specific or stroma-specific, are reported in Supplementary Figure 5A–C as proof of concept. In total, 236 stroma-specific (with logFC < −1.5) and 5679 tumour specific genes (with logFC > 0) were identified, candidate to be good reporter of tumour content. All the statistical estimates are available in Supplementary Data 1.From the candidate genes selected in each dataset (Supplementary Fig. 5D–G), we derived a consensus list of 5 tumour-associated (AP1M2, CDK5, PAFAH1B3, SLC25A10, SMG5) and 4 stroma-associated (CXCL12, IFFO1, MFAP4, TGFBR2) genes (Fig. 2). Their main known biological functions are summarised in Supplementary Data 2. This set of genes were used for a single-sample geneset enrichment analysis providing the BCPS, proportional to sample tumour purity, as characterised in the following section.BCPS evaluation of performance and comparison with ESTIMATE scoreTo evaluate the ability of the BCPS to estimate tumour sample purity and to compare its performance with the commonly used ESTIMATE score18, we interrogated five additional independent datasets: TCGA, Park, METABRIC, NeoTRIP (on-treatment samples) and Bianchini28 (Fig. 3 and Table 1).Fig. 3: BCPS validation and comparison with ESTIMATE score.Evaluation of the BCPS and comparison with ESTIMATE score in the TCGA, Park, NeoTRIP and Bianchini datasets. a Spearman’s correlation between the pathologist-estimated cellularity and either ESTIMATE or BCPS in TCGA (n = 1073). b Same analysis as in a for the Park dataset (n = 225). c ESTIMATE score and BCPS values measured in samples with high/medium tumour content or low/no tumour content in the NeoTRIP dataset (n = 219, on-treatment biopsy); two-sided Student’s t test. d ESTIMATE score and BCPS ability to discriminate between the two classes in c quantified by AUC. e ESTIMATE score and BCPS values measured in core biopsies (CBX) and matched fine-needle aspirations (FNA, n = 37 pairs) from the Bianchini dataset; two-sided Student’s t test. f ESTIMATE Score and BCPS ability to discriminate between the two classes in e quantified by AUC. g Example of KRT18 gene expression correction using the BCPS and linear regression to remove the impact of tumour purity. The Bianchini dataset was used. h Volcano plots of differential gene expression analysis between FNA and CBX samples of the Bianchini dataset. The analysis was performed without any correction and after normalising gene expression using the BCPS or ESTIMATE scores.In the TCGA and Park datasets, BCPS had a significantly higher correlation with cellularity than ESTIMATE (ρ = 0.34–0.37 vs ρ = 0.23–0.27 respectively; correlation difference p < 0.001) (Fig. 3a, b). In the METABRIC dataset, where cellularity was grouped into low, moderate and high classes, BCPS better discriminated between the groups and AUC for the prediction of high vs low cellularity was 0.75 for the BCPS and 0.62 for ESTIMATE (Supplementary Fig. 6A, B). ESTIMATE was strongly inversely correlated with a gene expression-based immune score29 (ρ = −0.87 and ρ = −0.90 in TCGA and Park, respectively), while cellularity and BCPS showed weaker correlations (Supplementary Fig. 6C, D). In the TCGA dataset, Aran et al.5 computed multiple RNA, DNA and methylation based predictors of tumour purity. We compared their estimates with the BCPS in the set of samples where all the metrics were available. BCPS showed the highest correlation with the pathologist’s cellularity (Supplementary Fig. 6E).The association with tumour content was also evaluated in the on-treatment biopsies of the NeoTRIP dataset (Fig. 3c, d). Samples were grouped into high/medium tumour content or low/no tumour content, as established by expert pathologists, and the ability of the two scores to discriminate between these two classes was quantified by Student’s t test and ROC curve analysis. BCPS better discriminated between the two groups, achieving a significantly higher AUC (BCPS AUC = 0.77, ESTIMATE AUC = 0.65, p < 0.001).We then evaluated the BCPS and ESTIMATE score in the Bianchini dataset. This dataset contains matched primary breast cancer samples obtained by either fine-needle aspiration (FNA) or core-biopsy (CBX). The first sampling procedure is known to enrich for tumour cells while CBX better preserves the stromal content. Indeed, both BCPS and ESTIMATE score were higher in the FNA samples (paired Student’s t test p = 1.1 × 10−8 and p = 9 × 10−6, respectively). ROC curve analysis highlighted a significantly higher AUC for the BCPS compared to ESTIMATE score (BCPS AUC = 0.85, ESTIMATE AUC = 0.78, p = 0.046) (Fig. 3e, f).In the Bianchini dataset, differences between matched FNA and CBX are expected to be primarily related to sampling differences affecting tumour content. Consequently, it represents a relevant setting where to quantify such bias and evaluate the ability of a purity score to correct for it. To this aim, we performed a paired differential analysis between FNA and CBX introducing either no correction or normalising the data using the BCPS or ESTIMATE score. BCPS-normalized data were obtained by taking the residuals of the linear regression models evaluating the relationship between BCPS and each gene. An example of gene expression before and after correction is shown in Fig. 3g. Differential analysis between FNA and CBX before correction identified 60 up-regulated genes and 409 down-regulated genes. Data correction using the ESTIMATE score reduced the number of differentially expressed genes to 89 up-regulated and 160 down-regulated, but only 40 up-regulated and 65 downregulated genes were observed in data corrected using the BCPS, showing the higher performance of the BCPS in adjusting for differences due to tumour content in clinical samples (Fig. 3h). This analysis indicated that in all situations where extrinsic factors are expected to largely overweigh intrinsic factors in affecting tumour purity, BCPS is a useful and effective tool to take tumour purity into consideration and correct for it.BCPS recapitulates cellularity associations with clinico-pathological factorsFor the TCGA and Park datasets, not used for the BCPS development, we evaluated the association with available clinico-pathological factors, as reported for the pathologist cellularity in Fig. 1.BCPS was significantly associated with the same variables that were significantly associated with cellularity in both the TCGA and Park datasets (Figs. 4a and 1a). VCA analysis was quantitatively different but qualitatively similar (Figs. 4b and 1b), and the BCPS was significantly lower in on-treatment compared to pre-treatment samples in the Park dataset, as observed for the pathologist’s cellularity (Figs. 4d and 1d). Finally, the BCPS was not significantly associated with survival in any subtype in the TCGA dataset (Figs. 4c and 1c).Fig. 4: Association of the BCPS with clinico-pathological factors in the TCGA and Park datasets.a Landscape of association of available molecular and clinico-pathological variables with cellularity in the TCGA (n = 1073) and Park (n = 112) datasets. (* association between purity and continuous variables was assessed by Spearman’s correlation, association with two categorical groups was assessed by Student’s t test, and association with multiple categorical groups was evaluated by one-way ANOVA). b Variance component analysis (VCA) for each dataset computed for samples with no missing information (TCGA, n = 690; Park, n = 72). The analysis estimated the proportion of total variance explained by the provided variables. c Forest plot of Cox regression univariate analysis evaluating cellularity association with overall survival in the TCGA (n = 1082) dataset. Samples were evaluated overall and stratified by subtype (426 Luminal, 162 HER2+, 113 TN). d Cellularity changes in on-treatment biopsy compared to pre-treatment in the Park dataset (T1 = 112, T2 = 86). The impact of the timepoint on tumour purity was evaluated by Student’s t test and VCA.While such associations would require external validation, these results support the validity of the BCPS in providing an estimate of tumour purity leading to conclusions similar to what could be drawn using the pathologist’s evaluation.Use of the BCPS for prediction of prognosis and response to treatmentA typical goal in transcriptomic data analysis of clinical samples is the identification of genes or signatures associated with specific clinical endpoints. Quantification of expression could be affected by tumour content7. We previously developed and validated an ER-related and a proliferation-related metagenes as predictors of long-term outcome in ER+/HER2− breast cancer13. Here we applied the metagenes to 2277 ER+/HER2− breast cancer samples from the Brueffer dataset (Table 1). A multivariate Cox model with interactions, including the ER- and proliferation metagenes, and the BCPS explained better the survival data than a bivariate model without the BCPS (likelihood ratio test p = 0.035). The improvement was confirmed by a higher c-index and higher 7-years AUC when the BCPS was included in the model (Fig. 5a, b).Fig. 5: Use of the BCPS in breast cancer prognostication.a, b In ER+/HER2− samples from the Brueffer dataset (n = 2277) 7-years overall survival was predicted using a multivariate Cox model with interactions including an ER and a Proliferation metagene with or without the BCPS. C-index (a) and 7-years AUC (b) were computed for the two models highlighting a performance improvement when the BCPS was included. c Association of the BCPS with pCR in on-treatment biopsies from the NeoTRIP dataset (n = 219). d BCPS quantified in the surgical samples of the NeoSPHERE trial was associated with DEFS. Two groups based on the BCPS median were identified and represented by Kaplan–Meier curves; differences were evaluated by log-rank test.Next, we focused on BCPS values in on-treatment biopsies of the NeoTRIP dataset, where the BCPS was significantly lower (p = 1.3 × 10−3) in cases eventually achieving pathological complete response or pCR (Fig. 5c). This indicates that, despite a possible sampling bias, evaluation of cellularity in on-treatment biopsies could help early prediction of patients responding or not to neoadjuvant treatments. Remarkably, this mimicked the predictive power of the pathologist’s evaluation in the same cohort, where 66.3% of patients with low/no cellularity had a pCR, versus only 35.3% of patients with mid/high cellularity in the on-treatment biopsy (p = 1.8 × 10−5).A third example where evaluation of the BCPS could provide valuable information is in surgical samples obtained after neoadjuvant treatment. Cellularity was included as one of the factors determining the Residual Cancer Burden score (RCB) and positively contributes to increase the score30. High RCB associates with worse prognosis. Here we evaluated the BCPS in post-treatment surgical samples from the neoadjuvant clinical trial NeoSPHERE22. In line with the prognostic role of pathologist’s cellularity as part of RCB score, the BCPS significantly stratified long-term patients’ risk by quantifying the amount of residual tumour left. The association remained significant either identifying 4 groups based on BCPS distribution quartiles (Supplementary Fig. 7B) or two groups using the median value (Fig. 5d).
