Differential Survival and Therapy Benefit of Patients with Breast Cancer Are Characterized by Distinct Epithelial and Immune Cell Microenvironments

Patients with breast cancer can respond differently to therapeutics even for tumors with the same subtype. Although multiple mechanisms may explain heterogeneity in therapy benefit, preclinical work has shown that the cellular composition of a tumor can influence metastasis and treatment efficacy (1–8). For instance, clinical outcomes upon different interventions often correlate with immune infiltrate (9–12). However, the role of different epithelial and nonimmune cell types for the survival and therapy outcome of breast cancer remains understudied, mainly because of the difficulty to detect these cell types in the samples that are available from clinical trials.

Deconvolution algorithms, such as CIBERSORT (13) and DWLS (14), which disentangle bulk mRNA sequencing data into contributions from individual cellular components, provide a good estimate of the cellular composition of tumors (13, 15). These methods have been used to confirm the relevance of tumor-infiltrating immune cells for breast cancer outcomes (16–18). However, the currently available reference profiles for deconvolution (13, 19) were mostly obtained from cells isolated from blood. Therefore, the reference profiles for other cell types in a tumor (e.g., cancer cells, epithelial cells, endothelial cells, muscle cells, and fibroblasts) are largely unknown. In addition, the expression profile of immune cells isolated from the blood and tumors can be different (20), and this can lead to less reliable estimations of the abundance of different cell types (15). Recently, reference profiles obtained from tumor-derived single-cell mRNA sequencing data have been used to deconvolve the cellular components of melanomas and ovarian tumors (14, 15, 19). However, good reference profiles for breast cancer cell types are not yet available.

In this study, we created references profiles of all the cellular components of human breast tumors through single-cell sequencing. We optimally detangle the cellular composition of tumors for which bulk mRNA sequencing data is available, a technique that from here on we refer to as tumor cell deconvolution (TCD). We applied TCD to determine the cellular composition of patients with breast cancer in two clinical trials for which both bulk mRNA sequencing data and clinical outcome are available: (i) 3,000 patients in the multicenter prospective trial Sweden Cancerome Analysis Network - Breast (SCAN-B; ref. 21), and (ii) 528 patients in the multicenter, prospective randomized trial Microarray Analysis in breast cancer to Tailor Adjuvant Drugs Or Regimens (MATADOR; refs. 22, 23). We identify and subsequently validate specific microenvironments that distinguish differential survival and therapy benefits of patients with breast cancer, which include immune and nonimmune cell types (see flow charts for the experimental setup in Supplementary Fig. S1).

We obtained the SCAN-B mRNA sequencing data and anonymized patient data through Gene Expression Omnibus (GEO) GSE81538.

The MATADOR (ISRCTN61893718) study is an open-label, multicenter randomized, phase III trial conducted in 29 centers in the Netherlands. Six hundred and sixty-four female patients with pT1–3, N0–3, M0 breast cancer were recruited onto the trial. The inclusion criteria were described in detail elsewhere (22, 23).

The study protocol and amendments were approved by the ethical committee of the Netherlands Cancer Institute. The study was conducted in agreement with Good Clinical Practice guidelines and with the Declaration of Helsinki. All patients provided written informed consent to participate in the trial and to use the tumor tissue removed at surgery for translational research. The Reporting Recommendations for Tumor Marker Prognostic Studies (REMARK) criteria were used to report this study (38).

Patients were randomly assigned (1:1) to 6 cycles of doxorubicin 60 mg/m² and cyclophosphamide 600 mg/m² every 2 weeks [dose dense (dd) AC] or docetaxel 75 mg/m², doxorubicin 50 mg/m², and cyclophosphamide 500 mg/m² every 3 weeks (docetaxel-doxorubicin-cyclophosphamide; TAC) by means of the automated ALEA system (FormsVision BV) using Pocock minimization technique (24). Randomization was performed centrally at the Netherlands Cancer Institute. G-CSF (pegfilgrastim) was given to all patients. Radiotherapy and endocrine therapy were given according to the contemporary Dutch guidelines.

Out of 664 patients, tumor tissue was available for gene expression analysis for 604 (90.9%) patients. Library preparation failed in 7 patients. Six of 597 sequenced samples did not meet the quality checks showing up as extreme outliers based on read counts per gene and in a principal component analysis (PCA). Sixty-two patients had a clinically low risk of recurrence according to the modified Adjuvant! Online (25) and would not receive adjuvant chemotherapy according to current guidelines. In this group, no events occurred. For one patient, survival data was missing. The remaining 528 patients were used to develop delta treatment score (DTS). For the analysis, we excluded the Her2⁺ patients.

The primary objective of the trial was to identify a gene expression profile for recurrence-free survival (RFS) benefit of either dose-dense or taxane-containing chemotherapy and to assess its predictive performance. RFS is defined as the time from randomization to locoregional recurrence, distant metastasis, or death by any cause, whichever occurred first. The secondary objective was to directly compare RFS, overall survival (OS), and toxicity of the two treatment arms and has been described elsewhere (22, 23).

