Recent translational research: microarray expression profiling of breast cancer – beyond classification and prognostic markers?

A fundamental variable to consider is the quality of the data that can be obtained from microarray expression profiling of complex, heterogeneous epithelial tumors. Specifically, are the data sufficiently quantitative to allow for the recognition of coordinated patterns of gene expression indicative of a particular signaling pathway? To determine what we might expect under the best circumstance, we examine selected genes whose expression should be particularly well coordinated in breast cancer cells.

ERBB2 is amplified and pathologically overexpressed in about 25–30% of breast cancers [9] along with the neighboring gene GRB7 [10]. The log ratios or intensity values have been downloaded for these two probes from each of four publicly available primary breast cancer microarray data sets [3–5, 11]. High positive correlation coefficients for ERBB2 and GRB7 co-expression ranging from 0.633 to 0.910 (Table 1 and Fig. 2) were found in all four data sets. For each study, the corresponding graph in Fig. 2 provides a good indication of which tumors are amplified at the ERBB2 locus.

	Data
Genes compared a	van't Veer and colleagues b (n = 98)	Sørlie and colleagues c (n = 122)	Sotiriou and colleagues d (n = 99)	West and colleagues e (n = 49)
ERBB2 and GRB7	0.910	0.745	0.633	0.887
ERBB2 and STARD3	0.791	0.705	0.828	0.895
KRT5 and KRT14	0.745	N/A f	0.546	N/A
KRT5 and KRT17	0.837	0.784	N/A	N/A
KRT8 and KRT18	0.886	N/A	0.712	0.770
ESR1 and LIV-1	0.728	0.670	0.674	N/A
ESR1 and TFF1	0.624	N/A	N/A	0.077
ESR1 and PGR	0.493 g	N/A	N/A	N/A h
ESR1 and AREG	0.500	N/A	N/A	N/A
ESR1 and TGFβ2	-0.212 i	N/A	-0.160	-0.091
ESR1 and NDRG1	-0.669	N/A	-0.015	-0.329

^a If more than one probe was present for a gene, the one with the highest correlation and the fewest missing data values was used. ^b Correlation coefficients calculated from log₁₀ intensity ratios from van't Veer and colleagues [11]. The 78 sporadic tumors and 20 BRCA1/BRCA2 tumors were used to compute the correlations. ^c Correlation coefficients calculated from log₂ intensity ratios from Sørlie and colleagues [4]. ^d Correlation coefficients calculated from log₂ intensity ratios from Sotiriou and colleagues [5]. ^e Correlation coefficients calculated from intensity values from West and colleagues [3]. The R ² value for the plot of ERBB2 and GRB7 in Fig. 2a is based on log₁₀(intensity). ^f N/A, data for one or both the two genes not included in the publicly distributed data set. ^g The strength of the positive correlation is reduced by the wide range of values for the progesterone receptor (PGR) in the strong estrogen receptor message (ESR1) subgroup (Fig. 5c). ^h PGR values for 39 of 49 samples were negative. ⁱ The strength of the negative correlation is reduced by the wide range of values for TGFβ2 in the weak ESR1 subgroup (Fig. 5e).

We also looked for co-expression of cytokeratins as another measure of data quality. Cytokeratins are abundant proteins that form the intermediate filaments of epithelial cells. The basic units of the fibers are heterodimers of one type I cytokeratin and one type II cytokeratin subfamily member [12], and distinct 'expression pairs' have been identified including KRT5/KRT14 and KRT8/ KRT18 [13]. Thus, these genes should show a high degree of co-expression. In every case, when the probes were present, the correlation coefficients were high for the co-expression of KRT8 with KRT18 and for the co-expression of KRT5 with KRT14 (Table 1). These correlation coefficients were in the similar range of high significance that was observed for genes co-amplified with ERBB2 (GRB7 and STARD3).

These and other examples confirm that microarray platforms have, in fact, generated high-quality gene expression data with a strong quantitative character for RNA isolated from human tumor samples. In general, the correlations and gene coverage were highest in the data of van't Veer and colleagues [11], which used a 60-mer oligonucleotide array platform representing approximately 25,000 genes. Also, these data had the fewest missing values. The reference in this study was a pool of RNA extracted from all 78 sporadic tumors. Interpreting the data intensity values in this case is simplified since it is intuitive to think of zero as the average expression of a particular gene in this breast cancer cohort [11]. We have focused on the van't Veer and colleagues data for tumor comparisons in the remaining discussions.

