Classification of Breast Cancer Using Genetic Algorithms and Tissue Microarrays

Purpose: A multitude of breast cancer mRNA profiling studies has stratified breast cancer and defined gene sets that correlate with outcome. However, the number of genes used to predict patient outcome or define tumor subtypes by RNA expression studies is variable, nonoverlapping, and generally requires specialized technologies that are beyond those used in the routine pathology laboratory. It would be ideal if the familiarity and streamlined nature of immunohistochemistry could be combined with the rigorously quantitative and highly specific properties of nucleic acid–based analysis to predict patient outcome.

Experimental Design: We have used AQUA-based objective quantitative analysis of tissue microarrays toward the goal of discovery of a minimal number of markers with maximal prognostic or predictive value that can be applied to the conventional formalin-fixed, paraffin-embedded tissue section.

Results: The minimal discovered multiplexed set of tissue biomarkers was GATA3, NAT1, and estrogen receptor. Genetic algorithms were then applied after division of our cohort into a training set of 223 breast cancer patients to discover a prospectively applicable solution that can define a subset of patients with 5-year survival of 96%. This algorithm was then validated on an internal validation set (n = 223, 5-year survival = 95.8%) and further validated on an independent cohort from Sweden, which showed 5-year survival of 92.7% (n = 149).

Conclusions: With further validation, this test has both the familiarity and specificity for widespread use in management of breast cancer. More generally, this work illustrates the potential for multiplexed biomarker discovery on the tissue microarray platform.

A multitude of breast cancer mRNA profiling studies has stratified breast cancer and defined gene sets that correlate with outcome (1). These studies have resulted in plans for prospective application of nucleic acid–based tests to select patients that do not need further therapy after their primary resection (2, 3). However, the number of genes used to predict patient outcome or define tumor subtypes by RNA expression studies is variable, nonoverlapping, and generally requires specialized technologies. Immunohistochemical studies can be done with many fewer markers (4, 5), but suffer from the inherent flaw of subjective analysis and variable reproducibility. Thus, it would be ideal if the familiarity and streamlined nature of immunohistochemistry could be combined with the rigorously quantitative and highly specific properties of nucleic acid–based analysis to predict patient outcome. The AQUA method has the potential to combine the best properties of each of these assays. It is similar to immunohistochemistry in its format and slide preparation, but it allows rigorously quantitative and reproducible assessment of protein expression. These results require a specialized slide-reading platform that is commercially available from HistoRx (New Haven, CT).

Tissue microarrays provide a method for high-throughput protein expression analysis of large cohorts of cancer patients on a single slide (6) that results in standardization of many variables and capacity for embedded controls (reviewed in ref. 7). This technology has been slow to be used for discovery because the analysis of these arrays has generally been subjective, thus invalidating many of the algorithms used for discovery in nucleic acid array experiments. The recently developed automated quantitative analysis technology allows for immunofluorescent quantification reflecting the biological continuum of protein expression (8). Studies from our laboratory have shown that AQUA scores are highly correlated to protein content as measured by ELISA (9).¹

Thus, the combination of tissue microarrays and AQUA provides follow-up and extension of technology for cDNA microarray-based target discovery by allowing an in situ protein assay of markers of interest on large cohorts of tumors with the inclusion of spatial subcellular localization information and multiplexed analysis.

The oldest and most valuable prognostic and predictive marker in breast cancer is the estrogen receptor (ER). Estrogen not only regulates differentiation and growth in the mammary gland, but it is important in breast cancer development and progression. A multitude of studies have investigated the relationship between gene expression in ER-positive and ER-negative breast cancers to generate a list of genes that are able to discriminate between the two classes. There is overlap between ER status discriminating genes and estrogen-responsive genes, but the sets are substantially different. These studies, and the results of the subtype analyses (10–12), were the basis for our selection of 35 ER and estrogen-related markers as a starting point for the tissue microarray (TMA) marker discovery studies toward classification and prognostication of breast cancer (further justification of marker selection is discussed in the Materials and Methods section and in Table 1).

