Interactions between immunity, proliferation and molecular subtype in breast cancer prognosis

1Expression profiling studies in human tumors have enabled new insights into the genes and pathways that contribute to tumorigenesis and spurred the development of gene expression signatures prognostic of patient outcomes. Genes comprising prognostic signatures often provide clues to the pathobiological mechanisms that drive cancer progression. With the aim of discovering genes with statistical associations with breast cancer recurrence, we and others have identified a number of genes with roles in cellular proliferation, including multi-gene proliferation signatures that directly reflect tumor proliferative capacity. These signatures are highly significantly associated with poor patient outcomes, consistent with the view that uncontrolled cell proliferation is a central feature of neoplastic disease and, ultimately, a contributing factor in metastatic progression. Indeed, proliferation-associated genes are common components of many previously reported prognostic gene signatures, including Genomic Health's 21-gene Oncotype Dx test (Genomic Health, Inc., Redwood City, CA, USA), and frequently account for the majority of the prognostic power driving the performance of these signatures. Thus, a clear biological understanding of how prognostic genes relate to different aspects of tumor pathobiology is imperative to both the optimal construction of prognostic models and the elucidation of key regulators of cancer behavior.

2In recent years, we and others have observed that elevated expression levels of many genes involved in immune response pathways are associated with reduced risk of breast cancer recurrence. These observations support the view that cancer-leukocyte interactions in the microenvironment of established tumors may function to limit the growth and metastatic progression of breast cancer. However, the extent to which these genes reflect different effector cell populations, or contribute to patient prognosis in the presence of other predictive biomarkers such as proliferation, remains unclear.

3In this report, we investigate the biological origins of coordinately expressed genes in breast cancer that exhibit statistical associations with patient distant metastasis-free survival (DMFS). We identify gene clusters indicative of tumor-immune cell interactions that organize into three distinct immunity-related gene signatures, or metagenes, and shed light on their prognostic implications for tumors of differing proliferative capacity with an emphasis on highly proliferative breast cancers and the most aggressive intrinsic molecular subtypes in particular.

4To characterize prognostic gene modules, we created a multi-study microarray database of 2,116 breast tumor expression profiles of which 1,954 were annotated with corresponding clinicopathological data including DMFS (See Additional file 1 for clinical details). To facilitate gene discovery, we randomized the dataset across study groups and clinical features into two equivalent patient subpopulations, termed patient groups 977A and 977B (Table 1). In each patient group, Cox proportional hazards regression was conducted to identify genes with statistically significant associations with DMFS while controlling for false discoveries (q < 0.1). The analysis identified 3,094 significant gene probe sets in 977A and 3,304 in 977B (gene details provided in Additional file 2). In parallel, the DMFS-associated genes identified in each patient group were hierarchically clustered to enable analysis of gene correlation structure (Figure 1 and Additional file 3). As anticipated, a proliferation gene cluster was readily identifiable in both patient groups. This cluster of genes has been previously described in multiple studies as being significantly associated with patient survival, and consists of the highly correlated group of cell cycle genes associated with markers of tumor cell proliferation. In a subset analysis, we examined the correlation between this proliferation gene cluster and clinical markers of proliferation. As expected, we observed a strong positive correlation between the average expression of the genes comprising this cluster and Ki67 staining (by MIB1 antibody) and mitotic index (Additional file 4), consistent with the notion that these genes quantify tumor proliferative capacity.