It has been often reported that the gene expression patterns associated with ER status in breast cancer are remarkably distinct and that the set of ER classifiers is comprised of up to several hundred genes [3, 11, 14–17]. We believe that many or even most of these 'ER predictors' primarily distinguish tumors according to cell-type origin (i.e. those tumors with predominantly 'luminal features' from those tumors with predominantly 'basal features' [1]) rather than according to ER regulation. To make this argument, the expression ratios of estrogen receptor message (ESR1) and ERBB2 from the van't Veer and colleagues data [11] were used to divide the 78 sporadic tumors into five groups. The samples are arranged from highest to lowest ESR1 level in Fig. 3a, with the ERBB2 tumors grouped separately. The ERBB2 tumors were identified by positive values for both ERBB2 and GRB7 (Fig. 3b). A sixth group was defined based on BRCA1/BRCA2 mutation status (Fig. 3).

Figure 3 highlights several important features of the data. The first is that ESR1 is expressed as a continuous variable whereas ERBB2 and GRB7 have essentially a binary expression pattern due to gene amplification. Second, none of the ERBB2 samples have above average (positive log ratios) values for ESR1. This is consistent with other larger data sets measuring ERBB2 and ER protein levels as continuous variables, where it has been suggested that low ER levels contribute to the reduced anti-estrogen sensitivity of ERBB2 amplified tumors [18]. Another feature is that the assignment of -0.5 as the cutoff for 'true' ER negativity results in 17/18 of the BRCA1 tumors being classified as ER-negative, as has been confirmed elsewhere [19]. Only the two BRCA2 (Fig. 3a) samples and a single BRCA1 tumor appear to express any ESR1.

We also classified the tumors from the breast cancer data sets [3–5] based on the level of ESR1 alone and compared these results with the subgroups generated by the various clustering methods (see Additional file 1). All of the tumors defined as basal using clustering methods were in the lowest ESR1 category, and nearly all luminal A or luminal 1 tumors were found in the highest ESR1 groups. It is clear that the tumor groups defined by the extremes of ESR1 expression are the most easily recognized and consistently observed subtypes of breast cancer. Tumors with mid-range ESR1 expression (luminal B/C, luminal 2) will require further analysis with larger numbers of tumors to define the dominant molecular signals driving their pathogenesis.

As previously reported, the signatures of luminal cell types versus basal cell types are dominant in the expression profiles obtained from primary breast cancer [1–5]. One explanation for these strikingly different patterns is that the cell type in which the oncogenic transformation took place is fundamentally different in these two tumor groups (Fig. 1). The use of even a small number of markers illustrates this distinction. We show the co-expression of the prototypical luminal cytokeratins (KRT8/KRT18) and the basal cytokeratins (KRT5/KRT14) from the van't Veer and colleagues study [11] in Fig. 4.

The samples clearly show coordinated expression of KRT8 and KRT18 (Fig. 4a), and a striking pattern is apparent. Four different groups comprise the majority of the samples having positive values for KRT8/KRT18: the ESR1 expressing groups (strong ESR1, moderate ESR1 and weak ESR1) and the ERBB2 amplified samples. It should be noted, however, that within these four groups the expression of these luminal cytokeratins does not correlate with the level of ESR1 expressed. The ESR1-negative tumors, including both the sporadic and BRCA1 mutant samples, stand out as having especially low levels of KRT8/KRT18.

The reverse is true for the expression of KRT5 and KRT14 (Fig. 4b), as well as for KRT5 and KRT17 (data not shown). It is the 'truly' ER-negative tumors that have high basal cytokeratin expression, and the four groups of luminal tumors cluster together in the negative region of the plot. There are hundreds of genes that divide these tumor samples into these two main groups (luminal and basal), and these genes are often a major component of the various ER discriminator gene sets [3, 11, 14–17]. However, the role of these genes in ER signaling and their regulation by either E₂ or the ER remain unconfirmed.

The ligand-dependent genomic action of the ER is relatively well understood and the in vitro analysis of E₂-responsive genes in breast cancer cells has been actively pursued for many years [20]. Part of the Cunliffe and colleagues study [7] was to characterize the dynamic transcriptional response of two different breast cancer cell lines (MCF-7, T47D) to 17β-estradiol and anti-estrogens (ICI 182,780 and 4-hydroxy tamoxifen) using a custom-made 10 K cDNA array. This resulted in the identification of 386 hormone-responsive genes.

The study by Frasor and colleagues [8] was undertaken to better understand the transcriptional activities of selective ER modulators in breast cancer cells. They report the transcriptional changes induced in MCF-7 cells by E₂, and classify the genes according to their response to the pure anti-estrogens (ICI 182,780) and trans-hydroxytamoxifen and raloxifene using the Affymetrix Hu95A array. Their analysis identified a highly focused E₂-responsive gene signature of 129 genes.