Tissue microarray construction. Formalin-fixed, paraffin-embedded breast cancer tumors from the archives of the Yale University Department of Pathology were used in the construction of the tissue microarrays. Briefly, representative areas of invasive ductal carcinoma were selected from whole tissue H&E sections by pathologists (R.L.C., D.L.R.). Cores sized 0.6 mm in diameter were spaced 0.8 mm apart in a grid layout using a manual Tissue Microarrayer (Beecher Instruments, Silver Spring, MD). The resulting tissue microarray blocks were cut to 5-μm sections with a microtome, placed on slides with an adhesive tape-transfer method (Instrumedics, Inc., Hackensack, NJ), and UV cross-linked.

Cohort information. The Yale breast cancer cohort consists of 688 samples of invasive ductal carcinoma selected from the Yale University Department of Pathology archives as available from 1961 to 1983. The mean follow-up time of this cohort is 12.8 years and the mean age of diagnosis is 58.1 years. The median follow-up time is 8.9 years and the median age of diagnosis is 58.0 years. Of the 656 patients with follow-up time, 328 were censored at 20 years and 276 were uncensored at 20 years. Of the 328 censored patients, their median follow-up was 21.4 years, with the minimum at 4.2 months. This cohort contains approximately half node-positive specimens and half node-negative specimens, which have been used previously in other studies from our laboratory (9, 13–18). Complete clinical treatment information was not available for this cohort; however, some limited information was available for the node-positive patients. Most node-negative patients were treated with local radiation, and none received Herceptin. About 15% of the node-positive patients were treated with chemotherapy (primarily Adriamycin, cytoxan, and 5-fluorouracil), and subsequently ∼27% received tamoxifen (post-1978). The node-negative patients were largely treated only by surgical resection.

The Swedish breast cancer cohort consists of 564 premenopausal patients with invasive breast cancer, Unio Internationale Contra Cancrum stage II, enrolled in a randomized prospective clinical trial from 1986 to 1991 at two centers in Sweden (19). The patients were randomized to receive adjuvant tamoxifen (n = 276) for 2 years or no adjuvant treatment (n = 288) after the primary surgical treatment. Locoregional radiotherapy was delivered to all node-positive patients. Less than 2% of all of the patients received other adjuvant treatment. The median follow-up time was 13.9 years. The arms were well balanced, both having a median age of 45 years old. Between 27% and 30% of each arm showed no nodal metastasis at time of diagnosis.

The analysis of the 35 markers was all done on a 250-case cohort (9) of half node-positive and half node-negative specimens; the associations with clinical and pathologic variables with patient outcome are displayed in Table S1A. There were 105 classified as ER negative (44.9%) and 129 cases classified as ER positive (55.1%) for those tumors in which the clinical ER status (as determined by immunohistochemical analysis) were available. Previous analysis from our laboratory has shown a high correlation between our AQUA scores and pathologist-designated ER status (8). A larger cohort, used for further analysis of ER, GATA3, and NAT1, consists of a total of 675 tumors, including cores from the 250-case tissue microarray tumors, with the additional samples consisting of 193 ER-negative cases (48.3%) and 207 ER-positive cases (51.8%) for those cases with available ER status information. The associations with clinical and pathologic variables with patient outcome for the entire cohort are displayed in Supplementary Table S1B. This entire cohort contains approximately half node-positive specimens and half node-negative specimens, and has also been used previously in other studies (8, 13, 14, 16, 18, 20, 21) and is described in detail in Supplementary Data. The cohort used for validation is from the control arm of a study from a randomized trial of tamoxifen conducted in Sweden from 1986 to 1991 (19). Details of this cohort are described in Supplementary Data. In all cohorts, tissue microarray spots without sufficient breast tumor epithelium were excluded from the analysis.

