Our study design was approved by the medical ethical committee of the Erasmus MC Rotterdam, the Netherlands (MEC 02.953). To evaluate the predictive value of gene expression profiling in relation to tamoxifen treatment in patients with recurrent breast cancer, 112 fresh frozen ER-α–positive (≥ 10 fmol/mg of protein) primary breast tumor tissue specimens of patients with primary operable (invasive) breast cancer diagnosed between 1981 and 1992 were included. The median age at time of primary surgery (breast-conserving lumpectomy, 33 patients; modified mastectomy, 79 patients) was 60 years (range, 32 to 89 years). In this retrospective study all patients were selected for development of disease recurrence (14 with local or regional relapse, 86 with distant metastasis) that was treated with tamoxifen (40 mg daily) as first-line treatment. At the start of tamoxifen treatment, the median age was 63 years (range, 33 to 90 years), and 27 patients (24%) were premenopausal. None of the patients had received endocrine (neo)adjuvant systemic therapy, nor were they exposed to any hormonal treatment at an earlier stage.

Eighteen patients (16%) received adjuvant chemotherapy. Of these patients, seven were postmenopausal and 11 were premenopausal at time of surgery. At the start of tamoxifen monotherapy, eight patients were still premenopausal, whereas three patients changed to the postmenopausal status before recurrence. Two of these three patients showed objective response to tamoxifen. Therefore, chemical castration as prior endocrine therapy could not have had a significant impact on our results. The median follow-up of living patients was 94 months (range, 21 to 165 months) from primary surgery, and 53 months (range, 2 to 131 months) from the start of tamoxifen treatment. Tumor progression after tamoxifen occurred in 103 (92%) of the patients. During follow-up, 94 patients (84%) died. After tumor progression during first-line tamoxifen treatment, 69 patients were treated with one or more additional endocrine agents, whereas 64 patients were subsequently treated with one or more regimens of chemotherapy (mainly cyclophosphamide, methotrexate, and fluorouracil, or fluorouracil, doxorubicin, and cyclophosphamide) after the occurrence of hormonal resistance.

Criteria for follow-up, type of response, and response to therapy were defined by standard International Union Against Cancer criteria of tumor response,¹² and PFS results were described previously by us.¹³ Complete and partial responses were observed in 12 and 40 patients, respectively, resulting in 52 patients with an OR; PD within 3 to 6 months from start of treatment was observed in 60 patients. Median PFS time of OR was 17 months, whereas the median PFS time of patients with PD was 3 months.

Total RNA was isolated from 30-μm frozen sections (approximately 20 to 50 mg tumor tissue) with RNABee (Campro Scientific, Veenendaal, the Netherlands). The percentage of tumor cells was determined in two frozen 5-μm sections stained with hematoxylin and eosin, which were cut before and after sectioning for RNA isolation. The tumor samples had a median tumor content of 65%. A T7dT oligo primer was used to synthesize double-stranded cDNA from 3 μg total RNA and to generate antisense RNA by in vitro transcription with T7 RNA polymerase (T7 MEGAscript High-Yield Transcription kit; Ambion Ltd, Huntingdon, United Kingdom). Two micrograms of aRNA was labeled with Cy3 or Cy5 (CyDye, Amersham Biosciences, Roosendaal, the Netherlands) in a reverse transcription reaction. The labeled cDNA probes were purified using Qiagen polymerase chain reaction (PCR) cleanup columns (Westburg BV, Leusden, the Netherlands). Similar to the Stanford protocol, a cell line pool of 13 cell lines derived from different tissue origins was used as reference for all microarray hybridizations (details are available at MIAMExpress [www.ebi.ac.uk/miamexpress/]). Probes of the cell line pool were always labeled with Cy5.

Total RNA isolated for the microarray analysis was used to verify the quantity of specific messengers by real-time PCR. The RNA was reverse transcribed and real-time PCR products were generated in 35 cycles from 15 ng cDNA in an ABI Prism 7700 apparatus (Applied Biosystems, Foster City, CA) in a mixture containing SYBR-green (Applied Biosystems) and 330 nmol/L primers for differentially expressed genes (ie, CASP2, DLX2, EZH1, CHD6, MST4, RABEP, SIAH2, and TNC). SYBR-green fluorescent signals were used to generate cycle threshold values from which mRNA ratios were calculated when normalized against the average of three housekeeping genes (ie, hypoxanthine-guanine phospho-ribosyltransferase [HPRT], porphobilinogen deaminase [PBGD], and beta₂-microglobulin [B2M]).¹⁴

