Search for a Gene Expression Signature of Breast Cancer Local Recurrence in Young Women

Breast-conserving treatment (BCT) is established as the preferred treatment for patients with early-stage breast cancer (1). It offers both local control and improved survival (2), as well as superior psychosocial outcomes compared with patients treated with modified radical mastectomy (3, 4). However, the occurrence of an ipsilateral breast cancer can be traumatic and responsible for cancer-specific death (2). Several risk factors for developing an ipsilateral breast cancer after BCT with combined surgery and radiotherapy have been identified. The fact that (young) age is the single most important risk factor, independent of all others (5), suggests that some unrevealed biologic characteristics of the primary breast cancer may be associated with a higher risk of developing a true local recurrence, that is, regrowth of clonogenic cells that were not removed by surgery or killed by radiotherapy (6). Gene expression profiling of primary breast carcinomas using microarray technology has already been used to try and discover a signature associated with a higher risk of developing true recurrence after BCT (7–10). Results have been contradictory. The most recent and largest study, by Kreike and colleagues (10), is for a series of 165 young (≤50 years old) premenopausal Dutch patients using the Human Genome Oligo Set version 3.0 array (Operon) and did not confirm previous classifiers, namely, those based on the wound response signature (8) and that defined by Niméus-Malmström and colleagues (9). Kreike and colleagues were able to construct a local recurrence classifier based on the expression of 111 genes and validate this profile on an independent data set of 161 consecutive patients with breast cancer treated by BCT who used a different microarray platform (10). This classifier was mostly characterized by proliferation but did not bring a significant independent additive value, as young age remained the sole predictive factor of local recurrence in multivariate analysis.

This study reports on a series of 343 patients with primary breast cancers selected from the fresh-frozen tissue bank of 4 Dutch hospitals (Netherlands Cancer Institute, Amsterdam; Radboud University Nijmegen Medical Centre, Nijmegen; Erasmus Medical Center, Rotterdam; and Ziekenhuis Amstelland, Amstelveen, the Netherlands) and one French cancer center (Institut Curie, Paris, France). The 148 Dutch patients belonged to a larger series of 165 patients with breast cancer already reported (10) were combined with the 195 French patients. Patients were included in the study if they met the following criteria: invasive primary breast carcinoma diagnosed at or before the age of 50; all patients were premenopausal; and no previous history of cancer (except for one nonmelanoma skin cancer). All patients were treated between January 1984 and November 2002 by initial breast-conserving surgery (216 pT1, 124 pT2, 1pT3, 2 missing data), including either axillary lymph node dissection or sentinel node procedures (215 pN0, 99 pN1, 21 pN2, 3pN3, 6 missing data). All patients underwent postoperative radiotherapy to the breast (median dose of 50 Gy; range, 45–55) with, in 248 patients, a boost to the tumor bed (external beam radiotherapy or brachytherapy with a median dose of 15 Gy; range, 6–26) and/or the regional lymph node–bearing areas. One hundred and forty-one patients received adjuvant systemic treatment with chemotherapy only (110 patients), hormone therapy only (8 patients), or a combination of both (23 patients). Chemotherapy regimens (4 missing data) were either anthracycline-based (79 patients) or contained cyclophosphamide, methotrexate, and fluorouracil (CMF, 50 patients). Patients were selected if they had either experienced a local recurrence (ipsilateral breast recurrence with clinical, histologic, and/or genomic features consistent with true recurrences) within the first 10 years after primary treatment or remained free from local recurrence (any type of ipsilateral breast cancer) for at least 10 years after primary treatment. To ensure relevant data, we conducted gene expression analyses in tumors with at least 50% of cancer cells as assessed by hematoxylin and eosin staining of histologic sections of the snap-frozen samples.

This research was approved by the Institutional Review Boards of both the Netherlands Cancer Institute and the Institut Curie.

Clinical and tumor characteristics of the 343 patients are reported in Table 1. Three pathologists (M.J. van de Vijver, B. Sigal-Zafrani, and R. Hadhri) reviewed histopathologic characteristics of all tumors. Immunohistochemical (IHC) analysis was used to determine the status of estrogen receptor (ER; clone 6F11, 1:200 dilution; Novocastra) and progesterone receptor (PR; clone 1A6, 1:200 dilution; Novocastra). Tumors were considered to be positive for these receptors if at least 10% of the invasive tumor cells in a section showed nuclear staining (11). HER2 overexpression was determined according to American Society of Clinical Oncology (ASCO) guidelines.