Description of selection of genes. We mined the published results (and where possible, the full data sets) of 30 breast cancer RNA expression profiling studies to identify an overall list of candidate biomarkers to screen using tissue microarrays and AQUA analysis. The list of biomarkers were prioritized based on a combination of the following information: (a) appearance in multiple independent breast cancer profiling studies (references listed in Table 1); (b) classification strength/ranking of biomarkers in each study; (c) mapping to regions of breast cancer genetic alterations, amplification, or deletions based on a breast cancer cytogenetic meta-analysis of 15 papers that studied an overall total of 644 breast cancer tumors and/or cell lines (22–36); and (d) whether an antibody was available (either commercially or from individual laboratories) with validation in the literature by Western blotting and/or successful immunohistochemical analyses. This list was then further focused to only those who appeared in estrogen-related studies and that were successfully used for immunofluorescence, resulting in 35 markers successfully screened on tissue microarrays by AQUA in this study.

Immunofluorescence staining. The tissue microarrays were deparaffinized by two xylene rinses of 30 minutes each followed by two rinses with 100% ethanol for 1 minute each and a rinse in water. Antigen retrieval was done by boiling the slides in a pressure cooker in a sodium citrate buffer at a pH of 6.0. After rinsing briefly in 1× TBS, a 30-minute incubation with 2.5% hydrogen peroxide/methanol was used to block endogenous peroxidases. To reduce nonspecific background staining, slides were incubated with 0.3% bovine serum albumin/1× TBS for 1 hour at room temperature, followed by a series of 2-minute rinses in 1× TBS, 1× TBS/0.01% Triton, 1× TBS (TBS washes). Slides were incubated overnight at 4°C with either a monoclonal mouse anticytokeratin antibody (clone AE1/AE3; DAKO, Carpinteria, CA; 1:100) when using a rabbit or goat target antibody, or a rabbit anti-cytokeratin antibody (1:100; DAKO) when using a mouse target antibody. The sources, dilutions, and incubation times of the target antibodies are listed in Supplementary Table S1. Antibodies were selected that had either been generated and validated by individual laboratories, or were well-characterized antibodies frequently used for immunohistochemistry in the literature. Slides were washed in 1× TBS rinses and incubated with secondary antibodies for 1 hour at room temperature as follows: Alexa 488 goat anti-mouse or Alexa 488 goat anti-rabbit (1:100, Molecular Probes, Eugene, OR) for detecting cytokeratin, and species-specific horseradish peroxidase with a dextran-polymer backbone (Envision, DAKO) along with 6-diamidino-2-phenylindole (1:100, DAKO) for visualization of nuclei. For goat primary antibody detection, slides were incubated for 1 hour with biotinylated anti-goat (1:200, Vector, Burlingame, CA) and Cy-2-donkey anti-mouse (1:50, The Jackson Laboratory, Bar Harbor, ME). This was followed by a 1-hour incubation with streptavidin horseradish peroxidase (1:200, Perkin-Elmer, Wellesley, MA). All slides were washed with TBS rinses followed by a 10-minute incubation with Cy-5 tyramide (1:50 dilution in Amplification Diluent, Perkin-Elmer). The slides were mounted in 0.6% n-propyl gallate (an antifade mounting medium) and coverslipped. Examples of immunofluorescent staining are shown in Fig. 1.

AQUA analysis. AQUA software linked to a fluorescence microscopy system allowed for quantification of the protein of interest within the tumor region of each tissue microarray core, as described in detail by Camp et al. (8). Image acquisition begins with a low-resolution (64 × 64 pixels) image capture with a ×4 objective. The rows and columns of the tissue microarray are defined to form a grid on which the histospots are placed. Monochromatic, high-resolution (1,024 × 1,024 pixels), in-focus and out-of-focus images are then captured for each relevant wavelength for each histospot with a ×10 objective. The pan-cytokeratin antibody separates the epithelial breast cancer tumor component by binary masking from the surrounding stroma, and other fluorescent tags designate subcellular compartments (i.e., 6-diamidino-2-phenylindole for nuclei). Coalescence of cytokeratin staining was used to define the nonnuclear compartment for analysis of cytoplasmic and/or membranous staining. The images of the markers were captured at the Cy-5 wavelength because it is outside the range of tissue autofluorescence. Two algorithms, rapid exponential subtraction algorithm and pixel-based locale assignment for compartmentalization of expression, were then used for AQUA analysis. RESA improves subcellular compartment assignment by taking into account overlapping of cells due to the thickness of the tissue microarray section on the slide by subtracting the in-focus information from the out-of-focus information. The second algorithm, pixel-based locale assignment for compartmentalization of expression, assigns target staining in Cy-5 to the predefined subcellular compartments and quantifies its intensity within each compartment to give a quantitative measurement of the expression level for each histospot. The resulting AQUA scores are the measurements of the biomarker pixel intensity within a compartment divided by the total area of tumor (to normalize for differences in tumor area in each spot), and were used as their raw values in the statistical analyses described below, unless otherwise noted. Tumor samples with more than one AQUA score for a given case were averaged.

