Abstract
Purpose: Several prognostic gene expression profiles have been identified in breast cancer. In spite of this progress in prognostic classification, the underlying mechanisms that drive these gene expression patterns remain unknown. Specific genomic alterations, such as copy number alterations, are an important factor in tumor development and progression and are also associated with changes in gene expression.
Experimental Design: We carried out array comparative genomic hybridization in 68 human breast carcinomas for which gene expression and clinical data were available. We used a two-class supervised algorithm, Supervised Identification of Regions of Aberration in aCGH data sets, for the identification of regions of chromosomal alterations that are associated with specific sample labeling. Using gene expression data from the same tumors, we identified genes in the altered regions for which the expression level is significantly correlated with the copy number and validated our results in public available data sets.
Results: Specific chromosomal aberrations are related to clinicopathologic characteristics and prognostic gene expression signatures. The previously identified poor prognosis, 70-gene expression signature is associated with the gain of 3q26.33-27.1, 8q22.1-24.21, and 17q24.3-25.1; the 70-gene good prognosis profile is associated with the loss at 16q12.1-13 and 16q22.1-24.1; basal-like tumors are associated with the gain of 6p12.3-23, 8q24.21-22, and 10p12.33-14 and losses at 4p15.31, 5q12.3-13.1, 5q33.1, 10q23.33, 12q13.13-3, 15q15.1, and 15q21.1; HER2+ breast show amplification at 17q11.1-12 and 17q21.31-23.2 (including HER2 gene).
Conclusions: There is a strong correlation between the different gene expression signatures and underlying genomic changes. These findings help to establish a link between genomic changes and gene expression signatures, enabling a better understanding of the tumor biology that causes poor prognosis. Clin Cancer Res; 16(2); 651–63
This study reports the specific chromosomal and regional aberrations related to clinical and pathologic characteristics and prognostic gene expression signatures. Using gene expression data (mRNA) from the same tumors, we identified genes in the altered regions for which the expression level is significantly correlated with the copy number. We found that gains and losses varied between prognostic gene expression signatures and clinical and pathologic features. For example, the previously identified poor prognosis, 70-gene signature is associated with gain of 3q26-27, 8q22-24, and 17q24-25; the 70-gene good prognosis profile is associated with loss at 16q11-24. These findings help to establish a link between genomic changes and gene expression signatures that determine tumor behavior in breast cancer, enabling a better understanding of the tumor biology that causes poor prognosis. These correlations will also help in the development of improved prognostic tests that may be used in patient-tailored therapy for breast cancer.
In recent years, high-throughput technologies such as gene expression profiling using microarray analysis have improved the possibilities to assess the risk of developing distant metastases in individual breast cancer patients. Gene expression profiling studies have identified various prognostic gene expression signatures in breast cancer: a 70-gene prognosis signature (1), a 76-gene prognosis signature (2), the “molecular subtypes” (3–5), the “wound response signature” (6), the chromosomal instability signature (7), and the genomic grade index (GGI; ref. 8). These signatures are independently derived prognostic factors and show a partial agreement in the identification of patients at increased risk of developing distant metastases (9).
Differences in gene expression between breast tumor samples can be considered to be the result of genetic (10, 11) and/or epigenetic changes (12). Association between copy number change and gene expression levels has been shown, and ∼12% of gene expression variation can be explained by differences in DNA copy number (10). Regulators of molecular processes may be targeted by these alterations, and the expression changes are mostly the result of downstream targets of these regulators reacting to the changes. By combining gene expression and copy number data, these interactions can be revealed and therefore better pinpoint the involved processes.
Several studies have correlated DNA copy number changes with the mRNA expression variation of individual genes (13–16) and, more recently, gene expression signatures (17–22). Adler et al. (23) showed that an activated wound response signature in breast carcinomas is strongly associated with the amplification of MYC (8q24) and CSN5 (8q13), and in vitro experiments have shown that these genes can induce the activated wound signature. Other prognostic signatures, including the basal-like (4) and 70-gene (1) signatures, were not induced by MYC and CSN5 (23). Several additional specific associations have also been reported: basal-like tumors were found to be associated with the amplification of 6p21-25 (17) and 10p14 (24); luminal-type tumors have been found to be associated with several recurrent DNA copy number changes.
Integration of DNA copy number alterations and gene expression profiling may also result in improved classification and prediction of prognosis in breast cancer. For example, Chin et al. (18) who found that the accuracy of risk stratification according to the outcome of breast cancer disease can be improved by combining analyses of gene expression and DNA copy number. They measured both DNA copy number and gene expression profiles in 101 primary breast tumors and showed that high-level amplification and/or overexpression of genes at 8p11, 11q13, 17q12, and or 20q13 were strong predictors of reduced survival duration. In another study, by Chin et al. (20) also measured DNA copy number alterations in 171 breast tumors using high-resolution genome-wide profiling. Together with matched array expression data, they identified genomic regions showing strong coordinate expression changes and frequently amplified regions on 8q22.3 (EDD1, WDSOF1), 8q24.11-13 (THRAP6, DCC1, SQLE, SPG8), and 11q14.1 (NDUFC2, ALG8, USP35) associated with significantly worse prognosis (20).
More recently, Andre et al. (22) used high-resolution oligonucleotide comparative genomic hybridization (CGH) arrays and matching gene expression array data to identify dysregulated genes. DNA was extracted from 106 pretreatment fine-needle aspirations of stage II to III breast cancers that received preoperative chemotherapy. A high number of molecular targets (HER-2, epidermal growth factor receptor, PTEN, VEGFA, AURKA, CHEK, AKT, and FGFR1) showed a significant correlation between DNA copy number alterations and mRNA levels. The region between 8q24.14 and 8q24.22 was gained in 58% of the tumors and the genomic region at 8p11 and 8p12 showed a high level of amplification in 10% of the samples. The most striking finding was that triple-negative tumors showed gain at 6p21, and Andre et al. (22) identified at least three (VEGFA, PTEN, and EGFR) genes that showed a high frequency of DNA copy number alterations that was correlated with gene expression. Even the VEGFA gene was gained in 34% of triple-negative tumors.
We have previously developed an algorithm, Supervised Identification of Regions of Aberration in aCGH data sets (SIRAC; ref. 25), to perform supervised identification of genomic regions enriched for probes correlated with sample labeling (25). Here, we analyzed DNA copy number changes within labeled data sets using SIRAC and integrated these findings with gene expression (mRNA) profiles in a series of 68 primary breast carcinomas. These tumors were labeled based on well-known clinical and pathologic features and gene expression profiles. These included the 70-gene prognosis signature (1), the wound signature (6), GGI (8), and the breast cancer molecular subtypes (3–5). The associations of these subclasses with specific DNA copy number changes and gene expression were investigated.
Materials and Methods
Tumor samples
We selected 68 breast carcinomas, from the fresh-frozen tissue bank of the Netherlands Cancer Institute/Antoni van Leeuwenhoek Hospital, for DNA isolation and DNA copy number analysis. These 68 tumors were part of a series of breast tumor specimens of 295 consecutive women with breast cancer. For each of the 295 samples, gene expression (26) and clinicopathologic data (27) were previously available. Several studies have correlated DNA copy number changes with one-gene expression signature, especially the molecular subtypes of breast cancer (17–21). We wanted to be more comprehensive and these 68 tumors were selected to represent the following 4-gene expression profiles: (a) 70-gene prognosis signature: poor prognosis profile, good prognosis profile (1); (b) wound signature: activated, quiescent (6); (c) molecular subtypes: basal like, HER2 like, Luminal A, Luminal B, and normal epithelial like (3, 4); and (d) GGI: GGI 1 and GGI 3 (8). We selected these 68 tumors based on the following criteria. Our first aim was to select, from the 295, the 15 tumors per molecular subtype that are most similar to the centroid representing that subtype. Our second aim was to have equally distributed subgroups of patients with good and poor prognosis signatures (70-gene prognosis signature, wound signature, and GGI). The Institutional Review Board of the Netherlands Cancer Institute approved this study. All 68 patients had no prior malignancies and did not receive any systemic therapy before surgery (33 had lymph node–negative disease and 35 had lymph node–positive disease). One patient with lymph node–negative disease received adjuvant chemotherapy. Thirty-one of the 35 who had lymph node–positive disease had received adjuvant systemic therapy, of which 25 patients received adjuvant chemotherapy (a regimen of cyclophosphamide, methotrexate, and fluorouracil), 7 received adjuvant hormonal therapy (tamoxifen), and 2 patients got both modalities.
Immunohistohemistry
For all 68 breast carcinomas, estrogen receptor (ER), progesterone receptor (PR), and HER2 status were assessed by analyzing the expression of these proteins using immunohistochemistry on tissue microarrays (TMA) produced by a manual tissue arrayer (Beecher Instruments) according to standard procedures with primary antibodies against ER-α (1D5+6F11; dilution, 1:50; NeoMarkers, Lab Vision Corp.), PR (PgR; R-1; dilution, 1:500; Klinipath), HER2 (3B5; dilution, 1:3,000; ref. 1). Staining for ER and PR was scored as negative when <10% of tumor cells showed nuclear staining and was scored positive when >10% of tumor cells showed staining. HER2 status was evaluated according to commonly accepted pathologic guidelines (28). A tumor was considered to be HER2 positive when a score of 3+ was found; it was considered HER2 negative when a score of 0 was observed; tumors with a score of 1+ or 2+ were evaluated by chromogenic in situ hybridization (CISH) for HER2 assessment; tumors with HER2 gene amplification (six or more spots per tumor cells) were scored as HER2 positive and all other tumors as HER2-negative (29).
Genomic DNA isolation and labeling
DNA was isolated from frozen tumor material by cutting 25 sections of 30-μm-thick sections. A verification of a tumor cell percentage above 50% was done on H&E slides at both ends of the cutting. Frozen tissue sections were digested in 15 mL TNE [100 mmol/L sodium chloride, 10 mmol/L Tris-HCl (pH 8.0), 25 mmol/L EDTA (pH 8.0), and 1% sodium dodecyl sulfate with the addition of proteinase K (50 μg/mL)] and were incubated at 55°C for at least 24 h. Phenol-chloroform extraction was then done followed by ethanol precipitation. DNA was diluted in 10 mmol/L Tris-HCl (pH 7.6)/0.1 mmol/L EDTA. For all cases, array CGH (aCGH) was done using 2 μg of genomic DNA. All labeling reactions were done with the Cy3 and Cy5 conjugates from the Universal Linkage System (ULS, Kreatech Biotechnology; ref. 30). Labeling efficiency for ULS-Cy3 and ULS-Cy5 were calculated from A260 (DNA), A280 (protein), A550 (Cy3), and A649 (Cy5) after the removal of nonbound ULS, on a NanoDropsND-1000 spectrophotometer (NanoDrop Technologies).
PCR, sequencing, and mutational analysis
TP53 mutations were identified by temporal temperature gradient gel electrophoresis followed by DNA sequencing for exons 2 to 11 as described in ref. (31). PIK3CA mutations were assessed using PCR-based amplification for exons 9 and 20 and sequence analysis as previously described (32). PCR products were purified over a QIAquick spin column (QIAGEN) and were sequenced using the BigDye Terminator Cycle Sequencing kit (Applied Biosystems) and an ABI 3730 automated capillary sequencer. For all PCR products with sequence variants, both forward and reverse sequence reactions were repeated for confirmation.
Microarray hybridization and data preprocessing
As previously described (33), hybridizations were done on microarrays containing 3.5 k BAC/PAC-derived DNA segments covering the whole genome with an average spacing of 1 Mb. The whole library was spotted in triplicate on every slide (Code Link Activated Slides, Amersham Biosciences). Data processing of the scanned microarray slide included signal intensity measurement using the ImaGene Software (BioDiscovery, Inc.) followed by median print tip normalization. Intensity ratios (Cy5/Cy3) were log2 transformed and triplicate spot measurements were averaged. This resulted in a 3,277 × 68 data matrix used for the analysis (Supplementary Data File S1, log2 aCGH data).
Gain and loss of frequency
Based on copy number levels, the frequency of gain and loss for all BAC clones was calculated using fixed log2 ratio thresholds of 0.2 and −0.2, respectively.
Identification of associated aberrant regions
We used SIRAC (25) to identify the chromosomal regions that are associated with the classes defined by prognostic gene expression signatures or clinical and pathologic characteristics of the tumor. This is achieved by explicitly incorporating the labels (e.g., 70-gene poor prognosis versus 70-gene good prognosis signature, wound signature–activated versus wound signature quiescent, GGI 3 versus GGI 1, or clinicopathologic variables such as ER positive versus ER negative) in the analysis. More specifically, as a first step, a SAM analysis (34) is used to select DNA probes that discriminate between the classes of interest at a chosen false discovery rate. A relatively high false discovery rate cutoff of 0.05 was used. This approach was taken to ensure that as many genes, which are potentially associated with the labeling of interest, as possible are included in further analyses. These significant DNA probes will be called the “relevant” probes. An illustrative result is shown in Lai et al. (Fig. 1, step 1 in ref. 25). As the next step, a P value is attached to the number of relevant DNA probes in a given genomic region to quantify the probability that the observed outcome could be a chance event. To this end, we used the hypergeometric test. We used a Bonferroni correction for multiple testing. Regions of aberration with a corrected P value smaller than 0.05 were considered significantly enriched for genomic aberrations. This step was repeated for different window sizes to detect both small and large aberrations [see, Fig. 1, step 2 in Lai et al. 25). Finally, regions of aberrations were identified based on a consensus between the results of the different window sizes [regions; see Fig. 1, step 3 in Lai et al. (25)]. Importantly, no discretization, smoothing, or segmentation algorithms were applied to the aCGH data before the SIRAC analysis. This avoids the parameter optimization step these models usually require, significantly simplifying the analysis.
DNA dosage–sensitive genes within identified aCGH regions
We ordered both expression probes (oligonucleotides) and DNA copy number probes (BAC clones) by position as, assigned by National Center for Biotechnology Information-Build327
on the genome, to find genes with elevated expression and gain/amplification of their corresponding genomic region within the identified aberrant regions found by SIRAC. Only probes for which a genomic location was found were used, resulting in an aCGH data set of 2,952 BAC clones and a gene expression data set with 10,986 genes. The association between DNA copy number and mRNA gene expression levels was calculated using the Pearson correlation between genes and the relevant probes in the DNA copy number data, revealing DNA dosage–sensitive genes. The correlation P value is the P value associated with the correlation between the level of DNA copy number and level of gene expression (mRNA). The “t test” P value is the P value of the t test for each gene that quantifies the ability of that gene to distinguish between the classes of interest (for example good versus poor outcome classes) based on gene expression data. To ensure that as many genes, which are potentially associated with the labeling of interest, as possible are included in further analyses, we filtered the genomic regions as follows: (a) they had to be aberrant in at least 10% or more of the class of interest; (b) the BAC clone aberrant region had a Pearson correlation P value of <0.01 with the corresponding genes; and (c) each individual gene had a t test P value of <0.01 to distinguish between the classes of interest.Comparison and independent validation of results
We compared our results on aCGH and mRNA expression data from independent breast cancer samples published in earlier reports (17, 18, 20, 24). Second, we validated, in a series of 89 primary breast tumors, using SIRAC on aCGH and gene expression data, the regions and genes that were associated with a 70-gene prognosis signature and molecular subtype. All data were acquired from the Supplementary Data of Chin et al. (18). Preprocessing and normalization was done as described by Chin et al. (18). Because no exact mapping information was available for all clones in the Chin data set, we gave the clones a 3-bp length centered on the mapping position as supplied by Chin et al. (18). We removed all clones with >50% missing values. We interpolated the missing values using the averaged values of their two positional neighbors. Probes mapped to the same genomic region were averaged and represented as a single BAC clone. This resulted in 2,149 unique BAC clones. Gene expression data were also acquired from Chin et al. (18). Probes not mapping to a single ENSEMBL ID were removed; probes mapping to Y chromosome genes were removed. This resulted in 21,339 unique Affymetrix probe measurements.
Results
Sixty-eight tumors are well-balanced over clinical/pathologic data
In a series of 68 invasive breast carcinomas, we analyzed DNA copy number changes, their locations, and their relationship with clinical and pathologic characteristics, prognostic gene expression profiles, mutations in the TP53 and PI3K genes, and prognostic markers.
Of these 68 tumors, 38 (56%) were ER positive, 36 (53%) were PR positive, and 14 (21%) were HER2 positive; 36 (53%) tumors were histologic grade 3, 17 (25%) were grade 2, and 15 (22%) were grade 1. Seventeen (25%) tumors showed a mutation in one of the hotspots of the PIK3CA gene and 31 (46%) tumors showed a mutation in one of the examined exons (2–11) of the TP53 gene. Detailed clinical and pathologic data are shown in Table 1 and Supplementary Data File S2, clincopathologic data.
Clinicopathologic characteristics . | n (%) . |
---|---|
Age at diagnosis (y) | |
<41 | 15 (22.1%) |
41-50 | 43 (63.2%) |
>50 | 10 (14.7%) |
No. of positive lymph nodes | |
0 | 33 (48.5%) |
1-3 | 26 (38.2%) |
>3 | 9 (13.2%) |
Tumor size (cm) | |
<2 | 24 (35.3%) |
>2 | 44 (64.7%) |
Local recurrence | |
No | 63 (92.6%) |
Yes | 5 (7.4%) |
Regional recurrence | |
No | 62 (91.2%) |
Yes | 6 (8.8%) |
Metastases | |
No | 43 (63.2%) |
Yes | 25 (36.8%) |
Distant metastases as first event | |
No | 47 (69.1%) |
Yes | 21 (30.9%) |
Death | |
No | 44 (64.7%) |
Yes | 24 (35.3%) |
Chemotherapy | |
No | 43 (63.2%) |
Yes | 25 (36.8%) |
Hormonal therapy | |
No | 61 (89.7%) |
yes | 7 (10.3%) |
Surgery | |
Lumpectomy | 32 (47.1%) |
Mastectomy | 36 (52.9%) |
ER | |
Negative | 28 (41.2%) |
Positive | 38 (55.9%) |
Missing | 2 (2.9%) |
PR | |
Negative | 31 (45.6%) |
Positive | 36 (52.9%) |
Missing | 1 (1.5%) |
HER2 | |
Negative | 52 (76.5%) |
Positive | 14 (20.6%) |
Missing | 2 (2.9%) |
Histologic type | |
IDC | 66 (97.1) |
ILC | 2 (2.9) |
Histologic grade | |
Grade 1 | 15 (22.1) |
Grade 2 | 17 (25) |
Grade 3 | 36 (52.9) |
Vascular invasion | |
None | 46 (67.7) |
Moderate | 6 (8.8) |
Extensive | 16 (23.5) |
Lymphocytic infiltrate | |
None | 50 (73.5) |
Moderate | 8 (11.8) |
Extensive | 10 (14.7) |
TP53 mutation analysis | |
Wild-type | 35 (51.5) |
Mutation | 31 (45.6) |
Missing | 2 (2.9) |
PIK3CA mutation analysis | |
Wild-type | 50 (73.5) |
Mutation | 17 (25.0) |
Missing | 1 (1.5) |
Clinicopathologic characteristics . | n (%) . |
---|---|
Age at diagnosis (y) | |
<41 | 15 (22.1%) |
41-50 | 43 (63.2%) |
>50 | 10 (14.7%) |
No. of positive lymph nodes | |
0 | 33 (48.5%) |
1-3 | 26 (38.2%) |
>3 | 9 (13.2%) |
Tumor size (cm) | |
<2 | 24 (35.3%) |
>2 | 44 (64.7%) |
Local recurrence | |
No | 63 (92.6%) |
Yes | 5 (7.4%) |
Regional recurrence | |
No | 62 (91.2%) |
Yes | 6 (8.8%) |
Metastases | |
No | 43 (63.2%) |
Yes | 25 (36.8%) |
Distant metastases as first event | |
No | 47 (69.1%) |
Yes | 21 (30.9%) |
Death | |
No | 44 (64.7%) |
Yes | 24 (35.3%) |
Chemotherapy | |
No | 43 (63.2%) |
Yes | 25 (36.8%) |
Hormonal therapy | |
No | 61 (89.7%) |
yes | 7 (10.3%) |
Surgery | |
Lumpectomy | 32 (47.1%) |
Mastectomy | 36 (52.9%) |
ER | |
Negative | 28 (41.2%) |
Positive | 38 (55.9%) |
Missing | 2 (2.9%) |
PR | |
Negative | 31 (45.6%) |
Positive | 36 (52.9%) |
Missing | 1 (1.5%) |
HER2 | |
Negative | 52 (76.5%) |
Positive | 14 (20.6%) |
Missing | 2 (2.9%) |
Histologic type | |
IDC | 66 (97.1) |
ILC | 2 (2.9) |
Histologic grade | |
Grade 1 | 15 (22.1) |
Grade 2 | 17 (25) |
Grade 3 | 36 (52.9) |
Vascular invasion | |
None | 46 (67.7) |
Moderate | 6 (8.8) |
Extensive | 16 (23.5) |
Lymphocytic infiltrate | |
None | 50 (73.5) |
Moderate | 8 (11.8) |
Extensive | 10 (14.7) |
TP53 mutation analysis | |
Wild-type | 35 (51.5) |
Mutation | 31 (45.6) |
Missing | 2 (2.9) |
PIK3CA mutation analysis | |
Wild-type | 50 (73.5) |
Mutation | 17 (25.0) |
Missing | 1 (1.5) |
Based on the gene expression profiles of the 68 tumors and a set of molecular classifiers, we assigned the tumors to the following subgroups: (a) 70-gene signature: 44 (65%) tumors with a poor prognosis signature and 24 (35%) with a good prognosis profile; (b) molecular subtype: 21 (31%) basal like of which 17 were triple negative (ER, PR, and HER2 negative), 10 (15%) HER2 like, 21 (31%) Luminal A, 12 (17%), Luminal B, and 4 (6%) normal epithelial like (3, 4); (c) core serum response: 41 (60%) tumors with an activated wound signature and 27 tumors (40%) with a quiescent wound signature (6); (d) GGI: 43 (63%) tumors were identified as grade 3 tumors and 25 (37%) were identified as grade 1 (8).
Frequency plots in agreement with previous studies
Figure 1 displays a frequency plot summarizing the distribution of DNA copy number aberrations in all 68 tumors. A fixed threshold of log2 > 0.2 for gain (green) and log2 < −0.2 (red) for loss was used. The frequency plot indicated the number of times the log ratio of a specific clone exceeds one of the thresholds. We observed frequent DNA copy number gains at 1q, 8q, 17q (22%), and 16p and losses at 16q, 13q, 23q, and 8p. High-level amplifications (log2 ratio > 0.5) were detected at 1q, 8q, 11q, and 20q. These observations are consistent with previously reported aCGH-based breast cancer studies (17–21).
Expression-based filtered SIRAC regions
Subsequently, we used the SIRAC algorithm to identify genomic regions enriched for BAC clones and coregulated genes, for each of the clinical and pathologic characteristics, and identify labeling induced by the prognostic gene expression signatures. Expression-based filtered SIRAC regions for each of the identified subgroups defined by gene expression signatures or clinical and pathologic characteristic are shown in Figs. 2 and 3.
Clinical and pathologic characteristics
We identified several aberrations specifically associated with clinical and pathologic characteristics. ER-negative and high-grade tumors showed a higher frequency of gains and losses than ER-positive and low-grade tumors. ER-negative and grade 3 tumors both showed gains at 6p21.1-33 and losses at 4p15.31. In addition, gains at 10p14-15.1 and losses at 5q33.1, 10q23.33-24.2, and 12q13.12-13.2 were identified in ER-negative tumors, whereas grade 3 tumors showed additional gains at 3q21.3-22.2, 8q24.13-3, 8q22.1-3, 12p13.31-32, and 17q24.3-25.1 and losses at 8p12, 8p21.1-8p22, and 15q21.1-2. Low-grade tumors and ER-positive tumors did not show regions with copy number gains, but did show losses, including 16q12.2-16q13 and 16q22.1-23.3 (Fig. 3). We identified a loss at 4p15.32 in TP53 mutant tumors and losses at 16q22.1 and 16q23.1-3 in tumors with a PIK3CA mutation, but no specific additional alterations (Fig. 2).
To further investigate the relationship between clinicopathologic parameters and mutation status, we performed an analysis of PIK3CA-mutated tumors in ER-positive tumors (n = 13) versus PIK3CA wild-type tumors in ER-positive tumors (n = 25). This analysis revealed a single region on chromosome 16q deleted in PIK3CA-mutated ER-positive tumors.
Seventy-gene poor prognosis signature
We found that the tumors with a 70-gene poor prognosis gene signature were associated with gains at 3q26.33-27.1, 8q22.1-24.21, and 17q24.3-25.1. Two hundred seventy-one of the 10,986 genes were located in genomic regions, identified by SIRAC, which were associated with a 70-gene poor prognosis signature. We filtered the genes (see Materials and Methods) to ensure that as many genes, which are potentially associated with A 70-gene poor signature, as possible are included: this resulted in 25 of 271 genes that showed DNA copy number gain correlated with upregulated mRNA expression (Supplementary File S3; Table 2). Three of these 25 genes [CCNE2, LYRIC (both 8q22.1), and EXT1 at 8q24.11) appear in the 70-gene prognostic signature. Eleven genes [BIRC5, TK1, EVER1 (17q25.1), EIF4G1, POLR2H (3q27.1), LAPTM4B (8q22.1), ZNF706, WDSOF1 (both 8q22.3), CML66 (8q23.1), SQLE, and ATAD2 (8q24.13)] belong to the intrinsic gene list that defines the molecular subtypes (3–5), indicating their importance in breast cancer (Fig. 2). Loss at 14q31 correlated with downregulated mRNA expression in association with a 70-gene poor prognosis signature was only found in three genes; nevertheless, none of them passed our filtering criteria as described in Materials and Methods.
Cytoband . | BAC clone . | DNA copy number* . | Ref Seq . | Gene name . | mRNA expression† . | % of samples . | |||
---|---|---|---|---|---|---|---|---|---|
Poor prognosis . | Good prognosis . | Poor prognosis . | Good prognosis . | Poor prognosis . | Good prognosis . | ||||
8q24.13 | RP11-495D4 | 0.31 | 0.11 | NM_003129 | SQLE | 0.12 | −0.17 | 50.00 | 16.67 |
8q24.11 | RP11-17E16 | 0.26 | 0.12 | NM_000127 | EXT1 | 0.06 | −0.07 | 43.18 | 8.33 |
8q22.3 | RP11-21E8 | 0.23 | 0.14 | NM_015420 | DKFZP564O0463 | 0.05 | −0.10 | 29.55 | 8.33 |
8q23.1 | RP11-21E8 | 0.23 | 0.14 | AF283301 | CML66 | 0.10 | −0.07 | 29.55 | 8.33 |
8q24.21 | RP11-3O20 | 0.22 | 0.12 | AB042405 | MLZE | 0.01 | −0.09 | 29.55 | 4.17 |
17q24.3 | RP11-478P5 | 0.21 | 0.05 | AL833387 | SLC25A19 | 0.06 | −0.11 | 29.55 | 8.33 |
8q22.1 | RP11-24H3 | 0.18 | 0.08 | AK000745 | LYRIC | 0.03 | −0.11 | 25.00 | 0.00 |
8q22.3 | RP11-24H3 | 0.18 | 0.08 | NM_016096 | LOC51123 | 0.06 | −0.09 | 25.00 | 0.00 |
8q22.3 | RP11-24H3 | 0.18 | 0.08 | NM_003406 | YWHAZ | 0.05 | −0.08 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_004702 | CCNE2 | 0.09 | −0.29 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_018407 | LAPTM4B | 0.12 | −0.22 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_012415 | RAD54B | 0.08 | −0.12 | 25.00 | 0.00 |
8q24.13 | RP11-16G11 | 0.15 | 0.08 | NM_014109 | ATAD2 | 0.03 | −0.16 | 20.45 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_001168 | BIRC5 | 0.19 | −0.39 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_003258 | TK1 | 0.11 | −0.27 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_004762 | PSCD1 | 0.04 | −0.11 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_007267 | EVER1 | 0.02 | −0.13 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_003255 | TIMP2 | 0.02 | −0.15 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | BC025951 | PGS1 | 0.05 | −0.08 | 15.91 | 0.00 |
3q26.33 | RP11-682A21 | 0.10 | 0.01 | NM_004301 | BAF53A | 0.03 | −0.07 | 13.64 | 0.00 |
3q26.33 | RP11-682A21 | 0.10 | 0.01 | NM_020409 | MRPL47 | 0.07 | 0.00 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_002808 | PSMD2 | 0.04 | −0.07 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_004953 | EIF4G1 | 0.09 | 0.01 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_004068 | AP2M1 | 0.11 | −0.04 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | U37689 | POLR2H | 0.07 | −0.04 | 13.64 | 0.00 |
Cytoband . | BAC clone . | DNA copy number* . | Ref Seq . | Gene name . | mRNA expression† . | % of samples . | |||
---|---|---|---|---|---|---|---|---|---|
Poor prognosis . | Good prognosis . | Poor prognosis . | Good prognosis . | Poor prognosis . | Good prognosis . | ||||
8q24.13 | RP11-495D4 | 0.31 | 0.11 | NM_003129 | SQLE | 0.12 | −0.17 | 50.00 | 16.67 |
8q24.11 | RP11-17E16 | 0.26 | 0.12 | NM_000127 | EXT1 | 0.06 | −0.07 | 43.18 | 8.33 |
8q22.3 | RP11-21E8 | 0.23 | 0.14 | NM_015420 | DKFZP564O0463 | 0.05 | −0.10 | 29.55 | 8.33 |
8q23.1 | RP11-21E8 | 0.23 | 0.14 | AF283301 | CML66 | 0.10 | −0.07 | 29.55 | 8.33 |
8q24.21 | RP11-3O20 | 0.22 | 0.12 | AB042405 | MLZE | 0.01 | −0.09 | 29.55 | 4.17 |
17q24.3 | RP11-478P5 | 0.21 | 0.05 | AL833387 | SLC25A19 | 0.06 | −0.11 | 29.55 | 8.33 |
8q22.1 | RP11-24H3 | 0.18 | 0.08 | AK000745 | LYRIC | 0.03 | −0.11 | 25.00 | 0.00 |
8q22.3 | RP11-24H3 | 0.18 | 0.08 | NM_016096 | LOC51123 | 0.06 | −0.09 | 25.00 | 0.00 |
8q22.3 | RP11-24H3 | 0.18 | 0.08 | NM_003406 | YWHAZ | 0.05 | −0.08 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_004702 | CCNE2 | 0.09 | −0.29 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_018407 | LAPTM4B | 0.12 | −0.22 | 25.00 | 0.00 |
8q22.1 | RP11-320N21 | 0.17 | 0.06 | NM_012415 | RAD54B | 0.08 | −0.12 | 25.00 | 0.00 |
8q24.13 | RP11-16G11 | 0.15 | 0.08 | NM_014109 | ATAD2 | 0.03 | −0.16 | 20.45 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_001168 | BIRC5 | 0.19 | −0.39 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_003258 | TK1 | 0.11 | −0.27 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_004762 | PSCD1 | 0.04 | −0.11 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_007267 | EVER1 | 0.02 | −0.13 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | NM_003255 | TIMP2 | 0.02 | −0.15 | 15.91 | 0.00 |
17q25.1 | RP11-141D15 | 0.11 | −0.01 | BC025951 | PGS1 | 0.05 | −0.08 | 15.91 | 0.00 |
3q26.33 | RP11-682A21 | 0.10 | 0.01 | NM_004301 | BAF53A | 0.03 | −0.07 | 13.64 | 0.00 |
3q26.33 | RP11-682A21 | 0.10 | 0.01 | NM_020409 | MRPL47 | 0.07 | 0.00 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_002808 | PSMD2 | 0.04 | −0.07 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_004953 | EIF4G1 | 0.09 | 0.01 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | NM_004068 | AP2M1 | 0.11 | −0.04 | 13.64 | 0.00 |
3q27.1 | RP11-682A21 | 0.10 | 0.01 | U37689 | POLR2H | 0.07 | −0.04 | 13.64 | 0.00 |
NOTE: Poor prognosis/good prognosis: tumors with a 70-gene poor/good prognosis signature.
*Average DNA copy number (log2).
†Average mRNA expression (log2).
Seventy-gene good prognosis signature
Tumors with a 70-gene good prognosis gene signature were associated with losses at 16q12.1-16q13 and 16q22.1-24.1 (Fig. 2). One hundred eighty-four of 10,986 genes were located in genomic regions, identified by SIRAC, that were associated with a 70-gene good prognosis signature. Fourteen of these 184 (Supplementary Data File S4) showed DNA copy number loss correlated with downregulated mRNA expression. It was suggested by Naderi et al. (35) that one of the genes, BM039 (CENPN), is part of a “core” prognostic signature because it was also identified as one of the 231 genes that correlated with prognostic categories in the study leading to the 70-gene prognosis signature (1).
Core response signature
In tumors with an activated wound signature, we found a copy number gain of 8q24.21 (MYC), in agreement with the findings of Adler et al. (23). However, we did not find a gain of 8q13.2, where CSN5 is located, but we did find a gain of the adjacent region (8q13.3). There were no specific DNA copy number losses associated with an activated wound signature (Fig. 2).
Genomic grade index
Tumors with a GGI of 3 showed gain at 3q26.33-3q27.3 and 8q22.1-8q24.22, and losses at 4p15.31, 8p12, and 8p21.1-2. Low-grade tumors showed regions with copy number gains at 16p11.2-13.13, and losses including 16q11.2-16q13 and 16q21.1-24.1 (Fig. 3).
Molecular subtypes
The previously identified molecular subtypes in breast cancer were also associated with specific DNA copy number changes (17).
HER2+
The HER2+ type tumors only showed gain on the q arm of chromosome 17 (17q11.1-17q12 and 17q21.31-23.2), covering the genomic position where the HER2 gene (17q12) is located. This is of course as expected, and the results suggest that this is the only aberration that differentiates this subtype from the other samples (Fig. 2). The region of 17q11.1-17q12 showed a median copy number log2 ratio of 0.82, and 2 of 22 genes within this region also showed high expression levels: the HER2 gene itself and GRB7, i.e., growth factor receptor–bound protein 7. This was expected because the HER2+ subtype is characterized by the amplification of both the HER2 and GRB7 gene, which is found to be coamplified and overexpressed with HER2 (36). Both genes were found to be both gained and overexpressed in 80% of the samples that were classified as HER2+ type.
Basal-like tumors
Basal-like tumors were associated with the largest number of aberrations. Copy number gains were found in chromosomal regions 6p12.3, 6p21.1-23, 8q24.21-22, and 10p12.33-10p14, and losses were found in regions including 4p15.31, 5q12.3-13.2, 5q33.1, 10q23.33, 12q13.13-12q13.3, 15q15.1, and 15q21.1.
Luminal A
For the Luminal A tumors, we found gain at 1q21.3-44 and 16p13.12-13, and losses at 16q11.2-16q13 and 16q22.1-24.1. Among the 277 genes associated with the Luminal A–type tumors, 40 genes (1q21.3-1qter, 16p13.12-13) showed gain correlated with upregulated gene expression, 38 of which were located at 1q21-44. KCTD3 and three others (RAB13, MUC1, and PPP2R5A) belonged to the intrinsic gene list that is used to define the molecular subtypes (3–5). Losses correlated with downregulated mRNA expression were found in 28 genes, all located at 16q; including NOC4, DC13, BM039, and CKLF, which belonged to the intrinsic gene list that originally defined the molecular subtypes (Fig. 3; refs. 3–5).
Luminal B
The Luminal B tumors had less pronounced aberrations with gains at 17q23.2, including genes TLK2, PSMC5, CCDC44, and SMARCD2, and losses at 1p31.3 and 8p21.2-23.1. When regions of loss and low expression were combined, Luminal B tumors showed a distinct loss at 8p21.2-8p23.1, including the following genes: MGC29816, DPYSL2, FLJ10569, XPO7, SH2D4A, LZTS1, PDLIM2, PDGFRL, TUSC3, C8orf16, LONRF1, and PRAGMIN (Fig. 2).
Discussion
Here, we analyzed DNA copy number changes within labeled data sets using SIRAC and integrated these findings with gene expression (mRNA) profiles in a series of 68 primary breast carcinomas. Our goal was to identify DNA copy number alterations that may be associated with prognostic gene expression profiles in breast cancer. The 68 breast tumors were labeled based on well-known clinical and pathologic features and gene expression signatures. These included the 70-gene prognosis signature (1), the wound signature (6), GGI (8), and the breast cancer molecular subtypes (3–5).
Using SIRAC, we identified associations between DNA copy number alterations, clinical and pathologic parameters, and prognostic gene signatures. We found that gains and losses varied between prognostic gene signatures and clinical and pathologic features. Gains and losses were more frequently found in tumors with unfavorable prognostic features (i.e., histologic grade 3, ER negative, HER2 positive, 70-gene poor prognosis, and basal-like type tumor). In particular, gains at 3q26-27, 6p12-22, 8q21-24, and 17q11-25 and losses at 4p15 and 5q33, 10q23, and 12q13 were associated with unfavorable prognosis. Another poor prognostic signature, high chromosomal instability signature (7)—was associated with gains at 3q26-27, 6p21-22, 8q21-24, and 12p13 and losses at 4p15, 8p12, 8p21-22, and 14q32.
Losses at 16q11-13 and 16q22-24 were found in tumors that share the following clinicopathologic characteristics: ER positive, histologic grade 1, PIK3CA mutated, and were associated with favorable prognosis signatures (70-gene good prognosis profile, GGI 1, chromosomal instability signature low, wound quiescent signature, and luminal A). This large concordance between DNA copy number alterations in these subgroups shows that there is a high degree of overlap in the patients identified by different breast cancer gene expression signatures and clinicopathological factors, as shown before (9, 37, 38).
Agreement between different breast cancer studies reconfirmed
We compared our results that are associated with the molecular subtypes of breast cancer (basal like, Luminal A, Luminal B, and HER2-like) with two recent reports of DNA copy number studies of these molecular subtypes of breast cancer (17, 24). We find that these studies have striking similarities regarding regions of aberrations in the basal-like, HER2-like, and the Luminal A–type tumors. Luminal A–type tumors seemed more homogeneous with respect to DNA copy number changes than the Luminal B–type tumor. The association between loss of 4p15 and 5q12-33 and gain at 6p12, 6p21, 10p12-13 and basal-like molecular subtype (defined by gene expression profiles) has now been found in at least three independent breast cancer aCGH data sets and exhibits more gains and losses and a lower frequency of loss at chromosome 16q than luminal-like breast tumors (Fig. 3). As shown here, the region of 6p21-23 is commonly amplified in basal-like tumors and harbors several oncogenes, including PIM1. Another gene identified as a DNA dosage–sensitive gene in basal-like breast cancer is VIM (vimentin) at 10p12.33, which is expressed in the myoepithelial layer of the glandular breast cells and, moreover, can be used to distinguish basal-like from luminal-type breast cancers. Interestingly, Andre et al. (22) found that triple-negative tumors showed gain at 6p21 and they identified at least three (VEGFA, PTEN, and EGFR) genes that showed high frequency of DNA copy number alterations, which was correlated with gene expression. In agreement with their findings, we find that basal-like breast carcinomas exhibit elevated expression and gain at 6p21.1 in 14% of the basal-like tumors and 17 tumors of 21 basal-like tumors are triple-negative tumors.
However, there are discordant regions between our and these two published studies that may be due to different patient selection criteria or differences in methods used to measure DNA copy number levels. We could only compare the Luminal A– and Luminal B–type tumor with the data of Bergamaschi et al. (17), as Adelaide et al. (24) did not subdivide the Luminal tumors into subtype A and B. Interestingly, the Luminal A–specific gain on chromosome 1q12-41 and 16p12-13 was also identified by Bergamaschi et al. (17) and (24) Adelaide et al. (24) also found strong association between these regions and luminal-type tumors. The HER2-like–type tumor with a known amplification at 17q12-21 was found in all independent data sets. Association between the loss of 5q21-25 and TP53-mutated tumors was not supported in our data set, as we only found a loss at 4p15.32.
Independent validation
We validated, in a series of 89 primary breast tumors, using SIRAC on aCGH and gene expression data, the regions and genes that were associated with a 70-gene prognosis signature and molecular subtypes. Together with matched array expression data, we identified common genomic regions in the two data sets [ours and Chin et al. (18)] showing strong coordinate expression changes.
No detectable ESR1 gene amplification
There are also several reports on the association between aCGH findings and clinical and pathologic characteristics. Most studies find that estrogen-negative breast tumors show loss of 5q21-35 and gain of 6p21-22 and estrogen-positive deletions at 16q (17, 19, 21, 24, 39, 40). In contrast, Holst et al. (41) reported that the ESR1 (6q25.1) gene is amplified in 21% of breast carcinomas. In view of our results and published aCGH data, we were surprised by this high frequency of ESR1 gene amplification/gain. To test whether we could correctly call a 2-fold gain spanning just 600 kb for each of the different array platforms, we carried out an experiment swapping aCGH data of equivalent-sized random X chromosome segments from a 48,XXXX versus 46,XX hybridization into normal female versus normal female data, thereby modeling 2-fold gain. Seventy-one of 100 randomly swapped segments were correctly called gained for the BAC array. However, given that the authors actually claimed that the 600-Kb ESR1 area was subjected to amplification, and not just a simply gain, the same analysis was run on 6× and 8× inserted into the 2× set. Eighty-five percent in the 6×/2× and 99.8% in the 8×/2× was called as amplified, indicating that the frequency of 5% for ESR1 copy number gain in breast cancer, as identified by us in several data sets, is correct and not flawed by technical problems in detecting ESR1 copy number gains.
None of the samples reported here harbored ESR1 gene amplification, and together with published aCGH data, we found only 12 of 533 (2.3%) breast tumor samples with amplification for the ESR1 gene and gain in an additional 14 (2.6%) breast tumors. It seems highly unlikely that the frequency of ESR1 gene amplification would be as high as 21%, but more likely less than 5% (42).
8q22-24 harbors potential prognostic markers
Most surprisingly, the cytogenetic regions of 8q22-24 were associated with more than one prognostic gene signature [70-gene poor prognosis (1)], GGI (8), wound signature (6), and molecular subtypes (3)]. Even more strikingly, CCNE2 (cycline E2) at 8q22.1 showed strong correlation between copy number and gene expression changes. CCNE2 is also one of three overlapping genes between the 76-gene prediction signature from Wang et al. (2) and the 70-gene signature from van't Veer et al. (1), both signatures being predictive signatures of distant metastasis-free survival in lymph node–negative patients. More recently, Hu et al. (43) showed that regional gain of 8q22.1 elevates expression of the metastasis gene metadherin (MTDH or LYRIC), which is overexpressed in >40% of breast cancers and is associated with poor clinical outcome. These results suggest that this DNA copy number information at 8q22.1 could act as potential prognostic markers in breast cancer.
In conclusion, using high throughput genetic techniques, it has been found that many breast cancers have increased or decreased copy numbers of large numbers of genes, and have a very high number of genes with high or low expression. The task of prioritizing them for analysis of their (biological and or) clinical significance has been extremely challenging. We have performed a combined analysis of DNA copy number changes and gene expression profiling to identify combinations of genes with potential biological relevance in breast cancer.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Acknowledgments
We thank Marleen Kok, Juliane Hannemann, Fabien Reyal, and Marc Bollet for critically reviewing the manuscript.
Grant Support: Supported by the Dutch Cancer Society (NKB2002-2575).
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.