Tumor material was snap-frozen in liquid nitrogen within 1 hour of surgery. Thirty 30-μm sections were used for isolation of RNA. The first and last sections (5 μm) were stained with hematoxylin and eosin; only samples containing an average of at least 50% tumor cells were used in this analysis. Two protocols of RNA isolation were used. For Dutch samples, total RNA was isolated with RNA-Bee (Tel-Test) and dissolved in RNase-free water. Then, total RNA was treated with DNase with use of the Qiagen RNase-Free DNase Kit and RNeasy spin columns and dissolved in RNase-free water. Detailed information on protocols can be found at the Central Microarray Facility website of the Netherlands Cancer Institute (http://microarrays.nki.nl/). For French samples, RNA was extracted from frozen samples using the TRIzol method (Invitrogen) according to the manufacturer's instructions and purified using miRNeasy kit (Qiagen). The concentration, integrity, and purity of each RNA sample (260/280, 260/230, 28S/18S, RIN) were measured using the RNA 6000 LabChip Kit with the Agilent 2100 Bioanalyzer (Agilent Technologies). The gene expression microarrays used in this study were Illumina Human Whole Genome V3.0 arrays (Illumina, Inc.) containing 48,803 reporters. Details of RNA amplification, labeling, and hybridization are available on the Illumina website (http://www.illumina.com). The arrays were processed in the Central Microarray Facility of the Netherlands Cancer Institute and are publicly available at Gene Expression Omnibus (GEO), series record GSE30682 (http://www.ncbi.nlm.nih.gov/geo/).

The population of 148 Dutch patients was first studied to validate the signature developed by Kreike and colleagues (10) and to assess cross-platform reproducibility. For these 148 patients, both Operon (Human Genome Oligo Set version 3.0) and Illumina (Human Whole Genome V3.0) data were available. A multiple factor analysis (MFA; ref. 12) was applied to the 111-gene signature reported by Kreike and colleagues using the R-package FactoMineR (v1.14; ref. 13). MFA aims to integrate the 2 groups of variables (111 genes expression profiles from both Operon and Illumina data), which are expected to describe the same observations (148 tumor samples). The goal is to check whether both data sets share a common geometric structure, thus indicating that the same underlying biologic signal can be extracted from the same samples and the same set of genes, with different technology. The local recurrence status was added as an illustrative variable.

The previous local recurrence signature developed at the NKI was tested in the French population along with 21 other signatures: mutation status (14), chromosomal instability (15), fibroblast serum response (16), wound serum response (17), hypoxia signature (18), radioresistance signature (19), invasiveness signature (20), local recurrence signature (9), 70-gene NKI signature (21), p53 signature (22), recurrence score (23), Rotterdam signature (24), PTEN loss signature (25), genomic grade (26), proliferation signature (27), molecular portraits (28), prognostic gene expression classifier for ER-positive breast cancer (29), pooled signature (30), markers of proliferation (31), 3 modules (32), and the wound signature local relapse signature (8).

A diagonal linear discriminant analysis (DLDA) was applied to classify the samples and the performances were evaluated using Monte Carlo cross-validation (10-fold—100 permutations; ref. 33). The area under the receiver operating characteristic (ROC) curve (AUC) was used as the prediction quality criterion.

The 343 arrays were normalized using the variance stabilization transformation and the robust spline algorithm suggested by Du and colleagues (lumi R-package v2.2.0; ref. 34). Probes were then filtered to discard those called as unexpressed (or undetectable) in all samples. This left 30,198 probes to be used for analyses. To ensure comparability between the Dutch and French samples, the data set was corrected for a population effect using the following one-way ANOVA model applied to each probe i, i ∈ [1, 30198]:

where Y_ij is the expression of probe i in sample j ∈ [1, 343], O_p is the effect of the population, p ∈ {Dutch, French}, μ_i is the expression of gene i when p = Dutch, ϵ_ijp ∼ N(0, σ_i²) is a residual term. The corrected expression level of probe i was computed as follows:

Ward linkage hierarchical clustering of Pearson correlation similarity matrices was then applied on the 5,000 most variant probes (highest interquartile range) of the 343 samples, using the R-package EMA (v1.2; ref. 35).

On the basis of the expression level of genes ESR1, PR, ERBB2, and AURKA, each sample was assigned to a molecular subtype. All tumors were assessed using the average expression ratio of the probes for each gene, respectively, ESR1 (1731003), PR (1710082), ERBB2 (1729826/1753944), and the AURKA (1710733/1755528) proliferation marker. ROC curves were computed for all tumors with known IHC data. The microarray-derived ER, PR, and HER2 statuses were determined from these tumors according to the optimal sensitivity–specificity cutoff point (Fig. 1; ref. 36) and then applied to all tumors without IHC data. Tumors were grouped into molecular subtypes according to their gene expression, as defined by Voduc and colleagues (37): luminal A (ER⁺ and/or PR⁺, HER2⁻, low AURKA), luminal B (ER⁺ and/or PR⁺, HER2⁻, high AURKA), luminHER (ER⁺ and/or PR⁺, HER2⁺), triple-negative (ER⁻, PR⁻, HER2⁻), and HER2 enriched (ER⁻, PR⁻, HER2⁺).

In addition, the single sample predictor (SSP) provided by Hu and colleagues (28) was mapped onto our microarray platform and tested for class assignment. To do this, the expression data were centered on genes and the Pearson correlation of each sample with the centroids of the 5 molecular subtypes (luminal A, luminal B, HER2⁺, basal-like, normal-like) was calculated. Each sample was assigned to a subtype on the basis of the largest correlation. Samples with largest correlation less than 0.1 were considered unclassified.

A one-way ANOVA model, taking into account the subtype information and its interaction with relapse, was applied to all probes i, i ∈ [1, 30198] to identify differentially expressed genes (DEG) between local recurrence and no local recurrence:

where is the expression of gene i in sample j ∈ [1, 343] corrected from the population effect (A), T_t is the effect of the molecular subtype, t ∈ {HER2⁺, luminal A, luminal B, luminal HER, triple-negative}, LR_r is the effect of local recurrence, r ∈ {0, 1}, (T × LR)_tr is the interaction term between molecular subtype and recurrence event, μ_i is the expression of gene i when t = HER2⁺ and r = 0, ϵ_ijtr ∼ N(0, σ_i²) is a residual term.

The analysis was conducted using the R-package limma (v3.6.9) and the P values were adjusted for multiple testing using the Benjamini–Hochberg correction. All adjusted P values smaller than 5% were considered significant.

Overrepresentation of Gene Ontology categories in the list of DEGs was identified using a hypergeometric test (R-package GOstats v2.16.0). Categories represented by less than 10 genes were discarded.

The distinction of secondary primary cancers and true recurrences is an important issue, which can be used to validate patient selection. In the previous study, Kreike and colleagues have already assessed this question on the Dutch samples, using the overall gene expression pattern of 15 pairs of primary tumors and local recurrences. Regarding the French series, we analyzed 14 pairs of tumors using DNA copy number alteration (Affymetrix GeneChips Human Mapping 50K array Xba) profiles and the partial identity score proposed by Bollet and colleagues (38). Briefly, this score aims to estimate the clonal relatedness of the ipsilateral breast cancer and the primary tumor using the number of common DNA breakpoints among both profiles, weighted by their occurrence frequency in a control population.

After appropriate scaling of the set of clinical variables (Table 1) and expression data, we used Monte Carlo cross-validation (5-folds—100 permutations) and conducted binary classification using the support vector machines (SVM) algorithm (39) to predict outcome (local recurrence or no local recurrence). A radial basis function (rbf) kernel was used, with variance fixed at 10 for the clinical data and 1,000 for expression data; internal cross-validation was carried out inside each fold to select the SVM penalty parameter λ. We applied the algorithm separately for clinical and expression data, as well as combining the 2 using either the mean or the product of the individual rbf kernels to create a new valid kernel (39). In experiments, we calculated the balanced accuracy, that is, 0.5 (sensitivity + specificity). Analyses were first conducted on all patients, then for patients with each clinical subtype to see if this led to improvement in accuracy for certain subtype.

Furthermore, we applied the least angle regression (LARS) algorithm (40), which does both supervised classification and variable selection, and can help to suggest a ranking of important predictive variables.

We investigated the relationship between local recurrence and gene expression profiles of primary breast cancers of patients treated with BCT in a series of 343 patients. Of these, 119 patients had local recurrence (median time lapse of 3.3 years; range, 0.3–9.4 years) and 224 patients remained free from local recurrence for at least 10 years.

The clinical and pathologic characteristics of patients and tumors are displayed in Table 1. Patients who developed local recurrences were significantly younger at diagnosis than those who did not (P = 0.044; χ² test). Their tumors had significantly higher histologic grade according to Ellis and Elston (P = 0.012; χ² test), higher mitotic indexes (P = 0.011; ordinal χ² test), and more frequently contained an extensive ductal carcinoma in situ (DCIS) component (P = 0.041; χ² test).

The ER, PR, HER2, and proliferation statuses were determined from the gene expression data according to the optimal sensitivity–specificity cutoff from the ROC curves produced using IHC data (Fig. 1). These statuses were respectively defined from the probes ESR1 (1731003, cutoff = 8.687), PR (1710082, cutoff = 7.405), ERBB2 (1729826/1753944, cutoff = 9.481), and AURKA (1710733/1755528, cutoff = 8.622).

Of the complete data set of 343 patients, 148 were also used by Kreike and colleagues to define the 111-gene signature and studied using the Operon microarray platform.

MFA (Fig. 2) aims to integrate the 2 groups of variables (111-gene expression profiles from both Operon and Illumina data) which are expected to describe the same observations (148 tumor samples) to assess the robustness of this signature on another microarray platform, as well as its association with local recurrence. The 2 groups have coordinates that are close in the first dimension, which means that their contribution to the first principal component is similar. It also means that the first principal component of the MFA (the strongest signal) is common to the 2 data sets. The local recurrence event is well represented on the first component of the MFA. This observation is consistent with the performance of the 111-gene signature previously described.

Another way to assess the predictive power of the signature from Kreike and colleagues is to consider the French samples as an independent validation series. The class prediction (DLDA) on this population showed an AUC of 0.56, a specificity of 0.53 for a sensitivity of 0.63, and a misclassification rate of 43%. These results are worse than the initial performances published by Kreike and colleagues (sensitivity, 0.77 and specificity, 0.43). We assessed how other gene expression signatures conducted in pairing tumors to local recurrence status. The 111-gene classifier did not show convincingly better performance on the French samples. None of the tested signatures were able to classify the samples with an average balanced accuracy higher than 0.585 (Supplementary Table S1).

Hierarchical clustering of the 343 tumors and the 5,000 most variant probes resulted in the formation of 2 main clusters that could be further subdivided into 4 subclusters (Fig. 3). None of the tumor characteristic related to local recurrence was associated with these clusters.

These 4 clusters were closely related to the subtype definitions. We clearly identified the triple-negative tumors and 2 other clusters mainly composed of luminal A and luminal B–like tumors. The luminal A, luminal B, and the triple-negative clusters are clearly identified with, respectively, 72%, 76%, and 96% (Supplementary Table S2; hypergeometric tests, P = 2.86e-36, P = 1.60e-18, and P = 1.46e-53) of the tumors in the 3 groups. The last cluster is enriched in HER2 profiles and can be further divided in the HER2⁺ profiles and a mixture between luminal HER2 and luminal B profiles. The class assignments predicted from the study of Hu and colleagues (28), SSP, do not outperform the clustering based on the ESR1, PR, HER2, and proliferation subtype prediction (Supplementary Fig. S1).

Information on breast cancer subtypes is crucial and has to be taken into account in any supervised analysis to avoid confounding effects. The association between breast cancer subtype and relapse events (Supplementary Table S3; P = 0.033; χ² test) has already been reported in several previous studies (10, 37). To identify DEGs between patients developing a local recurrence and those who do not, we used the linear model defined in the work of Clarke and colleagues (2). P values lower than 5% after Benjamini–Hochberg correction were reported as significant.

No DEGs in the groups luminal A, luminal B, and luminal HER were found. A list of 182 probes was generated for the triple-negative tumors and a list of 115 probes for the HER2⁺ tumors. When the subtype factor was removed from the linear model, no genes were found to be significant.

To identify whether these 2 lists of probes were enriched for genes representing particular Gene Ontology categories compared with the entire list of probes, a hypergeometric test was used and the Gene Ontology categories with P values lower than 0.001 were reported. Overrepresentation of genes related to proliferation, immune response, and cell death were found in the triple-negative differentially expressed probes (DEG). No interesting category was found in the list of DEGs in the HER2 overexpressed population (Supplementary Table S4).

We applied the partial identity score proposed by Bollet and colleagues (38) on 14 pairs of tumors. Among the 14 pairs, all but one had a significant score and thus were classified as a true recurrence (Fig. 4). This indicates that the DNA copy number profile of the local recurrences and the primary tumors are highly similar. It also validates the inclusion of the French samples and the relevance of their recurrence status.

Using the SVM algorithm with an rbf kernel, we obtained an average balanced accuracy of 62% ± 2% using clinical variables and 59% ± 2% using expression data. Using a mean kernel to combine the 2 sets of variables gave 64% ± 2% average balanced accuracy whereas a product kernel gave 63% ± 2%. This indicates that combining the 2 types of variables in this way may give a small prediction accuracy improvement compared with treating them separately.

The same analyses with the SVM algorithm and rbf kernel were conducted by subtype (Table 2). There were too few HER2 patients to give meaningful results. Balanced accuracy results did not significantly improve with respect to treating all subtypes together, and in some cases, decreased.

Using LARS on the clinical data with 5-fold cross-validation and averaged over 1,000 trials gave a maximum best balanced accuracy of 63% obtained when including 32 variables in the predictor. Running LARS on all clinical data for variable selection purposes, we obtained a list of pertinent variables along with their weight in the prediction function at the maximum where 32 variables had been included. The first 10 selected variables are shown in Table 3.

Including the 5 subtypes as additional binary variables did not change the above ranking. Using LARS on the expression data with 5-fold cross-validation and averaged over 100 trials gave a maximum balanced accuracy of 58% obtained when including 99 variables in the predictor. Using LARS on the combined clinical and expression data with 5-fold cross-validation and averaged over 100 trials gave a maximum balanced accuracy of 61% obtained when including 73 variables in the predictor.

We have conducted gene expression profiling of 343 primary breast cancers that were treated with BCT. We chose to base our study on a series of young (≤50 years old), premenopausal women, not only because (young) age is recognized as one of the most important independent prognostic factors for ipsilateral breast recurrence (5, 41–46) but also to ensure a high level of homogeneity. In addition, all patients had undergone breast-conserving surgery followed by whole-breast radiotherapy for their initial breast cancers.

In a previous study, Kreike and colleagues (10) found that there were no significant differences in the overall gene expression profiles between primary tumors with or without local recurrence. With an increased number of patients and the use of a new microarrays platform (Operon), they were able to construct and to validate a predictive classifier for local recurrence with a specificity of 43% and sensibility of 77%. We have significantly extended the number of patients to include 343 primary breast cancers, using a different microarray platform (Illumina). From these 343 samples, 148 were common with the previous study of Kreike and colleagues. Among the 165 new French samples, the single-nucleotide polymorphism profiles of 14 pairs of primary tumors and ipsilateral breast cancers were analyzed. The relapse status of all but one patient was validated.

The first part of this study consisted in confirming the previous 111-gene signature. As expected, the MFA analysis on the 148 Dutch samples common to both studies and platforms (Operon and Illumina) on the 111 genes resulted in good correlation between the 2 data sets, and further, shows an association between local recurrence and the expression profiles of these 111 genes. However, when we applied this signature to an independent data set (French samples), we were not able to reproduce the reported performance, finding a sensitivity of 63% and a specificity of 53%. We also assessed a large number of alternative gene expression signatures and found that the 111-gene classifier did not show convincingly better performance than the others. With a large number of false positives (specificity = 53%) and many false negatives (sensitivity = 63%), this signature has limited value in the prediction of local recurrence.

To address this shortcoming, we investigated the gene expression profile of a significantly larger set of 343 tumors. Hierarchical clustering analysis of the 343 tumors showed clusters of specific breast cancer subtypes. The triple-negative, luminal A, and luminal B samples were well separated. The HER2⁺ samples were also well separated from the luminal HER2 which are grouped with some luminal B samples. The subtype prediction is a potential confounding factor that may have strong impact on the results of the study. Recently, Weigelt and colleagues (47) discussed the robustness of several microarray-based SSP and the need to establish a standard methodology for breast cancer subtype definition. Here, the combination of IHC status and gene expression data offers effective subtype discrimination. Moreover, the IHC status gives robust information that is mainly used in routine clinical practice.

We did not observe significant grouping of the tumors with respect to their relapse status. Previous studies (10, 37) have already described the association between molecular breast cancer subtypes and local recurrence. This observation is critical for further analysis, and particularly for the search of DEGs. On the basis of the linear model (2) which takes into account the molecular subtype and its association with local recurrence, we found significant DEGs only for the HER2⁺ and triple-negative subtypes. The DEGs found in the triple-negative cancers were associated with Gene Ontology categories such as proliferation, cell death, and immune response. The enrichment in genes involved in proliferation was also reported by Kreike and colleagues (10) and agrees with previously published classifiers. The fact that no interesting categories were found in the HER2⁺ overexpressed subtype might be explained by the small number of patients in this class (8 local recurrence and 6 no local recurrence). Further investigations and a larger number of patients with HER2⁺ tumors are necessary to optimize the search for DEGs in this tumor subtype.

More surprisingly, no DEGs were found if we did not take into account subtype information. This observation is remarkable and contradicts the previous results published by Kreike and colleagues (10). The 2 different microarrays platforms, the heterogeneity of the data and the different subtype proportions might have some influence on this finding.

On the basis of the 343 tumors profiles, we attempted to define a new classifier using clinical data and microarray gene expression profiling. Combining both variable types gave a small balanced accuracy improvement compared with treating them separately.

We conducted supervised classification (SVM and LARS) on both the clinical and expression data, to create classifiers to predict local recurrence or no local recurrence. Using SVM, balanced accuracy between 60% and 64% was achieved and remains unsatisfactorily low. Separation of the data into subtypes gave similar balanced accuracy, although in some case worse. This may be partially because of a loss of predictive power in smaller data sets when using supervised classification.

Using LARS on the clinical data confirmed that the total dose of radiotherapy to the tumor bed and age at diagnosis were the 2 most important prognostic factors associated with decreased risk of local relapse (5). The third and fourth factors were the surgical margin involvement by invasive component and the tumor proliferation, as already described by Kreike and colleagues (10, 48). Adding subtype information to the clinical variables did not improve significantly the algorithm's performances. Using LARS on expression data or on combined clinical and expression data did not give improved classification performance.

In conclusion, although we were unable to validate the 111-gene classifier proposed by Kreike and colleagues (10) and to define a strong classifier for local relapse after BCT from this large and well-characterized data set, it cannot be concluded that such a classifier does not exist. New studies, such as one that is currently under progress at the Netherlands Cancer Institute, the Institut Gustave Roussy (Villejuif, France), and the Karolinska Institutet (Stockholm, Sweden) on preoperative radiotherapy (image-guided preoperative accelerated partial breast irradiation (PAPBI): defining radiotherapy sensitivity RCB 2010-A00573-36) should help in finding new factors associated with radiosensitivity in breast cancer. As previously reported, our study also showed that proliferation, age, and tumor invasion are important criteria in predicting local recurrence. Finally, this study also highlights the difficulty of reproducing and validating gene expression signatures between competing microarray platforms and from one population to another.

No potential conflicts of interests were disclosed.

The authors thank F. Reyal for helpful discussion about this study and the manuscript.

This work was financed by the Institut Curie, the ‘Courir pour la vie, Courir pour Curie’ association, the ‘Odyssea’ association, and PHRC 2006.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.