The original statistical plan can be found in the MATADOR study protocol. During the course of the study, it became clear that the event rate was lower than expected. Therefore, an amendment was made to the protocol to use a cross-validation method instead of a separation in a training and a validation cohort.

Formalin-fixed paraffin-embedded (FFPE) material was collected between 2004 and 2012. RNA was isolated from FFPE tissue with a tumor cell percentage of at least 40% using the AllPrep DNA/RNA Mini Kit according to the manufacturer's instructions (Qiagen). Quantification and purity were measured using the NanoDrop 2,000 spectrophotometer (Thermo Fisher Scientific) and the 2,100 Bioanalyzer (Agilent Technologies). Libraries of cDNA were constructed with the TruSeq RNA Access Library Prep Kit (Illumina) and single-end sequenced using the HiSeq 2,500 (Illumina). Bulk mRNA sequencing data were analyzed using the RSEM package (26) and the reads were aligned to the hg38 reference genome using STAR (27) with default parameters. Transcripts per million values were used in further analyses.

To minimize batch effects that could arise during sample processing for RNA sequencing (RNA-seq), large sets of samples were processed within a relatively short time frame. This was a multicenter, randomized controlled trial, and samples were isolated and sequenced by mixing samples from different centers and samples that had been preserved in formalin for different lengths of time to control for possible batch effects. During the computational analyses, we also checked for batch effects caused by the isolation and sequencing runs by PCA, but none could be identified.

Following pathologic examination, tumor material was kept on ice in DMEM/F12 + GlutaMAX (GIBCO, Invitrogen Life Technologies), 1% HEPES, 1% penicillin–streptomycin from that point onwards. Tumors were mechanically dissociated using sterile scalpels and enzymatically digested with 25 μg/mL DNase I (Roche) and 15 Wünsch units TH Liberase/mL (Roche) diluted in PBS at 37°C for 35 minutes, followed by mashing through a 70-μm filter (BD Falcon) while adding DMEM/F12 + GlutaMAX (GIBCO, Invitrogen Life Technologies), 1% HEPES, 1% penicillin–streptomycin. Tumor cells were blocked in 80% FACS buffer [5 mmol/L EDTA in PBS supplemented with 5% FCS/20% serum mix (50/50 normal goat serum; monx10961, Monosan) and FcγII/III receptor blocking serum 2.4G2 (kind gift from L. Boon, Bioceros, the Netherlands)] for 10 minutes on ice before labeling with EpCAM PE (1B7, eBioscience), CD14-FITC (M-P-9, BD Biosciences) labeled antibodies for 40 minutes on ice. DAPI was added as a cell death marker, and cells were sorted either on a FACS AriaII Special Ordered Research Product (BD Biosciences), on a FACS Jazz (BD biosciences), or on a FACS Fusion (BD biosciences). A broad FSC SSC gate was followed by a gate excluding doublets, after which the DAPI-negative cells were selected. These DAPI-negative cells were sorted into 384-well plates and further processed for single-cell sequencing. To obtain enough epithelial (tumor) cells and macrophages, additionally EpCAM-positive and CD14-positive cells were selected and index-sorted into 384-well plates and further processed for single-cell sequencing. For 3 of 12 patients, only EpCAM-positive and CD14-positive cells were isolated and sequenced.

Single-cell mRNA sequencing was performed using Sort-seq as described in Muraro and colleagues (28). Single-cell libraries were sequenced on an Illumina NextSeq500 with 75-bp paired end reads. Read1 contains the cell barcode and Unique Molecule Identifier; read2 was mapped to the hg38 RefSeq transcriptome using Burrows-Wheeler Aligner with standard parameters. Single-cell analysis was performed using Seurat (29). For each of the sorts, the percentage of living cells varied between 1% and 5% of the total number of cells. In total, 14,592 cells from 12 different patients were sorted. Cells were required to have at least 2,000 unique protein coding transcripts and subsequently log-normalized to 10,000 transcripts per cell. With these filtering criteria, 6,580 were used for further analysis. The 5,000 most variable features were selected for clustering. Clusters were identified using the first 21 PCs with a resolution of 0.1; the cutoff for the number of components used in the clustering analysis was chosen based on the elbowplot form the Seurat package. The resolution of 0.1 yielded a relatively small number of clusters but with well-defined expression profiles. Uniform Manifold Approximation and Projections (UMAP) were also constructed using the first 21 PCs. We did not observe much variation in the expression profiles of the noncancerous cells between patients. This indicates that the batch effects between patients are not strong enough to cause separate clustering of the same cell type from different patients. Therefore, there was no need for correcting batch effects.

Bulk mRNA sequencing data for both the MATADOR cohort and the single-cell data were analyzed using the RSEM package (26); reads were aligned to the hg38 reference genome using STAR with default parameters (27). Transcripts per million values were used in further analyses. The count data for the SCAN B cohort was downloaded through GEO. The SCAN B data were mapped to the same reference genome (hg38) using the mapping algorithm Bowtie2. No batch correction algorithm was required nor used to match the MATADOR and SCAN B data sets.

