Triple-negative breast cancer (TNBC) is a molecularly heterogeneous cancer that is difficult to treat. Despite the role it may play in tumor progression and response to therapy, microenvironmental (stromal) heterogeneity in TNBC has not been well characterized. To address this challenge, we investigated the transcriptome of tumor-associated stroma isolated from TNBC (n = 57). We identified four stromal axes enriched for T cells (T), B cells (B), epithelial markers (E), or desmoplasia (D). Our analysis method (STROMA4) assigns a score along each stromal axis for each patient and then combined the axis scores to subtype patients. Analysis of these subtypes revealed that prognostic capacity of the B, T, and E scores was governed by the D score. When compared with a previously published TNBC subtyping scheme, the STROMA4 method better captured tumor heterogeneity and predicted patient benefit from therapy with increased sensitivity. This approach produces a simple ontology that captures TNBC heterogeneity and informs how tumor-associated properties interact to affect prognosis. Cancer Res; 77(17); 4673–83. ©2017 AACR.

In clinical practice, breast cancer is stratified into subtypes defined by differential protein expression of the estrogen receptor (ER) and progesterone receptor (PR), as well as expression and/or genomic amplification of human epidermal growth factor receptor 2 (HER2). Together with other clinical variables including patient age, tumor size, grade, and lymph node status, these factors form the basis for clinical decision-making with respect to treatment (1). The 10%–15% of breast carcinomas that lack expression/amplification of ER, PR, or HER2 form a subtype that is termed triple-negative breast cancer (TNBC; ref. 1). Tumors of the TNBC subtype are associated with earlier age of onset, higher grade at presentation, and overall poorer patient prognosis (2). Despite extensive efforts, this subtype still lacks targeted therapies and these tumors are generally treated with untargeted chemotherapy and radiation (3, 4).

Previous studies, including several high-throughput profiling efforts, have indicated that the TNBC subtype has higher levels of intertumoral (patient-to-patient) heterogeneity when compared with other subtypes with respect to both gene expression (5), and somatic genomic aberrations (6, 7). This heterogeneity may at least partially underlie why TNBC is a poor outcome subtype (8–11). Several efforts have investigated whether there are subtypes within TNBC with distinct cellular processes and responses (12–14). These studies, however, have used gene expression profiling of bulk samples enriched for epithelial cells of the tumor proper. Lehmann and colleagues (13) identified six subtypes within TNBC patients: BL1 (basal-like 1), BL2 (basal-like 2), M (mesenchymal), MSL (mesenchymal stem-like), IM (immunomodulatory), and LAR (luminal androgen receptor) – and demonstrated that these “TNBCType” subtypes are associated with differential responses to neoadjuvant chemotherapy (15). This scheme has since been translated into an assay for clinical implementation (INSIGHT TNBCTYPE; ref. 16).

It is well-established that the microenvironment of the tumor coevolves with the tumor proper, and that this coevolution is central in tumor progression, and ultimately patient outcome (17). To date the few studies that have investigated stromal gene expression profiles have utilized a pan-breast cancer cohort (18–20). Our effort here is the first to investigate stromal heterogeneity specifically within TNBC using microdissected stromal compartments of both the tumor proper and adjacent morphologically normal cells from 57 TNBC clinical samples. An unbiased class discovery approach identified four stromal axes with distinct molecular and clinical qualities associated with the TNBC microenvironment. Combining the stromal axis scores to generate a novel subtyping approach revealed interactions between the stromal axes with distinct prognostic associations.

Sample selection

Samples were collected from patients undergoing breast surgeries at the McGill University Health Centre (MUHC) between 1999 and 2012 who provided written, informed consent (MUHC REB protocols SDR-99-780 and SDR-00-966). All tissues were snap-frozen in O.C.T. Tissue-Tek Compound within 30 minutes of removal. Information regarding clinical variables was obtained through review of medical records. Samples used for this cohort (n = 57) were reported as negative for ER (based on a 10% threshold) and HER2 via IHC (ER and HER2) and/or FISH (HER2). All patients were PR-negative (IHC), with the exception of one case with weak expression. Hematoxylin and eosin-stained sections from each sample were evaluated by an attending clinical pathologist with expertise in breast tissue to identify representative areas of tumor and tumor-associated stroma, as well as histologically normal breast epithelium and stroma. Laser capture microdissected RNA was isolated from samples and hybridized to Agilent Technologies SurePrint G3 Human GE 8 × 60K Microarrays (See Supplementary Information).

Microarray dataset normalization

R/Bioconductor (vers 3.20; Bioconductor 3.1; ref. 21) was used for most analyses. Normalization was performed using the limma package (22) where loess was applied for dye bias correction, and quantile normalization was used across arrays. Replicates of noncontrol probes were aggregated by taking their mean value. To investigate technical error introduced in the LCM/microarray procedure, the stromal and matched adjacent normal sample from a single patient were repeated and found to be highly concordant (>0.8). Replicate expression profiles were then averaged for the remainder of the analysis. The most variable probe was chosen when there were multiple probes for the same transcript. Clear distinctions between normal and tumor stromal samples were also observed (See Supplementary Information; Supplementary Fig. S1). Raw and normalized microarray data have been deposited in the Gene Expression Omnibus database under accession number GSE90505.

Discovery of stromal axes