The issue of cross-platform consistency has been recently explored in detail [21, 22]. In one study, the correlation for gene expression data from breast cancer cell lines obtained using Affymetrix HG-U95v2 and a custom-made cDNA array was found to be in the range 0.66–0.76 [21]. The differences between platforms appear to result from errors in chip fabrication, ambiguities in gene annotation, specificity differences inherent in the hybridization of oligonucleotides versus cDNA clones and alternative methods for data filtering and normalization [21]. The authors found that the biological differences between the cell lines (e.g. BT-474 versus MCF-7) were more prominent than the variation between platforms. However, it appears probable that platform variability actually exceeds changes in expression induced by treating MCF-7 cells with E₂. In this light, the failure to find a robust consensus gene signature for the E₂ signaling pathway is not surprising.

Despite the disappointing overlap, some of the genes identified in these two studies are still likely to be valid targets of ER signaling. In order to investigate the feasibility of integrating experimentally generated pathway responses with data from primary tumors, we focused on a few specific genes identified in each of the in vitro studies. The expression values for each gene were plotted against the ESR1 level for the sporadic tumors (n = 78) in the data of van't Veer and colleagues (Fig. 5).

We compared LIV-1, AREG, TFF1 and PGR with the ESR1 expression levels in the van't Veer and colleagues data [11]. LIV-1 and AREG were identified as E₂-induced and ICI-repressed in the Frasor and colleagues study [8] and TFF1 was found to be E₂-induced and ICI-repressed in the Cunliffe and colleagues study [7]. PGR, the best known target of liganded ER, was not identified in either study. Each experimentally identified target gene showed a reasonably good correlation with ESR1 in the van't Veer and colleagues data (Table 1), and an interesting pattern is evident when plotted using the color-coded groups (Fig. 5a,5b,5c,5d). There is a trend for higher expression of these E₂-stimulated genes in the samples that express more ESR1: moderate ESR1 samples > weak ESR1 samples. We also examined the expression of several of the genes identified by Frasor and colleagues to be downregulated by E₂ and upon which ICI acted as an antagonist. Two of these genes, TGFβ₂ and NRDG1, have negative correlations to ESR1 in the data of van't Veer and colleagues [11], consistent with being targets of repression by ER signaling (Table 1 and Fig. 5e,5f).

The quantitative expression of ER has been shown to have clear clinical implications both in the adjuvant and the metastatic setting in terms of response to anti-estrogen therapy [23–25]. This is certainly consistent with the concentration-dependent occupancy of cis-regulatory sites being a fundamental aspect of DNA-binding transcription factors. It is interesting to note that about one-half of the tumors with the highest ESR1 expression do not show co-expression of PGR, AREG or TFF1 (Fig. 5b,5c,5d). The ER signaling pathway in these very highly ESR1 expressing samples may thus be fundamentally different. In a patient cohort very similar to this one (node-negative patients, younger than 60 years old at diagnosis), where ER protein was measured quantitatively, a paradoxical reduction in overall survival was associated with very high ER levels [26]. The clinical outcome data in the van't Veer and colleagues study also exhibits a trend towards a worse prognosis in the highest ESR1 subgroup (see Additional file 2).

The E₂ response in breast cancer cells would seem to provide a straightforward opportunity for microarrays to demonstrate their utility in identifying pathway signatures based on gene expression profiles. ER is a transcription factor with well-characterized cis-regulatory sites, and we have an abundance of hormone receptor agonists and antagonists at our disposal with which to modulate its activity. Decades of clinical and laboratory research with these compounds have validated the central importance of E₂ and ER in the pathogenesis of breast cancer. There is also strong evidence indicating that the currently employed microarray platforms can support the quantitative analysis of cell and tumor transcriptional programs.

So why is the global ER gene expression signature for E₂-responsive breast cancer cells still unclear? There are several factors in addition to the cross-platform comparison issues already discussed. One longstanding technical obstacle is the over-reliance on a single ER-positive breast cancer cell line (MCF-7) or on only a few ER-positive breast cancer cell lines (ZR-75-1, T47D) as experimental models. We have suggested that the true ER signature may be much smaller than originally proposed due to a failure to appreciate the close association between luminal cell differentiation and ER activity. The data in Figure 4 suggest that this association may be circumstantial. Alternatively, it is possible that breast cancers arise from a transforming event in multipotent progenitor cells and that acquired alterations (e.g. high autocrine ER stimulation or ERBB2 amplification) drive the differentiation of the malignant cells towards a luminal fate.

Global gene expression profiling has fundamentally changed the way we look at cancer cells by providing simultaneous measures of the activity of thousands of genes. We are optimistic that the quantitative ER pathway data generated by microarray experiments will better predict patient prognosis and clinical response to anti-estrogen therapy. Although many significant challenges remain, these technologies will undoubtedly play an important role as we identify new therapeutic targets and improve the efficacy of current breast cancer treatments.