To construct the reference profiles, all single cells belonging to the same Seurat cluster were averaged.

DWLS was used with default parameters using the single-cell derived reference profiles described above as reference profiles.

The cell type score (CTS_p) is a score based on the deconvolution-derived cell type abundances and is predictive for prognosis. The CTS_p is the hazard for every patient is calculated as follows.

where |{\beta _1}$ | is the collection of coefficients for the selected cell types, and |C{T_p}$ | is the abundance of the selected cell types in patient |p$ |⁠.

To find the optimal combination of cell types, a Cox proportional hazard model for each individual feature (⁠|{c_i}$ |⁠) is fitted. From this the |N$ | best features are selected, which are the |N$ | features that give the highest partial likelihood for the Cox proportional hazard model. For each of these |N$ | best features a second feature, |{c_2}$ |⁠, is selected by iteratively fitting a Cox model including the effects of |{c_1}$ | and |{c_2}.$ | Finally the |N$ | best combinations of two features are selected based on the partial likelihood of the Cox proportional hazard models. This process is repeated until all features, |{c_1}..{c_n}$ |⁠, are included in the model, resulting in many models each with a different combination of input features. Finally, the model with the optimal combination of features (⁠|x$ |⁠), |{c_1}..{c_x}$ |⁠, is selected, which is the model with the lowest cross validated P value. Association of the |CT{S^{}}$ | with RFS was tested in a Cox proportional hazard model corrected for the clinicopathologic features treatment, tumor size, lymph node status, pathologic subtype, histologic grade, age, and type of surgery.

The DTS is a score that is predictive for chemotherapy benefit and is calculated as follows. First, the optimal combination of features (cell types) is determined using a forward feature selection routine, where the complete set of possible features (⁠|n$ |⁠) are the 15 cell types predicted by DWLS. A Cox proportional hazard model for each individual feature (⁠|{c_i}$ |⁠), in combination with |T$ | and the interaction between the feature and |T$ | is fitted. From this, the |N$ | best features are selected, which are the |N$ | features that give the highest partial likelihood for the Cox proportional hazard model. For each of these |N$ | best features a second feature, |{c_2}$ |⁠, is selected by iteratively fitting a Cox model including the main effects of |{c_1}$ |⁠, |{c_2}$ |⁠, and |T$ |⁠, and the pairwise interactions between |c$ | and |T.$ | Finally the |N$ | best combinations of two features are selected based on the partial likelihood of the Cox proportional hazard models. This process is repeated until all features, |{c_1}\!.. {c_n}$ |⁠, are included in the model, resulting in many models each with a different combination of input features. Finally, the model with the optimal combination of features (⁠|x$ |⁠), |{c_1}\!..{c_x}$ |⁠, is selected, which is the model with the lowest cross validated P value. The coefficients from the Cox proportional hazard with the lowest cross validated P value as determined in the training set can be used to calculate a cell-type hazard score for each patient in the test set. The |H{S^{CT}}$ | is calculated as follows:

where |{\beta _1}$ | is the vector of coefficients for the main effects of the selected cell types, |C{T_p}$ | is the abundance of the selected cell types in patient |p$ |⁠, |\ {\beta _2}$ | if the coefficient for the main effect of treatment, and |{\beta _3}$ | is the vector of coefficients for the interaction effects between the selected cell types and treatment. Finally, the |DT{S^{CT}}$ | is calculated as follows:

where |HS_{p,\ T = 1}^{CT}$ | is the |H{S^{CT}}$ | for patient |p$ | when treated with treatment |{T_1}$ | (TAC) and |HS_{p,\ T = 0}^{CT}$ | is the |H{S^{CT}}$ | for patient |p$ | when treated with treatment T = 0 (dose-dense doxorubicin–cyclophosphamide; ddAC). The |DT{S^{CT}}$ | is therefore the log ratio in hazard for a patient when treated with TAC compared with when treated with ddAC. Predictive power of the |DT{S^{CT}}$ | was tested in a Cox proportional hazard model with the main effect of the |DT{S^{CT}}$ |⁠, the main effect of treatment, and the interaction between treatment and the |DT{S^{CT}}$ |⁠, corrected for clinicopathologic features tumor size, lymph node status, pathologic subtype, histologic grade, age, and type of surgery.

Both the CTS and the DTS are a continuum without clear cutoffs. To construct Kaplan–Meier graphs, cutoff points were chosen on the upper and lower quartile for CTS: 25% of patients with the lowest risk, 25% of patients with the highest risk, and 50% of patients with intermediate risk as determined by CTS. To construct Kaplan–Meier graphs for DTS, cutoff points were chosen on the upper and lower half of DTS: 50% of patients with the lowest risk and 50% of patients with the highest risk as determined by DTS.

For CTS, the patient cohort of the SCAN-B trial was split in three groups randomly with forced equal distribution of the clinical variables. Next, the first group was used for training and the last two groups were combined for validation. For DTS, the patient cohort of the MATADOR trial was split in two groups randomly with forced equal distribution of the clinical variables. Next, the first group was used for training and the last two groups were combined for validation.