5Further inspection of the cluster architecture revealed a large reproducible cluster of genes associated with immune cell functions that exhibited negligible correlation with the proliferation cluster (Figure 1 and Additional file 3). Gene ontology (GO) enrichment analysis of the genes within this large cluster showed highly significant enrichment for numerous immune cell processes including lymphocyte activation, antigen processing and presentation, positive regulation of immune system process, and other annotations specific for different immune cell lineages (P < 0.0001, false discovery rate (FDR)-adjusted; Additional file 5). Closer inspection of the nested correlation structure revealed distinct gene 'subclusters' that were highly reproducible between patient groups (Figure 1 and Additional file 3). While one predominant proliferation cluster emerged, three distinct immune gene subclusters (termed immune subclusters #1, #2 and #3) could be discerned in both 977A and 977B. To investigate the underlying biology associated with these subclusters, genes comprising each subcluster were selected from the dendrogram branches using a threshold of average Pearson correlation of 0.6 (see Methods). The number of gene probe sets per subcluster ranged from 20 to 59, and details regarding their subcluster membership are shown in Additional file 6. Although the genes comprising the subclusters were independently selected from 977A and 977B (based on correlation structure alone), we observed a high degree of probe overlap when comparing subclusters across the two groups (Additional file 7A). The majority of probes identified within a cluster of one patient group were also found within the cognate cluster of the other patient group, though three genes were observed to exhibit cluster inconsistency (associated with immune subcluster #2 in one patient group and immune subcluster #3 in the other). For a more decisive comparison of the expression patterns of the cognate clusters, we examined the correlation between cognate clusters of 977A and 977B. We observed near perfect correlations between cognate clusters with Pearson correlation coefficients (r) ranging from 0.97 to 0.99 (Additional file 7B). For the immune subclusters, this indicated that the hierarchical organization of the genes into three discernible expression vectors was a reproducible event.

6We hypothesized that the immune gene subclusters likely reflect the relative abundance of tumor-infiltrating immune cells. We employed several strategies to investigate this hypothesis. First, we investigated the expression patterns of the immune cluster genes in the microarray dataset of Abbas and colleagues, comprising a comprehensive collection of human leukocyte gene expression profiles. Strikingly, the nested correlation structure and gene composition of the immune gene subclusters, as observed in the breast tumors, remained largely unaltered in this pan-leukocyte dataset after hierarchical clustering (Figure 2; also see Additional file 8 for greater detail). Consistent with our hypothesis, we found that immune cluster #1 consisted of genes (mostly immunoglobulin-encoding genes) highly and exclusively expressed in B cell/plasma cell populations (hence termed the B/P Cluster). By contrast, expression of genes in immune cluster #2 (such as components of the T cell receptor-CD3 complex and granzymes) were found to be mostly restricted to T cells and natural killer cells (hence termed the T/NK Cluster), whereas the genes of immune cluster #3 (including major histocompatibility complex (MHC) class II (human leukocyte antigen; HLA) and myeloid-specific markers (for example, colony stimulating factor 1 receptor)) were most consistently expressed at highest levels in monocytes and dendritic cells (hence termed the M/D Cluster).

7Next, we examined the immune gene subclusters for gene-level enrichment of GO terms. Numerous highly significant biological annotations emerged that were consistent with our observations in the Abbas leukocyte dataset. Representative GO terms selected from among the top 10 most significant terms for each subcluster are shown in Table 2. The B/P cluster was highly enriched for variable region immunoglobulin genes involved in antigen binding - consistent with B cell/plasma cell biology. The T/NK cluster was enriched for terms consistent with the positive regulation of lymphocyte activation and differentiation, T cell signaling and natural killer cell functions. The M/D cluster was enriched for significant terms associated with MHC class II-mediated antigen processing and presentation - characteristic of macrophages and dendritic cells.

8We then tested for direct associations between the magnitude of expression of the immune gene subclusters and the relative abundance of tumor-infiltrating leukocytes. To reduce dimensionality of the gene expression data, we averaged the gene signal intensities within each gene subcluster according to the method of Dave and colleagues to generate a 'metagene' expression value for each breast cancer case. The immune metagene values were then compared to measurements of immune cell infiltrate assessed in tumor sections (n = 35) at another institution. Significant positive trends between metagene values and immune cell abundance were observed for each metagene (B/P, P = 0.08; T/NK, P = 0.02; M/D, P = 0.009; Additional file 9A). Additionally, we extrapolated the immune metagene concept to a more quantitative RNA analysis platform to investigate how the concept might be generalized to a diagnostic setting. Genes representative of the B/P and T/NK metagenes were profiled prospectively in a panel of estrogen receptor-positive (ER+), formalin-fixed paraffin-embedded (FFPE) breast tumor sections using the Panomics QuantiGene Plex 2.0 RNA assay system (Affymetrix, Santa Clara, CA, USA). Expression levels of the selected genes were found to be positively and significantly correlated with total leukocyte counts (B/P, P = 0.005; T/NK, P = 0.02; Additional file 9B-D). Taken together, these findings support the view that the immune gene subclusters reflect the relative abundance of infiltrating immune cell populations.

9We examined the prognostic relationships between the proliferation and immune metagenes. First, the metagene expression values were used to divide breast cancer cases into population tertiles. This procedure is illustrated in Figure 3 where patient group 977A is shown divided into proliferation metagene tertiles then further stratified into low (PL), intermediate (PI), and high (PH) expression tertiles by the B/P metagene. Kaplan-Meier plots of the DMFS of patients classified by the B/P tertiles are shown. Strikingly, we found that the prognostic power of the B/P metagene, while distinct from that of proliferation, is dependent on the proliferative status of the tumor. Specifically, we observed that its prognostic power resides exclusively in the highly proliferative tumors, as defined by the upper proliferation metagene tertile (Figure 3E). To investigate the robustness of this phenomenon and to reduce the potential for data overfitting, we used each patient group (977A and 977B) in both training and testing scenarios. For example, using group 977A as a training set, the gene content of the proliferation and the immune metagenes were defined and their corresponding expression tertile cut-points were determined. These metagenes and tertile cut-points were used to group 977B cases into low, intermediate and high expression tertiles for survival analysis. Shown in Figure 4 are the cross-group test results for each of the immune metagenes. Consistently, we observed that all three immune metagenes displayed a highly significant positive association with DMFS that is conditionally prognostic - dependent on the high proliferation phenotype defined by the upper tertile of the proliferation metagene (the PH tertile).

10Although the immune metagenes form distinct gene subclusters, their expression profiles are intrinsically correlated (that is, they are subcomponents of a larger immune gene cluster), suggesting the possibility that the metagenes could exhibit prognostic redundancy, as previously hypothesized. To address this question, we compared the prognostic significance of the immune metagenes to one another via multivariable analysis. We constructed Cox regression models inclusive of pair-wise combinations of the metagenes or all three metagenes, simultaneously (Table 3). In all pair-wise comparisons (models 1 to 3), the metagenes contributed significant independent prognostic information reflective of their non-redundant contributions to prognosis. In a fully combined model (model 4), the B/P and M/D metagenes exhibited the greatest non-redundant prognostic power. Additionally, we assessed the prognostic contributions of the immune metagenes in the presence of conventional variables including nodal status, T stage (tumor size), histologic grade, age, ER status and treatment status (Table 4). While the majority of variables showed moderately to highly significant associations with DMFS by univariable analysis, in the combined model, only the B/P metagene, nodal status, tumor size and treatment remained significant, with greatest significance observed for the B/P metagene (P = 0.0001). Together, these findings demonstrate that the immune gene signatures capture distinct aspects of patient prognosis with the B/P signature, in particular, imparting the most significant and additive prognostic power in highly proliferative breast cancer compared to the other immune metagenes and conventional prognostic markers.

11Next we investigated the impact of the immune metagenes on conventional clinical breast cancer subtypes (ER+ or ER-) and the Sorlie-Perou intrinsic molecular subtypes. First, we examined the distribution of the molecular subtypes as a function of the proliferation metagene (Additional file 10). As expected, the least aggressive subtype, luminal A (LumA,) was found predominately in the low and intermediate proliferation tertiles, whereas the more aggressive luminal B (LumB), Basal-like, and human epidermal growth factor receptor 2-enriched (HER2-E) subtypes were most abundant in the high proliferation tertile. When analyzed for associations with DMFS, all three immune metagenes retained significant prognostic power in the Basal-like, LumB and HER2-E subtypes (PH tertile). This is illustrated by the B/P metagene in Additional file 10, and described further in Table 5.

12In light of recent work illuminating prognostic roles for immune-related genes in specific pathological contexts such as ER- or HER2+ breast cancer, we asked whether the prognostic performance of the immune metagenes was exclusive to the high proliferation tertile in specific tumor subtypes (Table 5). In clinical ER- tumors and Basal-like breast tumors alike, all three immune metagenes were positively associated with DMFS in the high proliferation tertile (PH) and also the intermediate proliferation tertile (PI) - the latter observation indicating that tumor subtype modifies the proliferation dependency of the immune metagenes' prognostic impact. In LumB tumors, the T/NK and M/D metagenes (but not B/P) also trended towards or reached significance, respectively, in the PI tertile, whereas no metagene achieved significance in the PI tertile of the ER+, LumA or HER2-E tumor subtypes. The Claudin-Low (CL) subtype is a rare subset of Basal-like breast tumors with distinguishing features such as high immune cell infiltrate, stem cell-like features and properties characteristic of epithelial-to-mesenchymal transition. We identified 92 CL tumors in our dataset. Unlike other Basal-like tumors, which tend to be highly proliferative, we observed a fairly uniform distribution of CL tumors across the proliferation tertiles, and as expected, the CL tumors showed a bias towards belonging to the upper tertiles of the immune metagenes (data not shown). However, the immune metagenes were not found to be prognostic in CL tumors as a whole, nor in the intermediate or high proliferation tertiles (comprising only 22 cases and 39 cases, respectively). In the low proliferation tertile (PL), we observed an unexpected inverse survival association for some immune metagenes, such as the T/NK metagene, which achieved statistical significance in the low-proliferating ER-, LumB and CL subtypes. Together, these data suggest that the prognostic impact of the immune metagenes in breast cancer are both proliferation- and subtype-dependent, and may signal either good or poor outcome depending on the tumor's proliferative configuration and subtype context.

13Given the potential of the immune metagenes for additive prognostic effects, we asked if simple tertile-based metrics might shed light on the prognostic interplay between the metagenes and if such interactions could form the basis of an integrated model for patient prognosis in treatment-specific contexts. Focusing on the PH tertile (n = 657), we explored the prognostic attributes of specific combinations of low and high tertiles among the three immune metagenes, without applying mathematical optimization or weighting strategies. As shown in Figure 5A (left panel), we observed that patients having one, two or three low immune tertiles all had relatively poor outcomes, ranging from 41% to 52% DMFS at 8 years. No significant survival differences were observed between patients having one, two or three low immune tertiles. Conversely, having high tertiles for all three metagenes was significantly more favorable than having only two or one high tertile assignments (middle panel). Moreover, having high tertiles for all three metagenes was statistically significantly more favorable than having two high tertiles plus an intermediate tertile (that is, for the remaining metagene) or having two high tertiles plus a low tertile (right panel). These observations suggest that a tumor exhibiting a low tertile for any one of the three immune metagenes portends a poor survival outcome that trumps the benefit of having one, or even two, high immune tertiles among the other two metagenes.

14Next, we examined how this classification model might impact patient prognosis in specific therapeutic populations (Figure 5B). Lymph node-negative (LN-) patients who did not receive adjuvant therapy after surgery (left panel) exhibited a marked reduction in 10-year DMFS if their tumors displayed one or more low immune tertiles (green survival curve). In ER+, LN- patients who received tamoxifen monotherapy, the group with consistently high immune tertiles (red curve, middle panel) exhibited highly favorable outcomes (> 90% 10-year DMFS). A group with similarly favorable prognosis (identified by three high immune tertiles) was also observed in patients with highly proliferative ER- and Basal-like breast cancer who received adjuvant chemotherapy (red curve, right panel). The majority of these cases would be clinically classified as triple negative breast cancer - a particularly aggressive and treatment-limited form of the disease. These data suggest that the classification of patients according to immune metagenes can impact patient prognosis in ways that could influence treatment decisions for certain therapeutic subgroups.

15The immune contexture of human cancer, defined as the abundance, location and functional orientation of tumor-infiltrating immune cells, is gaining recognition as a principal determinant of the biological and clinical behavior of many cancer types. Although it is well-established that the immune contexture may elicit both pro- and anti-tumorigenic responses, a growing body of evidence indicates that the presence of abundant tumor-infiltrating leukocytes, within established tumors, foretells favorable prognosis. This association has been rigorously documented for a number of malignancies, most notably cancers of the skin, ovary, colon and breast, underscoring the broad protective effects of anti-cancer immunosurveillance.

16In this work, we investigated the prognostic relevance of transcriptomic footprints of the immune contexture of breast cancer and identified both immune and biological configurations of breast cancer with distinct prognostic attributes. Historically, immunohistochemical measures of the relative abundance of infiltrating immune cells in breast tumors, viewed as non-specific infiltrate or as specific leukocyte subpopulations (such as CD8+ T cells), have led to some controversy with regard to the role of the immune system in patient prognosis. However, prominent immune cell infiltrate observed within late-stage, high-grade, or lymph node-positive breast cancers has consistently been associated with recurrence-free survival of patients. More recently, we and others have employed bioinformatic strategies to investigate the biological underpinnings of genes associated with breast cancer outcomes. A common finding among these studies was the favorable prognosis associated with high expression of various immune-related gene cassettes representing admixed immune cell populations or B cell-enriched or T cell-enriched cell populations, specifically, among ER- or HER2+ breast cancer patients. In the current work, we demonstrate for the first time that a proliferation metagene reflecting tumor proliferative capacity can sharply demarcate breast cancer cases into proliferative subclasses (low, intermediate and high) where the prognostic attributes of immune gene signatures are differentially manifested.

17We identified three distinct expression vectors, or metagenes, within breast tumors that distinguish different tumor-infiltrating leukocyte populations: the B/P metagene (B cells/plasma cells); the T/NK metagene (T cells/natural killer cells); and the M/D metagene (monocytes/dendritic cells). While analysis at the population level revealed that the prognostic power of each immune metagene was uniformly restricted to tumors comprising the PH tertile, analysis by intrinsic subtypes further defined the prognostic orientation of the immune metagenes. For example, the immune metagenes were not associated with DMFS in CL and LumA subtypes of the PH tertile, although the small number of cases examined (n = 39 and n = 20, respectively) may have been statistically limiting. Conversely, the immune metagenes were significantly prognostic in both the PH and PI tertiles in the ER-, Basal-like and LumB subtypes. An interesting and unexpected finding was the prognostic implications of the immune metagenes in the PL ER-, Basal-like and LumB tumors comprising 13%, 9% and 7% of their respective populations. Not only were the immune metagenes not associated with favorable prognosis in the PL tertile, but statistically significant poor-outcome associations were observed, particularly for the T/NK metagene. That these observations were made in relatively small sample populations (ranging from 26 to 51 cases) necessitates caution when interpreting the results. However, that the poor-outcome association of the T/NK metagene achieved statistical significance in both the ER- and LumB (ER+) subpopulations suggests the possibility of the existence of a low proliferation-associated, ER-independent tumor phenotype where T cell and/or natural killer cell abundance may signify pro-metastatic rather than anti-metastatic behavior. By contrast, in the LumA tumors of the PL tertile (n = 347), all three immune metagenes trended towards associations with favorable DMFS, with the T/NK metagene achieving statistical significance. Together, these observations paint a complex picture of how tumor-immune cell interactions regulate malignant progression and suggest that the pro- or anti-tumorigenic properties of infiltrating immune cells vary not only with the proliferative status of the tumor, but are determined, in part, by factors associated with intrinsic subtype.

18How tumor proliferation rate relates to pro- or anti-cancer immune cell behavior has, to our knowledge, not been studied; however, it is plausible that proliferation status could act as a surrogate for one or more immunomodulatory pathological contexts. For example, it has been widely observed that, in breast cancer, the rates of proliferation and cell death are positively correlated. Apoptosis and necrosis are associated with both enhanced lymphocytic infiltrate in breast cancer and enhanced immunogenic response. Thus, cell death that increases available tumor antigen may attract antigen-presenting cells that in turn recruit and/or activate T and B cells in the tumor. Furthermore, increased angiogenesis, which supports increased proliferation, may allow better tumor access by immune cells. These possibilities may explain, in part, the high rate of recurrence-free survival observed in the high immune metagene tertiles of highly proliferative breast cancer. Furthermore, a reduced proliferative (and apoptotic) capacity could reflect a tumor microenvironment more conducive to immunosuppression, and subsequently, poorer survival outcomes. In such an instance, for example, CD4+/FOXP3+ T regulatory cells may predominate over CD8+ cytotoxic T cells. The abundance and location of tumor-infiltrating T regulatory cells, as well as their ratio with cytotoxic T cells, have previously been shown to associate with poor breast cancer outcomes.

19Most previously described immune gene signatures discovered in microarray analyses (including the metagenes described herein) trace back to a common origin for their discovery: a gene cluster of approximately 600 genes highly expressed by tumor-infiltrating leukocytes and whose expression patterns form a larger, diverse immune gene cluster when analyzed in bulk tumor tissues. However, the different gene selection methodologies used, and the variation in size and composition of patient populations examined, may together explain the diversity in gene make-up across the reported signatures, as well as conflicting observations regarding the prognostic performances of similarly derived immune gene cassettes. For example, we and others have observed that not all genes of this larger immune gene cluster are prognostic of breast cancer survival, with some carrying substantially more prognostic weight than others. Although the unsupervised gene selection methods used in previous studies have demonstrated predictive power of the immune genes, the supervised strategy we employed (that is, selecting genes with significant DMFS associations prior to metagene construction) enabled the parsing out of the immune genes with greatest prognostic strength. These genes may point not only to specific cellular components of the immune contexture, but also to immunological functional orientations required for tumor rejection. A more precise assessment of the cellular origins and functional attributes of these genes is needed.

20Rody and colleagues deconstructed the larger immune gene cluster into seven metagenes that appeared to reflect various components of the immune system. In multivariable analysis, only a T cell-related metagene was found to be significantly correlated with disease-free survival (P = 0.01) when considering a mixed population of 1,263 patients with breast cancer; consequently, this was the metagene carried forward for further analysis in ER and HER2 status-specific populations. By contrast, our multivariable analysis revealed that multiple immune metagenes may contribute additive prognostic information when considered in combination with one another. This is likely due in part to our focus on immune genes with a priori associations with DMFS, as well as the fact that our analysis was confined to the highly proliferative breast tumors as defined by our proliferation metagene, while that of Rody and colleagues was not restricted to the more proliferative cases. Furthermore, the composition of our immune metagenes also varied with those of Rody and colleagues. Although 80% of the probe sets comprising our B/P metagene overlapped with the Rody IgG (B cell) metagene, only 60% of our T/NK probe sets overlapped with the Rody LCK (T cell) metagene, and our M/D metagene comprised of novel probe sets not selected by their methods.

21Our observation that the prognostic attributes of the immune metagenes are largely non-redundant may reflect the importance of cooperative interplay among different immune cell types in metastasis-protective immunity. Indeed, proteins critical for such interactions are evidenced in the composition of our metagenes. For example, CD27, a component of the T/NK metagene, encodes a type I transmembrane protein of the tumor necrosis factor receptor superfamily that plays key roles in the expansion and memory of activated CD8+ killer T cells as well as B-cell activation and immunoglobulin synthesis. In natural killer cells, high expression of CD27 is associated with greater effector function and enhanced interaction with dendritic cells compared to CD27-low natural killer cells that exhibit a higher stimulation threshold and express inhibitory receptors. Moreover, the M/D metagene comprises a number of MHC class II (HLA) alpha and beta chain paralogs expressed by professional antigen-presenting cells. The products of these genes present extracellular antigens to T lymphocytes thereby stimulating expansion of T helper cells and, subsequently, the downstream activation of plasma B cells. If such interactions are essential to the maintenance of DMFS, this could explain our observation that patients having even a single low tertile immune metagene assignment are unlikely to achieve a durable remission (Figure 5A).

22Tertile-based cut-points were used in our analysis to characterize prognostic interactions in broad terms, not to develop optimized prognostic classifiers. Nevertheless, we observed that the tertiles could be used to identify significant therapy-relevant risk groups. ER+, LN- breast cancer is frequently treated with hormonal therapy alone or in combination with chemotherapy. In recent years, the decision to withhold chemotherapy from a fraction of these patients has been justified by the 21-gene Oncotype Dx test (Genomic Health, Inc.) which relies on a gene-based classification algorithm. Because proliferation genes carry the greatest prognostic weight in this algorithm, it is not surprising that virtually all of the highly proliferative cases are assigned to the high and intermediate risk groups where the use of chemotherapy is indicated. Interestingly, we observed that about 23% of the highly proliferative, ER+, LN- cases possessed high tertiles for all three immune metagenes, and subsequently exhibited excellent 10-year DMFS following tamoxifen monotherapy. This high survival rate of this group (> 90% at 10 years) is similar to the disease-free survival rate of the Oncotype Dx low-risk group. This indicates that the immune metagenes may have value in identifying a second low-risk fraction of patients from among the highly proliferative (high Oncotype Dx recurrence score) cases who might also be spared unnecessary chemotherapy. Furthermore, this same high tertile immune metagene profile identified patients with ER-, Basal-like breast cancer (predominantly triple negative breast cancer) that would have excellent 10-year survival following adjuvant chemotherapy. By contrast, cases associated with one or more low immune metagenes exhibited > 50% probability of distant metastasis before 5 years. Thus, the immune metagenes could aid in the selection of patients most in need of, and most suitable for the testing of, new therapeutic agents being evaluated in clinical trials. It should be noted that these prognostic observations may reflect some bias related to historical treatment standards such as the use of adjuvant CMF (cyclophosphamide, methotrexate and 5-fluorouracil), FAC (5-fluorouracil, doxorubicin and cyclophosphamide) and AC (cyclophosphamide and doxorubicin). In recent years, the addition of taxanes to chemotherapeutic regimens has reduced the rate of breast cancer relapseby 10% to 15%. Thus, the extent to which the immune metagenes would retain prognostic value in light of today's taxane-inclusive regimens warrants prospective evaluation. Furthermore, mounting evidence that anthracyclines and taxanes both possess immunomodulatory activity that impacts treatment efficacy provides further rationale for investigating the clinical importance of the immune metagenes in patients treated with chemotherapy. Finally, further work will be necessary to determine the optimal diagnostic platform for immune metagene assessment, as well as the precise metagene thresholds that provide maximal prognostic utility within specific breast tumor subtypes.

23Gene expression profiles that quantify immune cell abundance within breast tumors are prognostic of DMFS. The prognostic value of these signatures does not manifest in all tumor types, but rather can be stratified by tumor proliferative capacity in a manner that further depends on molecular subtype. Our findings suggest that multiple immune metagenes measured in combination may provide actionable prognostic information for the most aggressive breast cancer phenotypes for which prognostic assays remain lacking. This work sheds new light on the roles of tumor-infiltrating immune cells in safeguarding patients against distant metastasis, and suggests an important and quantifiable interplay between different immune cell populations in establishing long-term, metastasis-protective immunity.

24We assembled a multi-study microarray database of breast tumor expression profiles (n = 2,116) based on the Affymetrix U133 GeneChip microarray platform. The database encompasses 15 different breast cancer populations for which corresponding microarray data and clinical annotations were extracted from public data repositories including the Gene Expression Omnibus (National Center for Biotechnology Information, Bethesda, MD, USA), ArrayExpress (European Bioinformatics Institute, Hinxton, Cambridgeshire, UK) and caArray (National Cancer Institute, Bethesda, MD, USA) or by direct communication with study authors. Study population details and literature references are presented in Additional file 1. Previously unpublished breast tumor profiles from Belgium, England and Singapore have been deposited in the Gene Expression Omnibus and are accessible through GEO Series [GEO:GSE45255]. Raw array data (CEL files) were pre-processed and normalized using the R software package and library files provided by the Bioconductor project. In order to preserve a consistent normalization strategy across all study populations, raw data were MAS5.0 normalized on individual study populations using the justMAS function in the simpleaffy library from Bioconductor (no background correction, mean target intensity of 600). The specific array platforms employed were the HG-U133A, HG-U133 PLUS 2.0 and HG-U133A2 gene chips. To ensure equal information content from each chip type, only probe sets common to all chip types were utilized in subsequent analysis. This resulted in the use of 22,268 probe sets that were common to all microarrays in all study populations. Cross-population batch effects were corrected using the COMBAT empirical Bayes method. Of the initial 2,116 tumor profiles, 2,034 profiles represent primary invasive breast tumors sampled at the time of surgical resection, without exposure to neoadjuvant treatment. Of these, 1,954 cases were annotated with DMFS time and event. Other clinical annotation such as treatment type, ER status, nodal status, tumor size, histologic grade and patient age were available for the majority of cases. The pan-leukocyte expression profiles of Abbas et al. were downloaded from the Gene Expression Omnibus [GEO:GSE22886] and MAS5-normalized in the same fashion as the breast tumor datasets.

25Intrinsic breast cancer subtypes were assigned to samples using the Single Sample Predictor (SSP) algorithm described by Hu et al. and utilized by Fan et al.. Affymetrix probe sets were matched to the genes comprising the SSP centroids using UniGene annotation. Prior to batch-correction, the expression data for each gene were mean centered, and Spearman correlation was used to find the centroid most closely associated with each tumor sample. In cases where a correlation greater than 0.1 was not achieved with at least one centroid, a subtype was not assigned to that sample (n = 92 cases). Tumors representing the CL subtype were identified using the methods of Prat et al. and the supplementary information from. Briefly, CL centroids were generated using the Prat et al. microarray data set deposited in [GEO:GSE18229], with breast tumor samples assigned to the closest subtype centroid based on Euclidean distance.

26The 1,954 survival-annotated cases were dichotomized into training and testing sets comprising 977 cases each. Cases were iteratively randomized to two groups with monitoring of intergroup survival rates (log-rank test) and standardized differences for the variables: DMFS time, DMFS event, original study population, intrinsic subtype and ER status. The first randomization to achieve the following criteria was selected: log-rank test P-value for survival difference > 0.99, and < 10% standardized difference for each of the listed variables.

27Associations between gene expression and patient survival (DMFS) were assessed by Cox proportional hazards regression (likelihood ratio test) using the R survival package or by the Kaplan-Meier method (log-rank test, SigmaPlot 11.0). DMFS was defined as the absence of clinically confirmed distant relapse. Data were censored on date of last follow-up (disease-free), date of diagnosis of a second primary cancer, and local or regional relapse without evidence of distant recurrence. To minimize population-specific bias owing to patient follow-up duration, survival analyses were delimited to a 10-year window of patient follow-up. (Notably, only 1.5% of cases were annotated for post 10-year recurrence, and no convergence of survival curves post 10 years was observed.) Probe sets (genes) with likelihood ratio test P-values less than 0.01 and FDR-adjusted q-values less than 0.10 were identified and selected for supervised hierarchical cluster analysis. Univariable, multivariable and Kaplan-Meier analyses were performed in SigmaPlot 11.0. For statistical analyses low, intermediate and high tertile metagene designations were coded as 1, 2 and 3, respectively. All statistical tests were two-sided.

28Supervised hierarchical clustering and heatmap visualization of breast tumor and leukocyte gene expression profiles were conducted using Eisen's Cluster (v2.11) and TreeView (v1.60) software. Briefly, normalized log2 expression data were mean centered (on genes only), and genes and tumors were hierarchically clustered by average linkage using uncentered Pearson correlation as the distance metric. Clustered data were visualized by TreeView using default color saturation settings.

29Selection and construction of the immune metagenes was carried out essentially as previously described. Briefly, nested gene correlation structure was determined by the Pearson distance metric and average linkage clustering. The selected subclusters were then defined as the nested branches that generated an approximate average correlation of 0.6. This threshold was chosen to satisfy two primary goals: selection of genes with relatively high magnitude of correlation such that their correlation could not be considered a chance event; and selection of a reasonable number of genes (at least tens of genes) suitable for (GO) enrichment analysis. The metagene value for a given tumor was defined by averaging the signal intensities of the genes comprising each subcluster. In cases where two or more probe sets corresponded to the same gene identity within a subcluster, these probe sets were first averaged together, prior to cross-gene averaging, to guard against overrepresentation of any one gene with respect to its contribution to the metagene value. The metagene value (average signal intensity) thresholds defined by tertile cut-points in the training set are listed in Additional file 11.

30Thirty FFPE tissue blocks from surgically resected primary ER+ breast tumors were selected from a biobank of primary breast tumors representing recurrent cases and controls (that is, cancers that did not metastasize) maintained at Aarhus University Hospital, Denmark (SJHD and TLL). The blocks derived from patients with centrally confirmed ER+ stage II breast cancer diagnosed in 2000 or 2001. From each block, three 10-micron sections were cut to slides, the first of which was hematoxylin and eosin (H&E) stained to identify regions of highest tumor cellularity, which were then demarcated on the unstained slides. The H&E-stained slides were subsequently examined by a tumor pathologist at Wake Forest School of Medicine, NC, USA (MCW) for estimation of total leukocyte infiltrate within the demarcated regions of high tumor cellularity. Two slides were rejected based on evidence of poor preservation. For the remaining 28 samples, total leukocyte content was scored into four categories defined as low/absent, intermediate-low, intermediate-high, and high (irrespective of spatial considerations such as intra-tumoral or intra-stromal associations). As proof of principle and to test the technology platform, three genes were selected from each of two of the immune metagenes for analysis by Panomics QuantiGene Plex 2.0 (Affymetrix, Santa Clara, CA, USA)). Genes that displayed strong correlation with their cognate metagene and, simultaneously, large dynamic range of gene expression by microarray were selected for this purpose. Specifically, probe sets were designed to detect IGKC, IGLL5 and IGHA1 (the B/P metagene) and LCK, CD3E and CD27 (the T/NK metagene). Four housekeeping genes (ACTB, GAPDH, ACTG1, and EIF4G2) were included for normalization. Normalized expression ratios were generated by dividing the background-subtracted expression values by the geometric mean of the housekeeping genes. Lastly, the expression ratios were averaged to generate metagene values, and these values were correlated to total leukocyte abundance by Pearson correlation.