Unsupervised hierarchical clustering.Z-score transformation (37) was used to normalize between experiments to perform clustering analyses with this formula: [(AQUA score) − (mean of AQUA scores on tissue microarray X)] / SD, where X is a given tissue microarray experiment assessed for a single marker. This has previously been used for analysis of AQUA scores (38).¹ CLUSTER (39) was used to perform unsupervised average linkage hierarchical clustering on unweighted Z-score–transformed AQUA data and visualized with TREEVIEW (39).

Genetic algorithm for generation of a multiplex marker assay. A genetic algorithm was used to develop a multiplex marker assay for separation of prognostic groups based on expression of ER, GATA3, and NAT. A training set (n = 223) and validation set (n = 223) were randomly assigned to generate two groupings with equivalent baseline survival curves for all patients with values for all three markers and both censor and follow-up information. Data for each marker was normalized to control for differences in exposure time based on linear regression of redundant cases.

The genetic algorithm was used to determine the optimal cutpoints for the three selected biomarkers (40, 41). We generated software to produce random “solutions,” each containing three algorithms (one each for ER, GATA3, and NAT1), where each algorithm represents an equation that is a true/false statement. For each marker, a value is either greater than or less than or equal to a randomly selected cutpoint. The solutions consisted of three if/then statements that can be calculated as true or false. Each tumor was then classified on the basis of the number of true conditions where they are given one point for each criteria that is true, or zero points when false. Thus, each tumor has a final genetic algorithm solution score between 0 and 3. These scores are used to group the tumors as follows; group A (3 points), group B (2 points), group C (1 point), and group D (zero points). Then, the prognostic value of a solution was determined by assessing the risk ratio for disease-specific death of the D group (high risk) compared with the A group (low risk). The process is then repeated with new random cutpoint values for each marker and the new solution is compared with the previous solution in an iterative manner—selecting the solution with the highest risk ratio and discarding other solutions until the model converged onto a maximum predictive value (over 20 million iterations). As a further limitation, solutions were only evaluated if they resulted in at least 10% of the training set being assigned into the D group and at least 10% assigned into the A group (to prevent finding groups of less than 20 patients each in the low and high risk groups). The genetic algorithm resulted in a model with the best solution (i.e., highest relative risk) on the training set for stratifying 0-scoring tumors from 3-scoring tumors as follows: ER > 4.6, GATA3 > 33.8, NAT1 > 18.6. The tumors in the training set were assigned into one of the four groups: high risk (score of 0, group D), intermediate risk 1 (score = 1, group C), intermediate risk 2 (score = 2, group B), or low risk (score = 3, group A). This model was applied to the validation set.

Statistical analyses. Statview 5.0.1 (SAS Institute, Inc., Cary, NC) was used for the Cox proportional hazards analysis. Clusters from the unsupervised hierarchical clustering analyses and groups from the multiplex marker assay were assessed for relationship to patient prognosis by Kaplan-Meier survival curves with their significance analyzed by the Mantel-Cox log-rank test with Statview 5.0.1 or SPSS v12 (SPSS, Inc., Chicago, IL).