The single-cell mRNA sequencing data are available through GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE168410. The bulk mRNA sequencing data is available through GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE167977. The private sensitive patient data of the MATADOR study requires controlled access. Access to this dataset will be managed locally by the Translational Review Board of the Netherlands Cancer Institute by contacting Mark Opdam ([email protected]). The code for Treatment Interaction Predictor is made available through CodeOcean (https://codeocean.com/capsule/6732361/tree/v2).

To optimally deconvolve the cellular composition of a tumor from bulk mRNA sequencing data through deconvolution algorithms, it is necessary to generate reliable reference profiles from the various cellular components of the tumor. To obtain these reference profiles, we randomly sorted living DAPI-negative cells from freshly isolated breast tumors of twelve patients using flow cytometry and performed single-cell sequencing (Supplementary Table S1). These patients were not part of the SCAN-B or MATADOR clinical trials that are used later in the manuscript for deconvolution. We additionally isolated, respectively, EpCAM-positive and CD14-positive cells to acquire enough cancer cells and macrophages. Next, we processed the samples for single-cell mRNA sequencing and characterized the expression profiles of in total 6,580 cells, including cancer and noncancer cells. In addition to single-cell mRNA sequencing, for 7 patients we collected enough material to also analyze the tumors by bulk mRNA sequencing and the cellular composition by pathological examination of hematoxylin and eosin (H&E) sections (see below).

To identify all the different cell types in the 6,580 cells that were sequenced, we performed cluster analysis using the Seurat R package and identified 15 clusters of cells with distinct expression profiles (Fig. 1A). The gene expression profiles showed that most of the clusters comprising cells of multiple patients represent the various non-transformed cell types present in the tumor microenvironment including cancer-associated fibroblasts (CAF), myoepithelial cells, muscle cells, endothelial cells, macrophages, T cells, B cells, and mast cells (Fig. 1B and C; Supplementary Fig. S2). Importantly, we did not observe coclustering of microenvironmental cells depending on treatment status of the patient (Fig. 1B; Supplementary Fig. S3A). These data suggest that if expression profiles change upon neoadjuvant treatment, these changes were not large enough to lead to separate clusters of neoadjuvant treated and untreated microenvironmental cells in our analyses (Fig. 1B). Although some epithelial cell (EC) clusters (I, II, and V) contain cells isolated from both chemotherapy-naïve and neoadjuvant-treated tumors, other EC clusters (III, IV, VI, and VII) only contain cells isolated from either naïve or treated tumors (Fig. 1B; Supplementary Fig. S3A). Importantly, because these EC clusters are derived from individual patients, differential expression between these clusters may therefore simply indicate variations between individuals.

Next, to identify which clusters can contain noncancerous epithelial cells, we isolated and sequenced cells from the contralateral non-diseased breast tissue from patient 7. The noncancerous epithelial cells of patient 7 were present in EC-I and EC-II, suggesting that these clusters contain a combination of transformed and nontransformed cells (Supplementary Fig. S3B). This was confirmed by the low variance of the copy-number variation (CNV) across the genome in these clusters (Supplementary Fig. S3C) compared with the higher variance of the CNVs in the other clusters (Supplementary Fig. S3C). Of the nontransformed cell cluster, EC-I had significantly increased expression of the hormone receptors ESR1, AR, and PGR compared with cluster EC-II (Supplementary Fig. S3D), suggesting a hormone receptor–positive luminal origin. In addition, we performed the PAM50 analysis on the single-cell data and we mapped the expression of KRT14, KRT5, ERBB2, and Ki67 (Supplementary Fig. S4). This analysis shows that individual cells, assigned to different subtypes by the PAM50 classification or specific marker genes, are distributed throughout the different epithelial cell types that we identified. Consequently, the various EC cell types do not represent different known subtypes, but rather a fine-tuned classification of intratumoral heterogeneity. Combined, these data show that single-cell mRNA sequencing of breast tumors allowed us to identify the expression profiles of the most prominent cell types that are present in breast tumors.

Based on the cluster analysis, we identified 15 different cell types in human breast tumors. Next, we averaged the expression profile of all cells within a cluster to generate reference profiles for each of the 15 identified cell types. We used these reference profiles to optimally detangle the contribution of each cell type from bulk mRNA sequencing data using TCD. To test the robustness of this approach, in the tumors of the 7 patients for which enough material was collected to perform both bulk mRNA sequencing and pathological examining of H&E sections, we estimated side-by-side the cellular composition of the tumors by TCD and by a well-trained pathologist. Importantly, we observed a strong correlation between H&E-based cell type scoring by a pathologist and cell type scoring by TCD [Pearson correlation > 0.93; 95% confidence interval (CI) 0.84–0.97; Supplementary Data; Supplementary Fig. S5A], showing that TCD can estimate the cellular composition of a tumor with similar precision as a pathologist. To further assess the strength of our approach, we compared our method (that uses DWLS; ref. 14) with other deconvolution algorithms [CIBERSORT (13) and MuSiC (30)]. Considering that important cell types are missing from the default reference data of these algorithms, we used our reference profiles and we applied them to CIBERSORT and MuSiC. We observed a stronger correlation between the pathologist score and TCD deconvolution score when using DWLS (Pearson correlation = 0.93) than when using CIBERSORT (Pearson correlation = 0.83) or MuSiC (Pearson correlation = 0.54), thereby illustrating the robustness of our method (Supplementary Fig. S5A–S5C).

Next, we questioned whether we could further improve the robustness of the TCD-estimated cellular composition by refining the cell type reference profiles if we would increase the size of the single-cell sequencing patient cohort. We redetermined the reference profiles of the different cell types based on iterative sampling of data of different numbers of patients from the cohort of 12 patients. We observed that the correlation does not improve significantly when the cell type reference profiles are based on single-cell sequencing data of a patient cohort larger than 7 (Supplementary Fig. S5D). We next also calculated the correlation per individual cell type that was scored by the pathologist, and again increasing the number of patients from 7 to 12 did not increase the correlation between pathologist scored cell type contribution and deconvolution-derived contribution (Supplementary Fig. S5E). Finally, we found in patients of the MATRADOR trial a strong correlation between pathology TIL scores and TCD TIL scores, confirming the robustness of our method in an independent set of patients (Supplementary Fig. S5F). Combined, these data indicate that TCD can reliably estimate the cellular composition of tumors.

Based on histopathologic and molecular features, breast cancer can be classified into different subtypes with distinct outcomes and response to therapies, including ER/PR positive and HER2 negative, HER2 positive, and triple-negative tumors (31–34). We used TCD to detangle the cellular composition of over 3,500 tumors of breast cancer patients with different subtypes and treatment histories that were enrolled in the MATADOR and SCAN-B clinical trials, and for which bulk mRNA sequencing data were available (Figs. 2 and 3). However, considering that the tumors used to generate our reference profiles did not include Her2⁺ patients, we excluded these patients from further analysis. From here on, the contribution of a cell type that is determined by TCD will be denoted as cell type^TCD.

Although the distinct pathologic subtypes have distinct outcome and response to therapies (31–34), hierarchical clustering of the cell type contribution matrix of patients in the MATADOR and SCAN-B trials lead to clustering of patients with the same pathological subtypes (Figs. 2 and 3). Importantly, each pathologic subtype was divided in two or more subclusters with specific cellular compositions. For example, in the MATADOR trial, we observe two clusters of triple-negative patients, one of which was enriched for recurrence events and showed high contribution from EC-V^TCD and low contributions from immune cells including macrophages^TCD, B cells^TCD, and T cells^TCD (Fig. 2A, indicated in the left orange dotted box). The other cluster of triple-negative patients was depleted of recurrence events and had low contribution of EC-V^TCD and high contribution of immune cells (Fig. 2A, indicated in the right blue dotted box). Kaplan–Meier analysis confirmed the differential RFS probability of triple-negative patients with high and low total immune infiltrates^TCD (sum of T cells, B cells, macrophages, and mast cells; Fig. 2B). The same two clusters can be found in the SCAN-B data (Fig. 3, blue and orange dotted boxes). Although luminal and triple-negative patients do not completely fall into separate clusters, we did identify an enrichment for triple-negative patients in two clusters with cellular compositions with high contribution of EC-V^TCD and low contribution of immune cells, and cellular compositions with low contribution of EC-V^TCD and high contribution of immune cells (Fig. 3A, indicated in the dotted boxes). Moreover, in line with the MATADOR data, patients in the SCAN-B trial with triple-negative breast tumors with high levels of immune cells have better survival than patients with triple-negative tumor that contain low levels of immune cells (Fig. 3B). These data are in line with previous notions that high immune infiltrate is correlated with better clinical outcomes (9, 11, 12, 35–37) illustrating the power of TCD to find relevant microenvironments. Our data also illustrate that in addition to immune cells, specific epithelial cell types are part of specific microenvironments conferring differential clinical outcomes (e.g., the abundance of EC-V).

Considering the increasing amount of evidence that the cellular composition of a tumor is important both in the development of metastases and in any benefit from treatment (1–6), we aimed to investigate whether differential RFS is characterized by different tumor compositions (i.e., cell type contributions as determined by TCD). For each patient in the SCAN-B trial, we calculated the cell type contributions by TCD. Next, we split the trial in a training set (first 892 patients) and a validation set (last 1,779 patients) with equal distribution of clinical variables (Supplementary Table S2). On the training set, we calculated a cell type score (CTS), a Cox regression score based on cell type contribution that is prognostic for the RFS. In our training set, we calculated the CTS in a 10-fold cross validation setting, where within each of the 10 folds, the optimal combination of cell types is selected using forward feature selection. The model with the lowest cross validation P value was selected (see Methods) and subsequently fixed and used in the independent SCAN-B validation set. We observed that the CTS of this fixed model was significantly prognostic in a Cox regression model for survival in the validation set (Likelihood ratio test P < 1^–12 and P = 0.003 corrected for clinicopathologic variables; Supplementary Table S3). We assigned three subgroups: 25% of patients with the lowest risk, 25% of patients with the highest risk, and 50% of the patients with intermediate risk. These subgroups indeed had clear distinct RFS both in the training and validation cohorts as shown in the Kaplan–Meier analyses (Fig. 4A and B; Supplementary Table S4). To validate these findings in a completely independent patient cohort, we calculated the CTS score for the same fixed model on the MATADOR patients and again observed that the CTS was significantly prognostic for survival in a Cox regression model (Likelihood ratio test P < 1^–6 and P < 0.001 corrected for clinicopathologic variables; Supplementary Table S5). By analyzing the same three subgroups in the independent MATADOR cohort, we could validate the differential RFS for the different subgroups (Fig. 4C; Supplementary Table S6). Next, we compared the risk scores determined by CTS with the PAM50 subtypes. In the SCAN-B (Supplementary Table S4), but not in the MATADOR cohort (Supplementary Table S6), we found an over-representation of the basal subtype in the high-risk group as well as an over-representation of the luminal A subtype in the low- and intermediate-risk groups. To test the prognostic ability of CTS within individual subtypes, we next split the patient cohorts in the SCAN-B and MATADOR trials into triple-negative and luminal groups. Importantly, in the cohorts with either triple-negative or luminal patients, we observed that CTS also successfully assigned the different risk groups (Supplementary Fig. S6).

Next, we investigated the differences in the cellular composition of the different subgroups with distinct risks (Fig. 4D and E). The CTS, and therefore the different subgroups, are mainly based on different contributions of EC-II^TCD cells, EC-V^TCD cells, EC-VI^TCD cells, macrophages^TCD, T cells^TCD, B cells^TCD, and myoepithelial cells^TCD (Fig. 4D and E). The low-risk group is enriched for EC-VI^TCD cells, T cells^TCD, and B cells^TCD (Fig. 4D and E). Our findings that some immune cells are enriched in the tumors of the low-risk group is in line with previous findings that high immune infiltrate is correlated with better clinical outcomes (23, 35, 36), again showing that our method robustly identifies relevant microenvironments. Moreover, in line with the preclinical findings that macrophages can have both tumor-supporting and tumor-suppressing phenotypes, we find macrophages^TCD to be present both in a high- and low-risk group. In addition to identifying immune cells, we also find differential abundance for epithelial cell types (EC-II, EC-V, and EC-VI). For example, the high-risk group has high contribution from EC-II^TCD and EC-V^TCD cells. Combined, we conclude that TCD of the SCAN-B trial enabled us to identify and subsequently validate specific microenvironments that distinguish differential survival of patients with breast cancer, and show the relevance of various types of immune cells for these environments, but also the relevance of various epithelial cell types.

Because preclinical and clinical data suggest that the microenvironment of a tumor is important for therapy benefit (1–6), we sought to investigate whether differential outcome after treatment is characterized by distinct cellular compositions as determined by TCD. To do this, we further analyzed the TCD results of the MATADOR trial. In this randomized controlled trial, patients received adjuvant treatment with either ddAC or TAC with a median follow-up of 7 years. We hypothesized that for every cellular composition, there is a certain ratio between the hazard if that patient would be treated with TAC and the hazard if that patient would be treated with ddAC. By subtracting the linear predictors from a Cox model, we obtain a DTS that is the log of this HR and directly shows which treatment provides the better outcome: DTS < 0 indicates better outcome when treated with ddAC and DTS > 0 indicates better outcome when treated with TAC. We split the patient cohort of the MATADOR trial in a training set of 254 patients and a validation set of 256 patients with equal distribution of clinical variables (Supplementary Table S7). To determine the DTS for all cellular compositions (DTS^CT), we constructed a Cox proportional hazard model based on the cell type contributions and their interaction with treatment on the training set (see Methods). The optimal combination of cell types is determined through a 10-fold cross validation of the training set of patients (see Methods). We then tested the resulting model in the validation set of patients. There was a significant interaction between the DTS^CT and treatment in a Cox proportional hazard model (P_interaction < 0.001 and P_interaction = 0.017 corrected for clinicopathologic variables; Fig. 5A; Supplementary Table S8). Next, we assigned two subgroups: 50% of the patients with the lowest DTS^CT and 50% of the patients with the highest DTS^CT (Supplementary Table S9). For all patients from the MATADOR trial that had an event, we plotted the DTS^CT in a forest plot (Fig. 5B). This plot demonstrates that within the patient that had an event in the group with high DTS^CT, there is an overrepresentation of patients that were treated with ddAC (Fig. 5B). This is in line with the Kaplan–Meier analysis that shows that patients with a high DTS^CT have a better survival probability if treated with TAC than with ddAC. This data shows that DTS can distinguish differential treatment benefits (Fig. 5C). The interaction between patients with high DTS^CT and patients with low DTS^CT and treatment was statistically significant in a Cox proportional hazard model corrected for clinicopathologic variables (P = 0.047; Supplementary Table S10).

Next, we analyzed the differential cellular composition between the patient groups with high and low DTS^CT. DTS^CT is mainly based on the contribution of EC-V, myoepithelial cells, and T cells. This analysis showed that high contributions of T cells^TCD and myoepithelial cells^TCD and a low contribution of EC-V^TCD result in higher DTS^CT and therefore represents the specific microenvironment for better outcome on TAC, while high contributions from EC-V^TCD, and a low contribution of myoepithelial cells^TCD and T cells^TCD result in lower DTS^CT and therefore represents a microenvironment for better outcome on ddAC (Fig. 5D and E). From this data, we conclude that differential treatment benefits for ddAC and TAC is characterized by specific microenvironments containing various abundances of immune and epithelial cell types.

Breast cancer is a heterogeneous disease both in terms of cellular composition and therapy benefit. We used deconvolution of bulk RNA sequencing data to unravel the cellular compositions of tumors of patients with differential survival and therapy benefit. Although tumor cellular composition can directly be determined with single-cell mRNA sequencing, this technique cannot be routinely performed for individual patients. Deconvolution methods are an alternative to obtain similar information on the cellular heterogeneity of a tumor from bulk mRNA sequencing data (13, 15). Importantly, if bulk mRNA sequencing is available from clinical trials, it can be used to determine the cellular composition without the need of extensive pathology examination. Moreover, deconvolution can estimate the contribution of cell types that cannot be distinguished by standard pathology (e.g., contribution of different ECs).

To obtain accurate deconvolution results, we constructed 15 reference profiles from single cells isolated directly from breast tumors. Our reference profiles also have some limitations. First, the reference profiles were constructed using samples from relatively young patients (<57 years) and did not include Her2⁺ patients. We have therefore excluded the Her2⁺ patient group in our analyses. Future single-cell mRNA sequencing studies may further improve our reference profiles.

By deconvolution of the bulk mRNA data from patients in the MATADOR and SCAN-B trials, we identified cellular compositions that characterize subgroups of patients with differential survival and treatment benefit. Although breast tumors can be classified to specific IHC subtypes, our approach allows for further refinement of these subtypes based on specific cellular compositions. We observed that within subtypes, patients’ tumors have different cellular compositions which are associated with different survival and therapy benefit. For example, TNBC patients’ tumors with a high abundance of immune cells, such as T cells^TCD, B cells^TCD, and macrophages^TCD, have a much better prognosis than their counterpart with tumors with low abundance of immune cells. This is in line with previous work showing that high immune infiltrate in breast and other types of cancer is correlated with better clinical outcomes upon treatments (9, 11, 12, 35). More precisely, we found that high infiltration of immune cells correlates with better response to taxane-containing chemotherapy. It is therefore possible that T cells, B cells, and macrophages, which are abundant in good responders to TAC, may drive the efficacy of taxane in vivo. In addition to immune cells, we have identified that survival and therapy benefit are characterized by various contributions of distinct epithelial cell types. Future work should establish the biological differences and functions of these different epithelial cell types.

In conclusion, our study shows that TCD is a powerful tool to determine the cellular compositions of tumors enabling the characterization of specific microenvironments that are associated with differential survival and treatment benefit. This promising method may allow further tailoring of treatment in the future, ultimately leading to higher survival rates and/or better quality of life.

L. Kester reports a patent for “Method for determining cellular composition of a tumor” (18020578.3) pending. E.H. Lips reports grants from CRUK/KWF outside the submitted work. A.E. van Leeuwen-Stok reports grants from Roche, Novartis, Philips, Servier, and Agendia outside the submitted work. M.E.M.M. Bos is on the advisory board at AstraZeneca on Trastuzumab-deruxtecan (Europe). H.M. Horlings reports other support from Roche Diagnostics BV during the conduct of the study. J. Wesseling reports research not related to the research presented in this manuscript as funded by Cancer Research UK, KWF Dutch Cancer Society ZonMW, and the Antoni van Leeuwenhoek Foundation. E.E. Voest reports grants from Roche, Pfizer, Clovis, Novartis, AstraZeneca, BMS, Eli Lilly, MSD, Amgen, Bayer, and Sanofi outside the submitted work. L.F.A. Wessels reports grants from Genmab BV during the conduct of the study. M. Kok reports grants from BMS, Roche, and AstraZeneca, as well as other support from Daiichi outside the submitted work. H.M. Oosterkamp reports grants from Amgen, Sanofi, and Dutch Cancer Society during the conduct of the study, as well as personal fees from Roche, MSD, Daiichi, AstraZeneca, Lilly, and Pfizer outside the submitted work. S.C. Linn reports grants from Amgen, Sanofi, and Dutch Cancer Society during the conduct of the study. S.C. Linn also reports grants and nonfinancial support from AstraZeneca, Genentech/Roche, Tesaro (now owned by GSK), and Immunomedics; grants from Eurocept-pharmaceuticals, Novartis, and Pfizer; and other support from Cergentis, IBM, Daiichi Sankyo, and Bayer outside the submitted work. J. van Rheenen reports a patent for “Method for determining cellular composition of a tumor” (18020578.3) pending. No disclosures were reported by the other authors.

L. Kester: Conceptualization, software, formal analysis, validation, methodology, writing–original draft, writing–review and editing. D. Seinstra: Investigation, writing–original draft. A.G.J. van Rossum: Resources, validation, investigation, methodology. C. Vennin: Data curation, investigation, writing–review and editing. M. Hoogstraat: Resources, data curation, formal analysis. D. van der Velden: Resources. M. Opdam: Resources, data curation, formal analysis, validation. E. van Werkhoven: Data curation, methodology. K. Hahn: Data curation, supervision, investigation. I. Nederlof: Resources, data curation, validation. E.H. Lips: Resources, data curation, supervision. I.A.M. Mandjes: Resources. A.E. van Leeuwen-Stok: Resources. S. Canisius: Data curation, validation, methodology. H. van Tinteren: Resources. A.L.T. Imholz: Resources. J.E.A. Portielje: Resources. M.E.M.M. Bos: Resources. S.D. Bakker: Resources. E.J. Rutgers: Resources. H.M. Horlings: Resources, data curation, supervision. J. Wesseling: Resources, data curation, supervision. E.E. Voest: Resources, data curation, supervision, investigation. L.F.A. Wessels: Data curation, software, supervision, investigation, methodology. M. Kok: Resources, supervision. H.M. Oosterkamp: Conceptualization, resources, supervision, funding acquisition. A. van Oudenaarden: Conceptualization, data curation, supervision. S.C. Linn: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, writing–review and editing. J. van Rheenen: Conceptualization, supervision, funding acquisition, writing–original draft, writing–review and editing.

First, we would like to thank all patients and their families that were included in this study. Also, we would like to express our gratitude to the study teams of the participating centers, the Data Center of the Netherlands Cancer Institute for collecting the data, the Core Facility – Molecular Pathology and Biobanking of the Netherlands Cancer Institute for their help with the immunohistochemistry stainings, the Genomics Core Facility of the Netherlands Cancer Institute for their work on the RNA sequencing, the Hubrecht and NKI FACS facility, and the Dutch Breast Cancer Research Group (BOOG) for their role in coordinating the study. This work was supported by CancerGenomics.nl (Netherlands Organization for Scientific Research) program (to J. van Rheenen) and the Josef Steiner Cancer Research Foundation (to J. van Rheenen). The MATADOR trial was funded by the Dutch Cancer Society (CKTO 2004–04) and by unrestricted institutional research grants from Sanofi and Amgen (to S.C. Linn and H.M. Oosterkamp). C. Vennin was supported by a fellowship from the Human Frontiers in Science Program (LT000082/2018-L3).

The MATADOR trialists’ group is listed as follows: S.C. Linn: The Netherlands Cancer Institute, Amsterdam, The Netherlands; M. Soesan: Netherlands Cancer Institute, Amsterdam, The Netherlands; H.M. Oosterkamp: Haaglanden Medical Center, The Hague, The Netherlands; Dr. F. Jeurissen, Haaglanden Medical Center, The Hague, The Netherlands; N. Weijl, Haaglanden Medical Center, The Hague, The Netherlands; A.L.T. Imholz: Deventer Ziekenhuis, Deventer, The Netherlands; J.E.A. Portielje: HagaZiekenhuis, The Hague, The Netherlands; K.J. Beelen: Reinier de Graaf Gasthuis, Delft, The Netherlands; S. Bakker: Zaans Medisch Centrum, Zaandam, The Netherlands; G. de Klerk: Waterland Ziekenhuis, Purmerend, The Netherlands; S. Vrijaldenhoven: Noordwest Ziekenhuisgroep, Alkmaar, The Netherlands; A. van der Velden: Martini Ziekenhuis, Groningen, The Netherlands; H. de Graaf: Medisch Centrum Leeuwarden, Leeuwarden, The Netherlands; M. Smeets: Streekziekenhuis Koningin Beatrix, Winterswijk, The Netherlands; J. Meerum Terwogt: Onze Lieve Vrouwe Gasthuis, Amsterdam, The Netherlands; J. Schrama: Spaarne Gasthuis, Hoofddorp, The Netherlands; P. Kuijer: Spaarne Gasthuis, Hoofddorp, The Netherlands; H. Wilmink: Amsterdam University Medical Center, Location AMC, Amsterdam, The Netherlands; R. Hoekstra: Ziekenhuisgroep Twente, Hengelo, The Netherlands; J. Kroep: Leiden University Medical Center, Leiden, the Netherlands; J.F.M. Pruijt, Jeroen Bosch Ziekenhuis, Den Bosch, The Netherlands; L. van Gerven, LangeLand Ziekenhuis, Zoetermeer, The Netherlands; A.H. Vos, Bernhoven, Oss, The Netherlands; F. Erdkamp, Zuyderland Ziekenhuis, Sittard, The Netherlands; W.G. van Leeuwen-Breuk, BovenIJ Ziekenhuis, Amsterdam, The Netherlands; A. de Graeff, University Medical Center Utrecht, The Netherlands.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.