Probes with an interquartile range > 2.0 across all samples were used as features in hierarchical clustering (Ward's algorithm, Pearson correlation distance). Samples were mean-centered and scaled for each transcript across all patients prior to clustering. Pvclust (version 1.3-2; ref. 23), was used to measure cluster stability with 100,000 iterations and we selected clusters containing at least 12 genes with an Approximately Unbiased (AU) value of >85%. Gene clusters deemed unstable by Pvclust were not selected for further analysis, and no further characterization was performed before selection.

Linear order and assignment of signatures using ROI95

This unbiased approach, which is described in ref. 24, ranks samples based on their expression of a specific gene set. In our case, we use the characteristic genes for each stromal and Lehmann and colleagues axis. This method estimates each patient sample as either low, intermediate, or high using a random resampling technique with 1,000,000 iterations.

Differential gene expression and pathway analysis

To identify differentially expressed (DE) genes for each stromal axis, we fitted a linear model comparing the levels for each stromal axis using the R package limma (22) and corrected with Benjamini–Hochberg (P < 0.05). DE gene lists were examined using QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN; https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/) and compared against the Molecular Signatures Database (MSigDB v5.2; http://software.broadinstitute.org/gsea/msigdb) for pathway analysis.

Assignment of other subtyping schemes to patients across large patient cohort

We used our compendium of 5,901 bulk expression profiles of invasive breast cancer samples from 13 nonoverlapping datasets generated on different technologies (5). We define poor outcome as an observed distant metastasis within 5 years of diagnosis (where available) and used ER and HER2 status as reported for each dataset where available. As many datasets lacked information on PR status, we used ER and HER2 negativity to define TNBC patients. All patients within the compendium were also labeled with intrinsic subtype values via PAM50 (25), and TNBC patients were labeled by TNBCType via the web-based tool (26).

Generation of stromal and Lehmann axis scores in TNBC whole tumor samples

To assign scores along the stromal axes to TNBC patients across the compendium of bulk expression profiles, we performed ROI95 using the lists of differentially expressed genes for each stromal axis. Our software to generate stromal axis scores in a sample is available as a Bioconductor package entitled STROMA4 (http://bioconductor.org/packages/devel/bioc/html/STROMA4.html). STROMA4 was applied to each (of the 13) datasets of our TNBC compendium independently and the three classes (low, intermediate, high) were combined across all of these datasets. The Lehmann axes were similarly assigned using the characteristic genelist for each Lehmann subtype from the original publication. Assignments by ROI95 were compared to assignments by the TNBCType web-based tool (26). Only four disagreements for “high” classification by both tools were observed (Supplementary Table S1) confirming the accuracy of the ROI95 assignments.

Statistical analysis

Cohen kappa statistic was used to measure agreement between two ROI95-based categorizations into low, intermediate, and high (fmsb package vers. 0.5.1). Enrichment analyses were performed via a one-sided Fisher exact test in R. The minimum P value between each one-sided test was used to determine whether there was significant enrichment or depletion. For variables with more than two levels, multiple tests were performed to determine enrichment for each level individually against all other levels.

For comparisons between a ternary variable and a binary variable (for example, lymph node status versus high, intermediate, low stromal axis score), we removed the intermediate category and then used Cohen kappa for the two binary variables. To determine association with distant metastasis free-survival, we used a Cox proportional hazards regression model via the coxph function in R (vers. 2.38) using the three (ordinal) levels estimated by the ROI95.

Testing for enrichment of combinations of stromal axis scores

To determine whether the fraction of observed combinations of stromal axis scores was higher than expected, we used a one-sided binomial test (stats package in R; ref. 21). The observed number of patients with the combination being tested was used for the number of successes, the total number of patients was used for the number of trials, and the hypothesized probability of success was determined as the product of the fractions for the individual axis scores being tested as observed in the TNBC patient compendium. A P value less than 0.05 was deemed to be significant.

Expression profiling of microdissected tissue reveals four stromal axes in TNBC

To investigate stromal heterogeneity across TNBC tumors, 57 patient samples were selected based on negative ER, PR, and HER2 status according to clinical-pathological reports (Supplementary Table S2). Tumor nonepithelial (stromal) compartments were separately isolated by laser capture microdissection (LCM) and subjected to microarray-based gene expression profiling (see Supplementary Methods; refs. 18, 19, 27). Matched histologically normal stromal tissue was isolated for a subset of cases (n = 11). Following quality control and data normalization (see Materials and Methods), the success of the LCM procedure in isolating distinct tissue compartments was confirmed (see Supplementary Methods).

To establish whether differences are observed in TNBC stroma, the most variable genes in TNBC tumor stromal samples [interquartile range (IQR) > 2, n = 211 genes] were subjected to hierarchical clustering (Ward algorithm, Pearson correlation distance). Four distinct clusters were observed that contained a significant number of genes with strong pairwise gene-gene correlations of expression (Fig. 1A, colors along rows). These clusters were statistically stable and reproducible (pvclust, AU ≥ 85%). Genes within each cluster that exhibit strong coexpression across the patient cohort are termed the characteristic gene set.

Figure 1.

Hierarchical clustering identifies four stromal gene clusters in TNBC samples associated with disease outcome. A, Hierarchical clustering of tumor stromal gene expression profiles using genes with IQR > 2. Rows, transcripts; columns, samples. Stable clusters with A.U. > 0.85 and > 12 genes are indicated by colored bars at left. Values are centered and scaled per transcript across all samples and represented by the color key. PAM50 assignments and clinical outcome for each sample are represented under the heatmap. B, Assignment of samples into three classes [high, intermediate (ROI), or low] for each axis using ROI95 (classes demarcated by dashed lines in heatmaps). Patients with the smallest sum of expression (score) are ranked lowest and depicted in lightest color (at right) and those with the largest score are ranked highest and depicted in the darkest color (at left). Vertical colored bars at left of each heatmap correspond with the color assigned to samples high for that subtype. Rows, transcripts; columns, samples. Values are centered and scaled per transcript across all samples and are represented by the color key. C, Relationships between the assignments for each stromal axis. Rows, transcripts; columns, patients as in A above. Patient rankings for each cluster are denoted by colors as in B. Note that samples can be high for multiple stromal axes. D, Concordance between assignments of matched patients in LCM stroma and bulk expression datasets. Columns, patient assignments from LCM stroma profiles; bar color density, assignments to low (light), intermediate (medium), or high (dark) from bulk expression profiles. Colors denote clusters as in panels A, B, and C above. E, Kaplan–Meier survival analysis of the stromal axis scores for distant metastasis-free survival of TNBC patients in external TNBC bulk expression datasets (n = 1,098). Log-rank test P values are indicated at bottom left for each graph.

Figure 1.

Hierarchical clustering identifies four stromal gene clusters in TNBC samples associated with disease outcome. A, Hierarchical clustering of tumor stromal gene expression profiles using genes with IQR > 2. Rows, transcripts; columns, samples. Stable clusters with A.U. > 0.85 and > 12 genes are indicated by colored bars at left. Values are centered and scaled per transcript across all samples and represented by the color key. PAM50 assignments and clinical outcome for each sample are represented under the heatmap. B, Assignment of samples into three classes [high, intermediate (ROI), or low] for each axis using ROI95 (classes demarcated by dashed lines in heatmaps). Patients with the smallest sum of expression (score) are ranked lowest and depicted in lightest color (at right) and those with the largest score are ranked highest and depicted in the darkest color (at left). Vertical colored bars at left of each heatmap correspond with the color assigned to samples high for that subtype. Rows, transcripts; columns, samples. Values are centered and scaled per transcript across all samples and are represented by the color key. C, Relationships between the assignments for each stromal axis. Rows, transcripts; columns, patients as in A above. Patient rankings for each cluster are denoted by colors as in B. Note that samples can be high for multiple stromal axes. D, Concordance between assignments of matched patients in LCM stroma and bulk expression datasets. Columns, patient assignments from LCM stroma profiles; bar color density, assignments to low (light), intermediate (medium), or high (dark) from bulk expression profiles. Colors denote clusters as in panels A, B, and C above. E, Kaplan–Meier survival analysis of the stromal axis scores for distant metastasis-free survival of TNBC patients in external TNBC bulk expression datasets (n = 1,098). Log-rank test P values are indicated at bottom left for each graph.

Close modal

To assign values along each stromal axis to TNBC tumors, patients were linearly ordered on the basis of the sum of ranks of the characteristic genes for each stromal axis signature independently (see Materials and Methods). A rank-based permutation test (ROI95; ref. 24), was applied to each linear ordering to estimate boundaries of regions that delineate samples that are low, intermediate, or high for the characteristic gene set (Fig. 1B, black bars). Hence, each patient sample is independently scored for each of the four ternary axes (low, medium, high). This approach differs from traditional subtyping approaches that partition the patient cohort into distinct, nonoverlapping subtypes (Fig. 1C).

Each stromal axis is associated with markers of distinct cell types and processes

To characterize the molecular pathways and presence of specific cell types in each stromal axis, we identified differentially expressed genes between patients deemed low versus those deemed high for each stromal axis (LIMMA, FDR-adjusted P < 0.05 after ROI95, Supplementary Fig. S2; Supplementary Table S3). For the first axis, genes differentially expressed include both general (CD2, CD3D, IL2Rα, IL2Rβ, IL2Rγ), and lineage-specific (CD4, CD8A, CD8B) T-cell–associated markers, as well as markers of a Th1-mediated antitumor response including IL15 (28), granzymes (GZMBA, GZMB, GZMK, GZMH; ref. 29), markers of an IFN response (IFI30, IFIT5; ref. 30), transcription factors involved in Th1 differentiation (STAT1, STAT4; refs. 31, 32), and TNFα-induced genes (TNFAIP2, TNFAIP8; refs. 33, 34). These genes had greatest expression in patients with high scores for this axis (purple, Fig. 1B and C). Pathway analysis for the first axis identified signatures linked to the proliferation of T lymphocytes and activation of cytotoxic T cells, confirming the associations with T cells (purple, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4).

For the second axis, genes differentially expressed between low and high classes (magenta, Fig. 1B and C) include B-cell markers (CD19, CD79A, CD72), immunoglobulins (IGLL5, IGLL1, IGJ), and transcription factors associated with B-cell activation (POU2AF1, XBP1). These genes have greatest expression in samples scoring as high along this axis. Pathway analysis for the second axis identified signatures of B-cell proliferation (magenta, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4).

For the third axis, differentially expressed genes (teal, Fig. 1B and C) include keratins (KRT6B and KRT23), and metallothioneins. These genes are expressed by tumor epithelial cells (35) and thus may represent invasive tumor cells that have retained some of their epithelial characteristics due to tumor plasticity (36).

For the fourth axis, differentially expressed genes (orange, Fig. 1B and C) include multiple collagens (collagens 1A1/2, 3A1, 5A1/2, 8A1/2, 10A1, 12A1, 16A1), PDGFRB, FAP, in addition to collagen stabilizing and modifying enzymes (P4HA2, MMP2, LOXL1). All of these are factors associated with a desmoplastic reaction (37, 38). Analysis of the fourth axis identified a signature of desmoid-type fibromatosis (hypergeometric test, P < 0.05; orange, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4; ref. 39).

On the basis of these observations, we labeled the four stromal axes as T (T-cell, purple), B (B-cell, magenta), E (invasive epithelial cells, teal), and D (desmoplastic reaction, orange), respectively.

Stromal axes are associated with outcome in bulk expression profiles

To test associations between the stromal axis scores and clinical variables (e.g., tumor size, grade, stage, outcome, etc.) requires a larger cohort of TNBC patient samples. Because of the unavailability of TNBC stromal datasets, we developed and tested a statistical method to estimate the status along each stromal axis in bulk expression data (Fig. 1D; Supplementary Table S5). This method, entitled STROMA4, was applied to a large cohort of TNBC patient samples (n = 1,098) selected from 13 individual nonoverlapping publicly available breast cancer datasets (see Materials and Methods; Supplementary Methods; ref. 5). Stromal axis assignments were computed independently per dataset, and pooled across the constituent datasets (Supplementary Table S6). This enabled us to test if the low, intermediate, and high partitions of each axis stratified patients by clinical outcome. While the D score (orange) did not demonstrate significant association with outcome, the T, B, and E scores (purple, magenta, teal) were significantly correlated with outcome [log-rank test, distant metastasis free survival (DMFS) at 5 years all P < 0.05; Fig. 1E]. This demonstrates that scores along the T, B, and E stromal axes inform on clinical outcome for TNBC patients.

Stromal axes succinctly summarize TNBC heterogeneity

To establish whether the stromal axis scores are associated with subtypes of TNBC derived from bulk tumor gene expression profiles (TNBCType; ref. 13), we first applied our approach to the TNBCType subtypes. Using the data from Lehmann and colleagues (13), the gene sets that underlie each of the six TNBC subtypes were subjected to our methodology, estimating their activation as either low, intermediate, or high across the TNBC compendium (Materials and Methods). This procedure generates six “Lehmann axes” in a format that allows comparison with the four stromal axes.

To examine if any of the Lehmann axes interact (e.g., whether samples high for one subtype are also high for the second), associations between all possible states of all possible pairs of Lehmann axis scores were determined. Subtyping schemes partition patient samples into disjoint groups with the assumption that if a patient belongs to one subtype (for which they have the appropriate molecular profile), they do not belong to another subtypes (their profile is sufficiently distinct). However, we observe strong statistical correlations between 11 out of 15 pairs of TNBCType axes (Fig. 2A, Kappa test, all P values <0.01). While the mesenchymal (M) and immunomodulatory (IM) scores display a near perfect (P < 1e−10) anticorrelation that is consistent with distinct subtypes, we also observe an equally strong (anti-)correlation between the basal-like-2 (BL2), luminal androgen receptor (LAR), and mesenchymal stem-like (MSL) scores. Hence, for example, even though by our method many patients are high for both BL2 and LAR, each will be assigned to only one of these subtypes under the TNBCType approach. This correlation does not support their identification as distinct subtypes.

Figure 2.

Relationships between the axis scores for TNBCType and our stromal axis scores. Heatmaps depicted summarize ROI95 assignments for each TNBCType axis across the TNBC compendium. Samples are colored white, gray, and black to represent low, intermediate, and high subtype assignments respectively. Stromal axes are colored as before. Representative graphs are also drawn to summarize the observed associations. A, Agreement between the Lehmann axes as determined by the Kappa test. B, Agreement between the stromal axes as determined by the Kappa test. C, Patients are ordered by the T axis score. D, Patients are ordered by the B axis score. E, Patients are ordered by the E axis score. F, Patients are ordered by the D axis score.

Figure 2.

Relationships between the axis scores for TNBCType and our stromal axis scores. Heatmaps depicted summarize ROI95 assignments for each TNBCType axis across the TNBC compendium. Samples are colored white, gray, and black to represent low, intermediate, and high subtype assignments respectively. Stromal axes are colored as before. Representative graphs are also drawn to summarize the observed associations. A, Agreement between the Lehmann axes as determined by the Kappa test. B, Agreement between the stromal axes as determined by the Kappa test. C, Patients are ordered by the T axis score. D, Patients are ordered by the B axis score. E, Patients are ordered by the E axis score. F, Patients are ordered by the D axis score.

Close modal

To establish whether the four stromal axis scores show similar associations, we repeated this analysis with the stromal axes. Of the 6 pairs tested, only the T and B stromal axis scores were strongly correlated (P < 1e−10), while the D score showed a weak anticorrelation (P < 0.01) with the T and E scores (Fig. 2B).

To investigate associations between the Lehmann axis scores and our TNBC stromal axis scores, we performed a similar analysis comparing the Lehmann and stromal scores (Kappa test, Fig. 2C–F). Notably, the T stromal axis score (Fig. 2C), and to a lesser extent the B score (Fig. 2D), capture the inversely-correlated M and IM Lehmann axis scores (P < 1e−10). This observation is in line with the observation by Lehmann and colleagues that the IM score is correlated with lymphocytic infiltration (40). The stromal E score exhibits strong correlation with the BL1 score and anticorrelation with the LAR score (P < 1e−10; Fig. 2E). Patient samples estimated high for the D axis are almost always estimated high for the BL2, LAR, and MSL axes and low for the BL1 axis (P < 1e−10, Fig. 2F), which is in line with the observation that the MSL subtype is correlated to the presence of CAFs in the tumor (40). Together these observations highlight that the “Lehmann axes” are strongly associated with the stromal axes, and suggests that TNBC heterogeneity can be succinctly summarized by three distinct scores related to immune infiltration (B and T), androgen receptor signaling (inversely related to E), and a desmoplastic stroma (D).

The D axis is the stromal image of tumor proliferation

High expression of the MSL and BL1 Lehmann axes was observed to be negatively and positively correlated with the proliferative index of the tumor, respectively (13), and in the preceding subsection we show that our D stromal axis is also strongly positively and negatively associated with the MSL and BL1 Lehmann axes, respectively. This suggests that the D axis score reflects the stromal response to the proliferation of neighboring tumor epithelial cells.

We investigated this relationship in three ways. First, samples with high scores along the D axis had a significantly lower percentage of Ki67-positive tumor cells than samples with low D scores (two-sided t test, P < 0.05, Fig. 3A). Second, using a gene signature of proliferation (41) and the ROI95 to estimate proliferative states, we observed a strong statistical association between scoring along the D axis and expression of the proliferative signature (Fig. 3B and C; Kappa test, P < 0.01). Third, we observed that samples scoring low along the D axis are associated with higher grade (Fig. 3D). High-grade tumors tend to have a higher mitotic index (42). Together, these findings indicate that the D axis score is the stromal image of tumor proliferation from the neighboring epithelial cells.

Figure 3.

Increased levels of the D axis are associated with decreased proliferation and lower grade. A, Boxplot showing the association of the D axis score with Ki-67 staining as a marker for proliferation. *, comparisons that are significantly different (two-sided two-sample t test, P < 0.05). B, Ordering of patient-matched bulk expression dataset using a proliferation signature indicates an association with D status (via ROI95). C, Barplots indicating the fraction of high, intermediate, and low proliferation patients that are high, intermediate, or low for the D axis. Patients were assigned to classes by ROI95 using the proliferation signature (horizontal axis; B) and the characteristic gene set for the D axis (vertical axis). D, Barplots indicating the fractions of patients with low, intermediate, or high D axis scores with respect to grade.

Figure 3.

Increased levels of the D axis are associated with decreased proliferation and lower grade. A, Boxplot showing the association of the D axis score with Ki-67 staining as a marker for proliferation. *, comparisons that are significantly different (two-sided two-sample t test, P < 0.05). B, Ordering of patient-matched bulk expression dataset using a proliferation signature indicates an association with D status (via ROI95). C, Barplots indicating the fraction of high, intermediate, and low proliferation patients that are high, intermediate, or low for the D axis. Patients were assigned to classes by ROI95 using the proliferation signature (horizontal axis; B) and the characteristic gene set for the D axis (vertical axis). D, Barplots indicating the fractions of patients with low, intermediate, or high D axis scores with respect to grade.

Close modal

Stromal axis interactions induce 15 enriched subtypes with larger than expected populations

The many alternative subtyping schemes for breast cancer propose varying numbers of patient partitions that range from just four subtypes via classic approaches based on ER, PR and HER2 status to 10 subtypes identified by IntClust through joint DNA and mRNA analysis (43). Similarly, there have been different numbers of “sub-sub”-types proposed for TNBC (13, 14). Our four ternary stromal axes suggests that as many as 34(=81) such sub-subtypes could exist; each such subtype is described as a combination of the four axis scores (eg T-high, B-high, E-low, D-int). However, using our compendium of approximately 1,000 TNBC profiles, we observed that the number of patients assigned to each of the 81 possible combination subtypes varied significantly, with some subtypes populated by many samples and others with very few.

To measure this statistically, the enrichment or depletion of subtype populations were assessed relative to a background model (binomial-based test, see Materials and Methods). Only 15 of the 81 subtypes are significantly enriched beyond levels that would be observed solely by chance across our compendium (Supplementary Table S7). Many of the 15 enriched subtypes were either of the form T-high, B-high/B-intermediate, or of the form T-low, B-low/B-intermediate, suggesting a strong interaction between the T and B axis scores. Conversely, 17 of the 81 subtypes were significantly depleted, suggesting selection against some specific combinations of stromal axis scores (for example, the T-high, B-low/intermediate, D-high, E-high combination), in turn suggesting a complex interaction between these axes. The remaining 49 subtypes were populated as would be expected under a null binomial model.

The D axis is a master controller of the prognostic role of the T, B, and E axes

Above, we have investigated the prognostic capacity (DMFS at 5 years) of the stromal axis scores individually, of which the B, T, and E axes were statistically significant. To establish whether the interactions between the stromal axes have prognostic capacity in the TNBC compendium, we focused on the 15 stromal subtypes that had larger than expected populations. We observe that all subtypes where the D score is high were not significantly different in terms of the prognostic capacity (Fig. 4A, log-rank test, P > 0.05). However, when the D score is low or intermediate, the spectrum of subtypes (varying states of B, T and E) do have significantly different survival characteristics (Fig. 4B; log-rank test, P < 0.05). Additional survival analysis across the full TNBC patient cohort supported these findings (Supplementary Methods; Supplementary Fig. S3).

Figure 4.

The D axis score controls the prognostic effect of the B, T, and E axes. A, Kaplan–Meier curves showing lack of prognostic value of B, T, or E axis scores in D-high patients. Only the D-high subset of the 15 overrepresented combined classes is shown. ↑, |, and ↓ represent high, intermediate, and low scores for stromal axes, respectively. B, Kaplan–Meier curves showing prognostic value of B, T, and E axis scores in D-intermediate and D-low patients. Only the D-intermediate and D-low subsets of the 15 overrepresented combined classes are shown. ↑, |, and ↓ represent high, intermediate, and low assignments for stromal axes, respectively. C, Boxplot showing the association of the D axis score with difficulty to predict poor prognosis but not good prognosis. D, Heatmap depicting the subtype-specific performance of prognostic classifiers. Colors are proportional to the rank of the classifier within the specific patient cohort, with red representing the highest-performing classifiers relative to the remaining classifiers. Ticks represent the level of significance of the classifier (log-rank test, P < 0.05, 0.01, 0.001, respectively). Stromal axis classifiers (predictors) and the closest adjacent signatures are highlighted. E, Subtyping schema for TNBC patients based on observed associations with patient prognosis. The “low” category for each axis also includes patients scored as intermediate for that axis. The number of patients in the TNBC compendium assigned to each branch/leaf of the tree is indicated.

Figure 4.

The D axis score controls the prognostic effect of the B, T, and E axes. A, Kaplan–Meier curves showing lack of prognostic value of B, T, or E axis scores in D-high patients. Only the D-high subset of the 15 overrepresented combined classes is shown. ↑, |, and ↓ represent high, intermediate, and low scores for stromal axes, respectively. B, Kaplan–Meier curves showing prognostic value of B, T, and E axis scores in D-intermediate and D-low patients. Only the D-intermediate and D-low subsets of the 15 overrepresented combined classes are shown. ↑, |, and ↓ represent high, intermediate, and low assignments for stromal axes, respectively. C, Boxplot showing the association of the D axis score with difficulty to predict poor prognosis but not good prognosis. D, Heatmap depicting the subtype-specific performance of prognostic classifiers. Colors are proportional to the rank of the classifier within the specific patient cohort, with red representing the highest-performing classifiers relative to the remaining classifiers. Ticks represent the level of significance of the classifier (log-rank test, P < 0.05, 0.01, 0.001, respectively). Stromal axis classifiers (predictors) and the closest adjacent signatures are highlighted. E, Subtyping schema for TNBC patients based on observed associations with patient prognosis. The “low” category for each axis also includes patients scored as intermediate for that axis. The number of patients in the TNBC compendium assigned to each branch/leaf of the tree is indicated.

Close modal

To determine whether this effect was dependent on chemotherapy treatment, patients were stratified on the basis of whether they received chemotherapy and survival analysis was performed (Supplementary Figs. S4 and S5). As observed in the unstratified cohort, the D axis score acts as a master regulator of the other axes; none of the B, T, or E axes are prognostic in either chemotherapy-treated or -untreated D-high patients. However the B and T axes are prognostic among D-low/intermediate patients in chemotherapy-treated patients, while the B and E axes are associated with prognosis among D-low/intermediate patients not treated with chemotherapy. Thus, the D axis confounds the prognostic ability of the other stromal axes irrespective of adjuvant chemotherapy status.

The D axis score and the inherent prognostic difficulty of some patients

The vast majority of previously reported prognostic gene signatures for breast cancer are derived from bulk expression data from the tumor proper. These signatures measure a broad range of tumoral hallmarks and cancer-related processes (e.g., proliferation, genomic instability, immune response). In a previous effort (5), we identified a subset of patients whose observed outcome was consistently mispredicted by almost all reported prognostic gene signatures. As our stromal D axis score appears to be a master controller of the prognostic capacity of the T, B, and E axes, we hypothesized that patient samples estimated high for the D axis might have higher inherent prognostic difficulty, that is, gene signatures that predict prognosis will almost always incorrectly predict D-high patients. If this is true, then knowledge of the state of the D axis might provide a significant breakthrough that would increase the accuracy of prognostic classifiers.

Using the inherent difficulty score from Tofigh and colleagues (5), we observed a significant difference in inherent difficulty of observed poor-outcome patients contingent on D-high versus D-low or D-intermediate status (two sample, two-sided t test; P < 0.05; Fig. 4C, DMFS at 5 years with blue and red representing good and poor outcome, respectively). Thus, those poor-outcome patients that are high for the D axis are systematically mispredicted as good-outcome.

The TNBC stromal axes are generalizable to other patient cohorts

Although the stromal axes were identified within a TNBC cohort, a natural question is to ask whether they have prognostic capacity in other breast cancer subtypes. It is well-established that signatures related to cell cycle and tumor proliferation provide prognostic information especially within ER-positive related cohorts (5, 10, 44), suggesting potential efficacy for the D stromal axis. The prognostic capacity of immune-related signatures are also well established, particularly within ER-negative cohorts, suggesting efficacy for the T and B axes (5, 10, 45, 46). We asked whether there is evidence that the microenvironmental states captured by our TNBC stromal axes can be used universally across the disease to predict disease progression.

Following the methodology of Tofigh and colleagues (5), we trained prognostic predictors for each of the stromal axes and compared these predictors with previously reported predictors (n ∼ 120), trained in the same manner for each breast cancer subtype, across a large collection (n ∼ 5,000) of invasive breast cancer profiles (Fig. 4D; Supplementary Fig. S6). More specifically, for the gene set of each stromal axis and each prognostic signature, a Naive Bayes' Classifier was generated under statistical cross-validation within datasets while reserving several complete datasets for additional independent validation (Supplementary Fig. S7). The use of such classifiers allows, for example, the learning procedure to “weight” specific genes within a gene signature as more or less important in predicting patient outcome. DMFS at 5 years was used as the clinical end point, and performance was measured by the product of accuracy, via standard survival analyses (log-rank test) and via a random sampling–based approach.

We observed that the stromal axis scores, despite being discovered among TNBC patients, were associated with patient prognosis in other patient cohorts. The predictors derived from the stromal axes are observed to cluster closely with other predictors polling similar processes (Fig. 4D; Supplementary Fig. S6). Specifically the D predictor clusters with a predictor derived from comparing normal and tumor stroma (47), the T predictor clusters with other signatures linked to the presence of immune cells in the tumor (48), and the B predictor clusters amongst other B-cell–related signatures (49). Interestingly, the E predictor clusters among predictors of metastatic potential [including metastasis to brain (50) and lung (51)] and among signatures related to immune suppression/protumor immune signaling (52, 53). This suggests that the E axis may represent an early signature in the microenvironment of invasion or metastases. Future work could determine whether the molecular events underlying the E axis could be therapeutically targeted to prevent metastases.

As perhaps expected, the T predictor was one of the best predictors in the TNBC cohort (Fig. 4D). However, it was also significantly associated with prognosis in all subtypes except the ER-positive/HER2-positive cohort. While the D predictor was also significant in these cohorts, it was among the top predictors of prognosis in unstratified pan-breast cancer analysis. Although the D predictor is significant in the TNBC subtype (P < 0.001), it was not among the highest-performing classifiers. It has been previously observed that signatures measuring proliferation have their best performance in unstratified analyses, likely because there are large concomitant differences in the proliferative indices and rate of poor outcome between ER+ and ER patients (5). Given that the D axis score correlates strongly with tumor proliferation, it is perhaps unsurprising that the D predictor is one of the best predictors in unstratified analysis. The B predictor is significantly associated with prognosis in the same cohorts as the T predictor, although it never appears among the highest performing signatures. Interestingly the E predictor showed a stronger association with prognosis in unstratified and ER-positive cohorts, than within the TNBC subtype. Together these results suggest that the properties underlying stromal axis scores are nearly universal to breast cancer, and have different prognostic capacities in different subtypes.

The poor-outcome TNBC subtype remains the focus of much research as we still lack effective classification markers, prognostic signatures, and targeted therapies. This study represents the first large-scale effort to investigate the tumor microenvironment in TNBC patients. With respect to the TNBC subtype, previous studies have focused on gene expression profiling (12–14) or DNA sequencing (54) of bulk material enriched for tumor cells. Efforts to study the tumor microenvironment, including our own (18–20, 27), have used LCM to isolate stromal elements in a pan-breast cancer fashion, not restricted to TNBC.

Recently the spatial and temporal heterogeneity of tumor stroma has received much focus. It has previously been reported that the stromal profile of a single section is sufficient to capture the spatial heterogeneity of the tumor stroma (55). In contrast, it has been observed that the response of the tumor stroma to neoadjuvant treatment is varied between patients, and that the stromal response to treatment is predictive of patient outcome (56). However, the heterogeneity of the stroma within a subtype has not been investigated.

Here, we identify four stromal axes in TNBC patients. A key distinction between this study and previous work is that each patient is scored along all axes, rather than being assigned to individual distinct subtypes. We observe that this method better captures the heterogeneity of the TNBC stroma. Despite being discovered in LCM-derived material, these stromal axes are shown to hold true even when applied to matching bulk expression sample profiles (Fig. 1D; Supplementary Table S5). This allows us to assign levels of the four axes to a large bulk-derived compendium of TNBC patients, revealing that the activation state of the B, T, and E axes were associated with patient outcome, and that the D axis score is associated with tumor proliferation.

Using the infrastructure and concepts from Tofigh and colleagues (5), we show evidence that the vast majority of existing gene signatures for predicting patient prognosis fail for many D-high patient samples. In particular, these signatures often incorrectly predict D-high poor-outcome patients to have good outcome. D-high patients are the least proliferative among the TNBC, although all of these tumors are highly proliferative in comparison with non-TNBC tumors. This suggests that the approximately 40% of TNBC tumors estimated to be D-high and the least proliferative have also been problematic for existing prognostic signatures. Therefore, when our novel D signature assigns a high score to a sample, current prognostic predictors should not be utilized as they will likely fail. Conversely, it is only when a patient sample is deemed low (or intermediate) for D that the T, B, and E axis scores provide additional prognostic information. These observations allow us to generate a subtyping schema that ablates the complexity of having many (81) potential subtypes (Fig. 4E). The variable survival rates observed among D-low patients also explains why the D score did not have prognostic capacity when studied in isolation (Fig. 1E).

When the D score is high, the immune response within the microenvironment (estimated via the T and B scores), and the E score are insufficient to predict outcome. The desmoplastic stroma as encompassed by the D signature may suppress the tumor-antagonistic effects of a stimulated immune response. This is consistent with the observation that the cohort of samples with high D scores have moderate to poor prognosis when compared to the entire TNBC cohort. Hence, the prior or concurrent targeting of desmoplastic stroma may enhance the therapeutic benefit of immunomodulatory therapy in this patient cohort.

Using bulk expression profiles, the TNBCType scheme estimates that six subtypes capture the heterogeneity of TNBC patients. By applying our methodology to the genes that define their subtypes, we present strong evidence that essentially every patient sample belongs to multiple Lehmann and colleagues subtypes. However, as their methodology assigns each patient to exactly one subtype (e.g., MSL and not LAR), these patients would be treated according to the standard of care for the MSL subtype, and potential antiandrogen therapies suitable for the LAR subtype would be ignored despite the evidence that the patient is also LAR positive. Therefore application of our approach may improve the already positive findings from Lehmann and colleagues (13, 15).

Ng and colleagues (10) suggest that TNBC may be characterized by three or possibly four signatures related to androgen signaling, immune infiltration, and a desmoplastic stroma, all of which are captured by our four stromal axes presented here. This effort is the first to offer a quantitative estimate of the number of distinct, populated TNBC subtypes. If there are four core axes of TNBC and each axis is scored as high, intermediate, or low, then there are a total of 34(=81) possible subtypes. We found that 15 of these 81 subtypes had more patients than expected by chance, while 17 subtypes that had significantly fewer patients than expected by chance. The remaining 49 subtypes had populations within our datasets that are not significantly different than what we would expect by chance, and therefore larger TNBC cohorts are necessary to investigate their prognostic and predictive values. As the axes are associated with genomic subtypes, we note that development of an approach similar to AIMS (57) may be required in future.

Our approach has therefore identified four interacting stromal axes that can be combined to subtype TNBC patients. These axis scores have also displayed broader applicability among other cohorts of breast cancer patients, indicating that they play a key role in determining breast cancer outcome. The finding that immune-related (T and B) scores lose prognostic value in D-high cases is of special relevance in triple-negative breast cancer, where various immunotherapy approaches are entering clinical practice. The presence of a cohort comprising approximately 40% of TNBC cases, for which immune status is unrelated to prognosis, may underlie the lack of response to immune-directed agents observed in a majority of patients.

No potential conflicts of interest were disclosed.

Conception and design: A. Omeroglu, M.T. Hallett, M. Park

Development of methodology: S.M.I. Saleh, H. Zhao, A. Omeroglu, M.T. Hallett, M. Park

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): N.Bertos, M. Souleimanova, H. Zhao, A. Omeroglu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.M.I. Saleh, T. Gruosso, M. Gigoux, M.T. Hallett

Writing, review, and/or revision of the manuscript: S.M.I. Saleh, N. Bertos, T. Gruosso, M. Gigoux, A. Omeroglu, M.T. Hallett, M. Park

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N. Bertos, M. Souleimanova, H. Zhao, A. Omeroglu, M.T. Hallett

Study supervision: M.T. Hallett, M. Park

Other (generated figures): S.M.I. Saleh

We thank Vanessa Dumeaux and Daniel Del Balso for their advice during the preparation of this manuscript and for the careful reading of this manuscript. We thank Dr. Marie-Christine Guiot for her assistance with anti-Ki67 immunohistochemistry.

This study was supported by funding from CQDM (Consortium québécois sur la découverte du medicament/Quebec Consortium for Drug Discovery), and the NIH to M. Park and M. Hallett. S.M.I. Saleh is supported by the CIHR Systems Biology Training program and a McGill Faculty of Medicine Internal Studentship (awarded as Gershman Memorial and Hugh E. Burke Fellowships). The breast tissue and data bank at McGill University is supported by funding from the Database and Tissue Bank Axis of the Réseau de Recherche en Cancer of the Fonds de Recherche du Québec-Santé and the Quebec Breast Cancer Foundation (M. Park). The bulk expression profile dataset was generated thanks to financial support provided by Genome Quebec (M. Hallett and M. Park).

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.

1.
Esposito
A
,
Criscitiello
C
,
Curigliano
G
. 
Highlights from the 14(th) St Gallen International Breast Cancer Conference 2015 in Vienna: Dealing with classification, prognostication, and prediction refinement to personalize the treatment of patients with early breast cancer
.
Ecancermedicalscience
2015
;
9
:
518
.
2.
Boyle
P
. 
Triple-negative breast cancer: epidemiological considerations and recommendations
.
Ann Oncol
2012
;
23
:
vi7
vi12
.
3.
Palma
G
,
Frasci
G
,
Chirico
A
,
Esposito
E
,
Siani
C
,
Saturnino
C
, et al
Triple negative breast cancer: looking for the missing link between biology and treatments
.
Oncotarget
2015
;
6
:
26560
74
.
4.
Oakman
C
,
Viale
G
,
Di Leo
A
. 
Management of triple negative breast cancer
.
The Breast
2010
;
19
:
312
21
.
5.
Tofigh
A
,
Suderman
M
,
Paquet
ER
,
Livingstone
J
,
Bertos
N
,
Saleh
SM
, et al
The prognostic ease and difficulty of invasive breast carcinoma
.
Cell Rep [Internet]
; 
2014
.
Available from
: http://www.cell.com/cell-reports/abstract/S2211-1247(14)00765-7.
6.
The Cancer Genome Atlas Network
. 
Comprehensive molecular portraits of human breast tumours
.
Nature
2012
;
490
:
61
70
.
7.
Bergamaschi
A
,
Kim
YH
,
Wang
P
,
Sørlie
T
,
Hernandez-Boussard
T
,
Lonning
PE
, et al
Distinct patterns of DNA copy number alteration are associated with different clinicopathological features and gene-expression subtypes of breast cancer
.
Genes Chromosomes Cancer
2006
;
45
:
1033
40
.
8.
Andreopoulou
E
,
Schweber
SJ
,
Sparano
JA
,
McDaid
HM
. 
Therapies for triple negative breast cancer
.
Expert Opin Pharmacother
2015
;
16
:
983
98
.
9.
Lehmann
BD
,
Pietenpol
JA
. 
Clinical implications of molecular heterogeneity in triple negative breast cancer
.
Breast Edinb Scotl
2015
;
24
:
S36
40
.
10.
Ng
CKY
,
Schultheis
AM
,
Bidard
F-C
,
Weigelt
B
,
Reis-Filho
JS
. 
Breast cancer genomics from microarrays to massively parallel sequencing: paradigms and new insights
.
J Natl Cancer Inst
2015
;
107
:
djv015
.
11.
Jiang
T
,
Shi
W
,
Natowicz
R
,
Ononye
SN
,
Wali
VB
,
Kluger
Y
, et al
Statistical measures of transcriptional diversity capture genomic heterogeneity of cancer
.
BMC Genomics
2014
;
15
:
876
.
12.
Rody
A
,
Karn
T
,
Liedtke
C
,
Pusztai
L
,
Ruckhaeberle
E
,
Hanker
L
, et al
A clinically relevant gene signature in triple negative and basal-like breast cancer
.
Breast Cancer Res
2011
;
13
:
R97
.
13.
Lehmann
BDB
,
Bauer
JJ
,
Chen
X
,
Sanders
ME
,
Chakravarthy
AB
,
Shyr
Y
, et al
Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies
.
J Clin Invest
2011
;
121
:
2750
67
.
14.
Burstein
MD
,
Tsimelzon
A
,
Poage
GM
,
Covington
KR
,
Contreras
A
,
Fuqua
S
, et al
Comprehensive genomic analysis identifies novel subtypes and targets of triple-negative breast cancer
.
Clin Cancer Res
2014
;
21
:
1688
99
.
15.
Masuda
H
,
Baggerly
KA
,
Wang
Y
,
Zhang
Y
,
Gonzalez-Angulo
AM
,
Meric-Bernstam
F
, et al
Differential response to neoadjuvant chemotherapy among 7 triple-negative breast cancer molecular subtypes
.
Clin Cancer Res
2013
;
19
:
5533
40
.
16.
Ring
BZ
,
Hout
DR
,
Morris
SW
,
Lawrence
K
,
Schweitzer
BL
,
Bailey
DB
, et al
Generation of an algorithm based on minimal gene sets to clinically subtype triple negative breast cancer patients
.
BMC Cancer
2016
;
16
:
143
.
17.
Hanahan
D
,
Coussens
LM
. 
Accessories to the crime: functions of cells recruited to the tumor microenvironment
.
Cancer Cell
2012
;
21
:
309
22
.
18.
Finak
G
,
Sadekova
S
,
Pepin
F
,
Hallett
M
,
Meterissian
S
,
Halwani
F
, et al
Gene expression signatures of morphologically normal breast tissue identify basal-like tumors
.
Breast Cancer Res
2006
;
8
:
R58
.
19.
Finak
G
,
Bertos
N
,
Pepin
F
,
Sadekova
S
,
Souleimanova
M
,
Zhao
H
, et al
Stromal gene expression predicts clinical outcome in breast cancer
.
Nat Med
2008
;
14
:
518
27
.
20.
Boersma
BJ
,
Reimers
M
,
Yi
M
,
Ludwig
JA
,
Luke
BT
,
Stephens
RM
, et al
A stromal gene signature associated with inflammatory breast cancer
.
Int J Cancer
2008
;
122
:
1324
32
.
21.
R Core Team
.
R: A Language and Environment for Statistical Computing [Internet]
.
Vienna, Austria
:
R Foundation for Statistical Computing
; 
2015
.
Available from
: http://www.R-project.org/.
22.
Ritchie
ME
,
Phipson
B
,
Wu
D
,
Hu
Y
,
Law
CW
,
Shi
W
, et al
limma powers differential expression analyses for RNA-sequencing and microarray studies
.
Nucleic Acids Res
2015
;
43
:
e47
.
23.
Suzuki
R
,
Shimodaira
H
.
pvclust: Hierarchical Clustering with P-Values via multiscale bootstrap resampling [Internet]
; 
2014
.
Available from
: http://CRAN.R-project.org/package=pvclust.
24.
Paquet
ER
,
Lesurf
R
,
Tofigh
A
,
Dumeaux
V
,
Hallett
MT
. 
Detecting gene signature activation in breast cancer in an absolute, single-patient manner
.
Breast Cancer Res
. 
2017
;
19
:
48
.
25.
Parker
JS
,
Mullins
M
,
Cheang
MCU
,
Leung
S
,
Voduc
D
,
Vickery
T
, et al
Supervised risk predictor of breast cancer based on intrinsic subtypes
.
J Clin Oncol
2009
;
27
:
1160
7
.
26.
Chen
X
,
Li
J
,
Gray
WH
,
Lehmann
BD
,
Bauer
JA
,
Shyr
Y
, et al
TNBCtype: a subtyping tool for triple-negative breast cancer
.
Cancer Inform
2012
;
11
:
147
56
.
27.
Pepin
F
,
Bertos
N
,
Laferrière
J
,
Sadekova
S
,
Souleimanova
M
,
Zhao
H
, et al
Gene-expression profiling of microdissected breast cancer microvasculature identifies distinct tumor vascular subtypes
.
Breast Cancer Res
2012
;
14
:
R120
.
28.
Jabri
B
,
Abadie
V
. 
IL-15 functions as a danger signal to regulate tissue-resident T cells and tissue destruction
.
Nat Rev Immunol
2015
;
15
:
771
83
.
29.
Cullen
SP
,
Brunet
M
,
Martin
SJ
. 
Granzymes in cancer and immunity
.
Cell Death Differ
2010
;
17
:
616
23
.
30.
Zitvogel
L
,
Galluzzi
L
,
Kepp
O
,
Smyth
MJ
,
Kroemer
G
. 
Type I interferons in anticancer immunity
.
Nat Rev Immunol
2015
;
15
:
405
14
.
31.
Takeda
K
,
Akira
S
. 
STAT family of transcription factors in cytokine-mediated biological responses
.
Cytokine Growth Factor Rev
2000
;
11
:
199
207
.
32.
Chang
H-C
,
Han
L
,
Goswami
R
,
Nguyen
ET
,
Pelloso
D
,
Robertson
MJ
, et al
Impaired development of human Th1 cells in patients with deficient expression of STAT4
.
Blood
2009
;
113
:
5887
90
.
33.
Kumar
D
,
Whiteside
TL
,
Kasid
U
. 
Identification of a novel tumor necrosis factor-α-inducible Gene, SCC-S2, containing the consensus sequence of a death effector domain of fas-associated death domain-like interleukin- 1β-converting Enzyme-inhibitory Protein
.
J Biol Chem
2000
;
275
:
2973
8
.
34.
Sarma
V
,
Wolf
FW
,
Marks
RM
,
Shows
TB
,
Dixit
VM
. 
Cloning of a novel tumor necrosis factor-alpha-inducible primary response gene that is differentially expressed in development and capillary tube-like formation in vitro
.
J Immunol
1992
;
148
:
3302
12
.
35.
Tai
S-K
,
Tan
OJ-K
,
Chow
VT-K
,
Jin
R
,
Jones
JL
,
Tan
P-H
, et al
Differential expression of metallothionein 1 and 2 Isoforms in breast cancer lines with different invasive potential: identification of a novel nonsilent metallothionein-1H mutant variant
.
Am J Pathol
2003
;
163
:
2009
19
.
36.
Cheung
KJ
,
Ewald
AJ
. 
A collective route to metastasis: Seeding by tumor cell clusters
.
Science
2016
;
352
:
167
9
.
37.
Gandellini
P
,
Andriani
F
,
Merlino
G
,
D'Aiuto
F
,
Roz
L
,
Callari
M
. 
Complexity in the tumour microenvironment: cancer associated fibroblast gene expression patterns identify both common and unique features of tumour-stroma crosstalk across cancer types
.
Semin Cancer Biol
2015
;
35
:
96
106
.
38.
Gilkes
DM
,
Semenza
GL
,
Wirtz
D
. 
Hypoxia and the extracellular matrix: drivers of tumour metastasis
.
Nat Rev Cancer
2014
;
14
:
430
9
.
39.
Beck
AH
,
Espinosa
I
,
Gilks
CB
,
van de Rijn
M
,
West
RB
. 
The fibromatosis signature defines a robust stromal response in breast carcinoma
.
Lab Invest
2008
;
88
:
591
601
.
40.
Lehmann
BD
,
Bojana
Jovanović
,
Chen
X
,
Estrada
MV
,
Johnson
KN
,
Shyr
Y
, et al
Refinement of triple-negative breast cancer molecular subtypes: implications for neoadjuvant chemotherapy selection
.
PLOS ONE
2016
;
11
:
e0157368
.
41.
Nielsen
TO
,
Parker
JS
,
Leung
S
,
Voduc
D
,
Ebbert
M
,
Vickery
T
, et al
A Comparison of PAM50 intrinsic subtyping with immunohistochemistry and clinical prognostic factors in tamoxifen-treated estrogen receptor–positive breast cancer
.
Clin Cancer Res
2010
;
16
:
5222
32
.
42.
Barnard
NJ
,
Hall
PA
,
Lemoine
NR
,
Kadar
N
. 
Proliferative index in breast carcinoma determined in situ by Ki67 immunostaining and its relationship to clinical and pathological variables
.
J Pathol
1987
;
152
:
287
95
.
43.
Curtis
C
,
Shah
SP
,
Chin
S-F
,
Turashvili
G
,
Rueda
OM
,
Dunning
MJ
, et al
The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups
.
Nature
2012
;
486
:
346
52
.
44.
Dai
H
,
Veer
L van't
,
Lamb
J
,
He
YD
,
Mao
M
,
Fine
BM
, et al
A cell proliferation signature is a marker of extremely poor outcome in a subpopulation of breast cancer patients
.
Cancer Res
2005
;
65
:
4059
66
.
45.
Teschendorff
AE
,
Miremadi
A
,
Pinder
SE
,
Ellis
IO
,
Caldas
C
. 
An immune response gene expression module identifies a good prognosis subtype in estrogen receptor negative breast cancer
.
Genome Biol
2007
;
8
:
R157
.
46.
West
NR
,
Milne
K
,
Truong
PT
,
Macpherson
N
,
Nelson
BH
,
Watson
PH
. 
Tumor-infiltrating lymphocytes predict response to anthracycline-based chemotherapy in estrogen receptor-negative breast cancer
.
Breast Cancer Res
2011
;
13
:
R126
.
47.
Planche
A
,
Bacac
M
,
Provero
P
,
Fusco
C
,
Delorenzi
M
,
Stehle
J-C
, et al
Identification of prognostic molecular features in the reactive stroma of human breast and prostate cancer
.
PLOS ONE
2011
;
6
:
e18640
.
48.
Alexe
G
,
Dalgin
GS
,
Scanfeld
D
,
Tamayo
P
,
Mesirov
JP
,
DeLisi
C
, et al
High expression of lymphocyte-associated genes in node-negative HER2+ breast cancers correlates with lower recurrence rates
.
Cancer Res
2007
;
67
:
10669
76
.
49.
Watkins
NA
,
Gusnanto
A
,
Bono
B de
,
De
S
,
Miranda-Saavedra
D
,
Hardie
DL
, et al
A HaemAtlas: characterizing gene expression in differentiated human blood cells
.
Blood
2009
;
113
:
e1
9
.
50.
Bos
PD
,
Zhang
XH-F
,
Nadal
C
,
Shu
W
,
Gomis
RR
,
Nguyen
DX
, et al
Genes that mediate breast cancer metastasis to the brain
.
Nature
2009
;
459
:
1005
9
.
51.
Landemaine
T
,
Jackson
A
,
Bellahcène
A
,
Rucci
N
,
Sin
S
,
Abad
BM
, et al
A six-gene signature predicting breast cancer lung metastasis
.
Cancer Res
2008
;
68
:
6092
9
.
52.
Galon
J
,
Costes
A
,
Sanchez-Cabo
F
,
Kirilovsky
A
,
Mlecnik
B
,
Lagorce-Pagès
C
, et al
Type, density, and location of immune cells within human colorectal tumors predict clinical outcome
.
Science
2006
;
313
:
1960
4
.
53.
Tosolini
M
,
Kirilovsky
A
,
Mlecnik
B
,
Fredriksen
T
,
Mauger
S
,
Bindea
G
, et al
Clinical impact of different classes of infiltrating T cytotoxic and helper cells (Th1, Th2, Treg, Th17) in patients with colorectal cancer
.
Cancer Res
2011
;
71
:
1263
71
.
54.
Shah
SP
,
Roth
A
,
Goya
R
,
Oloumi
A
,
Ha
G
,
Zhao
Y
, et al
The clonal and mutational evolution spectrum of primary triple-negative breast cancers
.
Nature
2012
;
486
:
395
9
.
55.
Mani
NL
,
Schalper
KA
,
Hatzis
C
,
Saglam
O
,
Tavassoli
F
,
Butler
M
, et al
Quantitative assessment of the spatial heterogeneity of tumor-infiltrating lymphocytes in breast cancer
.
Breast Cancer Res
2016
;
18
:
78
.
56.
García-Martínez
E
,
Gil
GL
,
Benito
AC
,
González-Billalabeitia
E
,
Conesa
MAV
,
García
TG
, et al
Tumor-infiltrating immune cell profiles and their change after neoadjuvant chemotherapy predict response and prognosis of breast cancer
.
Breast Cancer Res
2014
;
16
:
488
.
57.
Paquet
ER
,
Hallett
MT
. 
Absolute assignment of breast cancer intrinsic molecular subtype
.
J Natl Cancer Inst
2015
;
107
:
357
.

Supplementary data