We obtained antibodies to 35 markers that were either well-characterized commercial antibodies or well-studied antibodies from individual laboratories (Table 1). We then assessed a 250-case breast cancer cohort by AQUA to quantify the expression of each marker in the epithelial tumor regions of the tissue microarray histospots. The resulting scores are measurements of the biomarker pixel intensity within a compartment divided by the total area of epithelium to account for differences in epithelial area in each spot (see ref. 8 for complete explanation of AQUA technology). The generation of continuous data for each marker allowed us to objectively assess protein expression and to use the data in rigorous algorithmic analyses.

To further analyze each marker, we used the survival and disease-specific death information for the cohort to investigate the relationship of each marker with prognosis by univariate Cox proportional hazards (Table 2). Six of the 35 markers were significant when assessed by their continuous protein expression levels for relationship with 5-year disease-specific survival, including BCL2 (P = 0.0400), COX6C (P = 0.0049), ER (P = 0.0115), GATA3 (P = 0.0011), HER2 (P = 0.0048), and NAT1 (P = 0.0206). Examples of immunofluorescence staining are shown in Fig. 1 for high protein expression of ER, NAT1, and GATA3 (Fig. 1A, B, and C, respectively), and for very low expression of ER, NAT1, and GATA3 (Fig. 1C, D, and E, respectively).

Although the individual prognostic value of markers is important, the classification strength from array studies derives from the power of multiplexing. We began by using the hierarchical clustering and visualization of a heat map to provide a context for exploration of the biological coexpression. Unsupervised average linkage clustering was applied to this set of markers as shown in the heat map in Fig. 2A. The cases shown include those with data for at least 80% of the markers to allow inclusion of a larger number of samples due to missing data. The total number of tumors clustered is 151 due to missing values. ER is in a small cluster with two other markers, GATA3 and NAT1. The larger cluster containing the GATA3, NAT1, and ER grouping also includes BCL2, HSP27, SLC9A3R1, and IGFBP4. When looking at the results of linear regression analysis for these markers individually with ER protein expression levels, they had statistically significant direct relationships (i.e., higher ER levels corresponding to higher marker levels) as follows: BCL2 (R = 0.423, P < 0.0001), GATA3 (R = 0.497, P < 0.0001), HSP27 (R = 0.327, P < 0.0001), IGFBP4 (R = 0.203, P = 0.0055), NAT1 (R = 0.487, P < 0.0001), and SLC9A3R1 (R = 0.187, P = 0.0182). There was a noticeable grouping of markers (CDH3, KRT7, GGH, and HER2) away from the rest of the markers in the heat map that generally have indirect relationships with ER levels (i.e., higher ER levels correspond to lower levels). Also to be noted is the close relationship in the clustering analysis of keratin 8 and 18, which are coexpressed (42).

The cluster with the highest ER levels, cluster C, had the best prognosis. Cluster A, with the lowest ER levels, had a poor prognosis (87.5% and 60.0% 5-year survival, respectively; Fig. 2A and B). When comparing the cluster groupings with traditional immunohistochemically determined ER status (>10% nuclei staining), there were differences in ER status positivity between the groupings. Cluster C consists of 81% ER-positive tumors, cluster B with 56% ER-positive tumors, cluster D with 53% ER-positive tumors, and cluster A with 12% ER-positive tumors. This analysis shows that multiplex expression analysis can reveal prognostic classes that are substantially different than that obtained with current immunohistochemical methods.

As the goal of this discovery was to find a minimal number of markers with maximal prognostic value, we focused on the small grouping of two markers, GATA3 and NAT1, with highest correlation with ER protein expression (highlighted in Fig. 2A). As shown in the univariate Cox proportional hazards analyses in Table 2, higher levels of ER, GATA3, and NAT1 are all individually related to better patient prognosis.