31The DAVID Bioinformatics Resource (Database for Annotation, Visualization and Integrated Discovery) version 6.7, sponsored by the National Institute of Allergy and Infectious Diseases, National Institutes of Health, was used to investigate the statistical enrichment of biological terms and processes associated with the genes comprising the immune gene clusters. Briefly, Affymetrix probe set unique identifiers were imported into DAVID and the functional annotation tools were utilized as described.

32B/P: B cell/plasma cell; CL: Claudin-Low; DAVID: database for annotation, visualization and integrated discovery; DMFS: distant metastasis-free survival; ER: estrogen receptor; FDR: false discovery rate; FFPE: formalin-fixed paraffin-embedded; GO: Gene Ontology; H&E: hematoxylin and eosin; HER2-E: and human epidermal growth factor receptor 2-enriched; HLA: human leukocyte antigen; LN: lymph node; LumA: luminal A; LumB: luminal B; M/D: monocyte/dendritic cell; MHC: major histocompatibility; PL, I, H: low, intermediate, or high proliferation tertile; T/NK: T cell/natural killer cell.

33The authors declare that they have no competing interests.

34SN and LDM conceived of the project and analytical strategy. MAB and LDM qualified and curated the breast cancer microarray datasets. MAB processed the microarray data (for example, normalization, batch-correction) and assigned cases to breast cancer subtypes. JWC contributed to data curation, quality control and subset analysis. JB and CS contributed microarray and clinical data, and contributed to manuscript development. TLL, SJHD and MCW provided paraffin blocks and pathological guidance, and contributed to methods development. SN, JR, PD, JPV and LDM vetted content and wrote the manuscript. All authors read and approved the final manuscript.