Microarray slides were manufactured at the Central Microarray Facility at the Netherlands Cancer Institute (NKI).¹⁵ Sequence-verified clones obtained from Research Genetics (Huntsville, AL) were spotted with a complexity of 19,200 spots per glass slide using the Microgrid II arrayer (Biorobotic, Cambridge, United Kingdom) The geneID list can be found at http://microarrays.nki.nl. Labeled cDNA probes were heated at 95°C for 2 minutes and added to preheated hybridization buffer (Slide hyb buffer 1; Ambion, Huntingdon, UK). The probe mixture was hybridized to cDNA microarrays for 16 hours at 45°C.

Fluorescent images of microarrays were obtained by using the GeneTAC LS II microarray scanner (Genomic Solutions; Perkin-Elmer, Milan, Italy). IMAGENE v5.5 (Biodiscovery, Marina Del Rey, CA) was used to quantify and correct Cy3 and Cy5 intensities for background noise. Spot quality was assessed with the flagging tool of IMAGENE, which in this study was set at R more than 2 for both Cy3 and Cy5. Fluorescent intensities of each microarray were normalized per subgrid using the NKI MicroArray Normalization Tools (http://dexter.nki.nl) to adjust for a variety of biases that affect intensity measurements (eg, color, print tips, local background bias).¹⁶ All ratios were log 2 transformed.

Microarray data analyses were performed with the software packages BRB Array Tools, developed by the Biometric Research Branch of the US National Cancer Institute (http://linus.nci.nih.gov/BRB-ArrayTools.html) and Spotfire. BRB was implemented for statistical analysis of microarray data, whereas Spotfire was used for cluster analysis. The class comparison tool in BRB combines a univariate F test and permutation test (n = 2,000) to find discriminating genes and to confirm their statistical significance. In the class comparison a significance level of .05 was chosen to limit the number of false-negative results.

Spotfire was used to perform hierarchical clustering. To analyze microarray data from different batches of slides, genes were Z score normalized per batch. The Z score was defined as (value − mean)/standard deviation. After normalization, microarray data were clustered via complete linkage. The similarity measure for clustering was based on cosine correlation and average value.

Sensitivity, specificity, positive predictive value, negative predictive value, and odds ratios were calculated and presented with their 95% CIs. The performance of the signature in the validation set was determined via the likelihood ratio of the χ² test. A supervised-learning approach has been applied to reduce our 81 differentially expressed genes to a smaller 44-gene predictive signature. First, all 81 genes were rank ordered on the basis of their significance as calculated with the BRB class comparison tool. Next, starting with the most significant gene, the Pearson correlation coefficient of expression with the other 80 genes was calculated. Succeeding genes were excluded from the signature as long as their expression correlated significantly (P < .05) with the most significant gene. The first gene of the 81-gene profile that did not correlate with expression of the most significant gene was added to the final signature, and the whole procedure of expression correlation analysis with this second gene was repeated with the remaining less significant genes. In this way, genes with overlap in their expression would be removed and the 44-gene predictive signature was derived.

The predictive score for the traditional factors–based model included menopausal status, disease-free interval (DFI > 12 v ≤ 12 months after primary surgery), dominant site of relapse (relapse to viscera or bone v relapse to soft tissue), log ER and log progesterone receptor (PgR) levels. In the analyses of PFS, the Cox proportional hazards model was used to calculate the hazard ratios and their 95% CIs. Survival curves were generated using the method of Kaplan and Meier, and a log-rank test for trend was used to test for differences. Correlation between microarray data and real-time PCR data was determined with Spearman rank correlation test. Computations were performed with the STATA statistical package, release 7.0 (STATA Corp, College Station, TX). All P values are two sided.

For the validation of our 44-gene signature, a classification algorithm (gene prediction tool) was developed that is comparable to the compound covariate predictor from BRB Array Tools. The gene prediction tool applies two cutoff values instead of the midpoint used in the compound covariate predictor tool for classification. The two thresholds are the median values of PD and OR, and are defined in the tumors of our training set. To obtain a robust classification algorithm, genes from the signature become only classifiers whenever the expression values are outside the two thresholds and as a result mainly represent one class, either PD or OR. When the expression level falls between the cutoff values, the gene is excluded as a classifier because the value can represent both response classes (ie, PD and OR). The gene classifiers from the predictive 44-gene signature are identified for each tumor from the validation set using the above-described algorithm. Finally, the ratio between the identified response-predicting genes and resistance-predicting genes determines the predicted signature-based response outcome.

To select discriminatory genes for the type of response to tamoxifen, a training set of 46 tumors was defined that comprised primary tumors of 25 patients with PD and tumors of 21 patients with OR (Fig 1). The tumor RNAs of this training set were hybridized in duplicate, and genes or expression sequence tags that had more than 10% absent calls over the experiments were eliminated. This resulted in 8,555 and 7,087 assessable spots, respectively. Using a significance level of .05 in the BRB class comparison tool, 569 and 449 genes, respectively, were differentially expressed between the PD and OR subsets. The overlap (ie, 81 genes) was designated as the differentially expressed signature. After supervised hierarchical clustering (Fig 2), this discriminatory signature correctly classified 21 of 25 patients with PD (84% sensitivity; 95% CI, 0.63 to 0.95) and 19 of 21 patients with OR (91% specificity; 95% CI, 0.68 to 0.98), with an odds ratio of 49.8 (P < .0001). The positive predictive value and negative predictive value for resistance to tamoxifen was 91% and 83%, respectively.

Additional analysis (rank ordering of genes on the basis of significance level followed by a step-up calculation of correlation coefficient of expression) reduced the initial set of 81 genes to a smaller 44-gene predictive signature with similar accuracy.

The mRNA expression levels of eight signature genes were also measured by quantitative real-time PCR (ie, CASP2, DLX2, USP9X, CHD6, MST4, RABEP, SIAH2, and TNC). The quantitative PCR data were correlated with the microarray data. Spearman rank correlations were positive for seven of eight genes.

The 81-gene signature consisted of 15 ESTs and 66 known genes (Fig 2). Functional annotation of the genes in our signature showed that they were involved in estrogen action (26%), apoptosis (14%), extracellular matrix formation (9%), and immune response (6%). The remaining genes had a function in glycolysis, transcription regulation, and protease inhibition.

Seventeen genes were regulated by or associated with estrogen (receptor) action, of which nine genes showed upregulation and eight genes showed downregulation in tamoxifen-resistant tumors (Fig 2). Clustering of the signature revealed a cluster of six genes that was associated with the extracellular matrix (ECM; TIMP3, FN1, LOX, COL1A1, SPARC, and TNC), which were overexpressed in patients with resistant disease. In an exploratory analysis to search for associations between genes, using Ingenuity Pathway Analysis tools, seven genes of the signature (IL4R, LDHA, MAP2K4, NPM1, SIAH2, CASP2, and TXN2) were associated with apoptosis, whereas two other signature genes (API5, BNIP3) were related to antiapoptosis processes. Interestingly, four apoptosis genes (API5, NPM1, LDHA, and BNIP3) were upregulated and five genes (IL4R, MAP2K4, SIAH2, CASP2, and TXN2) were downregulated in tumors with resistant disease. A cluster of four genes linked to the immune system (FCGRT, PSME1, HLA-C, and NFATC3) was downregulated in PD.

Notably, the 81-gene signature showed an over-representation of genes located on chromosome 17, but an under-representation of genes located on chromosomes 4, 15, 18, and 21 (Fig 4). In particular, genes localized to cytoband 17q21 to 17q22 seemed to be significantly (P = .03) over-represented (ie, five of 66 informative genes [APPBP2, COL1A1, EZH1, KIAA0563, and FMNL]) in the signature (6.5%) compared with 199 of 12,771 known genes (1.5%) for the whole microarray.