To more rigorously evaluate these three markers, AQUA analysis for ER, GATA3, and NAT1 was done on an expanded set of breast cancer tumors on tissue microarrays (cohort details in Supplementary Table S1B). In this cohort, pairwise evaluation of these markers by Cox multivariate analysis methods shows that each maintains independent predictive value. The results of unsupervised hierarchical clustering of ER, GATA3, and NAT1 for all of the breast cancer tumors (with data for all three markers) defined two main clusters, cluster 1 with different subsets of coexpression of NAT1, GATA3, and ER, and cluster 2 with predominantly low levels of the three markers (Fig. 2C). These two clusters had significant differences in prognosis as shown in Fig. 2D with 84.6% 5-year survival in cluster 1 (log-rank P < 0.0001) versus 65.9% in cluster 2. Although cluster 2 shows remarkably homogeneity in expression of the three markers, cluster 1 illustrates that although having overall higher levels ER, NAT1, and GATA3, there are groupings of marker coexpression that subclassify patients with good outcome (Fig. 2E). Specifically, one cluster (1B) shows a 97% 5-year survival compared with another (1C) with only 80%. Although the unsupervised clustering can define classes that represent dramatic differences in outcome, the method does not lend itself to prospective assignment of class, nor does it take advantage of outcome-based information in class discovery. We sought another method for analysis of this three marker set for potential prospective use.

Genetic algorithm methods have been used extensively in other fields (for review, see ref. 40) to discover a set of mathematical functions that best define properties of a population and more recently have been used on gene expression data (41). They are less popular in array analysis due to the large number of genes and, hence, lower likelihood of convergence. But the limited number of protein biomarkers abstracted from the clustering analysis lends itself well to the genetic algorithm approach. Iterative analysis of algorithms (fitness tests) on a training set results in convergence to optimal predictive values that can then be tested on a validation set. Thus, to optimally select patient subsets with the best prognosis by a multiplex assay, we divided our entire cohort into two randomly assigned, training and validation sets with equal baseline survival. We used a genetic algorithm to define a series of risk groups for disease-specific death at 5 years. Each algorithm in this case represents an equation that is a true/false statement, such that the marker value is either greater than or less than or equal to a randomly selected cutpoint. A “solution” thereby consists of three if/then statements that can be calculated as true or false. Each tumor is then assessed one point for each criterion that it meets when true, or no points when false. Tumors can thereby receive a 0, 1, 2, or a 3 cumulative score. These scores are referred to as A (3 points), B (2 points), C (1 point), and D (zero points). The predictive value of a solution was determined by assessing the risk ratio of the D group (high risk) compared with the A group (low risk). The solution that converged after over 10,000,000 generations (iterations) with the best fit (highest relative risk between good and poor prognosis groups) was selected for validation. The solution was ER > 4.6, GATA3 > 33.8, and NAT1 > 18.6, where the units are normalized AQUA scores.

The application of this solution to the training set showed that 27 patients were assigned to the good prognosis (low risk), with one event resulting in 96.3% survival at 5 years, and 23 patients assigned to the poor prognosis group (high risk) with eight events (34.8% survival at 5 years). We were also able to assign patients into two intermediate risk groups, with 78.6% and 62.3% 5-year disease-specific survival. The Kaplan-Meier survival curves of the four groups for the training set are shown in Fig. 3A (log-rank P < 0.0001). We then applied this algorithm to the internal validation set, which resulted in the Kaplan-Meier survival curves shown in Fig. 3B (log-rank P = 0.0073). In the validation set, the 24 patients in the low-risk group had only one event (95.8% 5-year survival), and the 24 patients in the high-risk group had 58.3% 5-year survival, showing again that our model was able to stratify patients into different prognostic groupings on the basis of a prospectively assignable score using only three protein expression measurements; this grouping remains independent in a multivariate Cox model, including age at diagnosis; nuclear grade; tumor size; nodal status; clinical immunohistochemical ER status; and the individual AQUA scores for ER, GATA3, and NAT1 (Table 3).

For external validation of these findings, we repeated this assay on a cohort of patients from a randomized clinical trial of tamoxifen conducted in Sweden from 1986 to 1991 (19). A tissue microarray from these tumors were analyzed for ER, GATA3, and NAT1, normalized for exposure times to the Yale cohort, and then grouped as defined by the genetic algorithm above. The data from 149 patients from the control arm of this study show a 92.7% 5-year survival in group A compared with 58.8% 5-year survival in the high risk (group D). Like the Yale cohort, this group of patients was treated only with surgery. Unlike the Yale cohort, this group was limited to premenopausal patients with low stage (I or II) breast cancer representing a more challenging cohort for classification.

The ability to multiplex markers allows for greater complexity in the assessment of multiple biomarkers that can contribute to predicting patient outcome. Several breast cancer studies have used such a strategy, such as the ratio between HOXB13 and IL17BR gene expression to predict response to adjuvant tamoxifen therapy (3) or the use of the 21 reverse transcription-PCR based test described by Paik et al. (2). There are several commercially available tests based on gene prognosis signatures developed from expression profiling studies (e.g., refs. 2, 43, 44), which could be valuable for subsets of breast cancer patients. This study extends the multiplexing concept to protein expression analysis by quantitative extension of the familiar immunohistochemistry platform.

Several recent breast cancer tissue microarray studies have used hierarchical clustering to immunoprofile breast cancer cohorts with semiquantitative scoring methods to identify prognostic subgroups. These publications include classification by binarized ER, HER2, COX-2, p53, and vascular endothelial growth factor expression (45); analysis of breast cancer expression of neuroendocrine markers (46); definition of the basal subtype (5); confirmation of cDNA expression data (47); and several others focused on prognosis prediction (48–50). The most recent work using this strategy is by Ring et al. (51) and, like this work, it results in patient stratification on the basis of outcome. Some statisticians argue that this type of clustering is invalid because the clustering is done on ordinal, rather than continuous, data. This study, as well as another recent AQUA-based study (from our laboratory) used hierarchical clustering on continuous data to identify a subset of patients with poor prognosis (52).

Multiplexing, as described in this work, is a combination of true multiplexing and serial section multiplexing. True multiplexing is when multiple markers are measured on the same tissue spot. Although we are multiplexing by measuring keratin, 6-diamidino-2-phenylindole, and our target all on the same tissue spot, we are only quantitatively measuring one marker (the target) per tissue section. We are able to do true multiplexing with multiple targets, but currently only two targets can be assessed per tissue spot. We clustered our data based on a serial section multiplexing. Because the data is quantitative and continuous, we believe it is statistically valid to do this type of clustering analyses.

This analysis was restricted only to genes with available antibodies and could thus be confounded by this selection. However, the ultimate goal of this type of effort is to discover the smallest set of biomarkers with the maximum prognostic value. Toward the same goal, other groups have come up with similarly small, but unique sets of antibodies that are highly prognostic (51). We anticipate that future studies with other antibodies may suggest the addition of new markers or even the removal of some of the three markers identified in this study. Thus, it is possible, or even likely, that the ultimate most prognostic algorithm will extend or otherwise change the three-marker set defined in this work. Furthermore, without any antibodies available or generated, semiquantitative analysis of in situ hybridization results could provide additional information on highly interesting targets as supplement to this type of study.

A limitation of this study is that the predictive value of this test is only valid at 5 years. When this same set of markers was tested for 15 years, the significance is lost. However, this is most likely due to the fact that the genetic algorithm is derived using data for 5-year disease-specific survival. The use of this or any test to exclude patients from chemotherapy, as proposed for the Oncotype Dx test (2), requires assessment of longer-term outcome. Work is under way on this data set to discover a small set of markers that predict low recurrence frequency over a long time period. However, the key value of an assay of this nature is to find the low-risk patients that are highly likely to recur. Table 3 shows that group D patients have a relative risk of nearly 14. Furthermore, the external validation set of young node-negative patients in the external validation cohort shows that group D can define a subset of nearly a quarter of the patient population that has a 40% chance of recurrence. With further validation, this technology has the potential to translate into a useful clinical test. It could be used to select a subset of breast cancer patients with a subclass of tumor that is highly likely to recur and thus patients could be strongly encouraged to undergo adjuvant therapy, even if tumor stage and other clinical variables may suggest otherwise. Further analysis is under way with other markers to increase the classificatory power of this type of assay using the AQUA platform.

We thank Laura Dillon, Thomas D'Aquila, Lori Charette, Matthew Neopolitano, and Jung Kang for their technical assistance in this project.