The molecular diversity of Luminal A breast tumors

Genomic data

Data from the TCGA study is accessible through the TCGA web portal at https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp. Data from METABRIC dataset was made available upon request and is now accessible through the European Genome–phenome Archive (EGA) with the accession number EGAS00000000083. Raw and pre-processed aCGH data for the Russnes et al. combined dataset can be accessed through the Gene Expression Omnibus (GEO) repository with accession numbers: GSE8757 [18], GSE20394 [19], and GSE19425 [20]. Mutation data for the Ellis et al. and Banerji et al. dataset are available as supplemental information within the respective publications. The datasets used in this manuscript can be explored using the cBio Cancer Genomics Portal at http://www.cbioportal.org/public-portal/ [21].

Recurrent mutations in breast cancer subtypes

In this work, we account for somatic mutations reported by multiple studies [10, 11, 16]. To estimate statistical significance of recurrent mutations across multiple datasets, we integrated the full set of reported mutations from these studies and used the binomial distribution to model the somatic mutation frequencies of genes. Given N samples, to estimate if a gene had a higher somatic mutation rate (mutations per nucleotide) than expected by chance, we evaluated if a gene with K observed non-silent mutations (summed over all tumor samples: $K={\sum }_{i=1}^N{k}_i$) and a global coding sequence of length L (summed over all tumor samples: L = l * N) had more mutations than expected, i.e., the average somatic mutation rate for all genes (G) with observed somatic mutations $p=\frac{1}{G}{\sum }_{i=1}^G\frac{K_i}{L_i}.$

Thus: $P(X\ge K)=1-P(X<K)=\sum _{i=0}^{K-1}\left(\begin{array}{c}L\\ i\end{array}\right)p^i{(1-p)}^{L-i}.$

From the set of p-values we estimated corresponding false discovery rates, or q-values, using the Benjamini–Hochberg procedure [22].

Copy number clustering

DNA copy number data were produced and processed for TCGA at the Broad Institute [10]. Briefly, copy number levels were inferred from Affymetrix SNP 6.0 CEL files by Birdseed and tangent normalization. Segmentation was then performed by Circular Binary Segmentation [23]. Copy number clustering was performed on normalized and segmented copy number data. A unified breakpoint profile (region by sample matrix) was derived by combining all breakpoints across all samples and determining the minimal common regions of change. Unified breakpoint profiles were computed using the R package CNTools [24] from Bioconductor [25] (http://bioc.ism.ac.jp/2.5/bioc/html/CNTools.html). Hierarchical clustering was done using the R function hclust, with manhattan distance and Ward’s agglomeration method [26].

Cluster centroids for each copy number cluster derived from the Luminal A TCGA dataset were computed by averaging cluster member features (unified breakpoints).

For each centroid we determined the most intrinsically variable breakpoints using the median absolute deviation (MAD) measure. Unified breakpoints with MAD > 0.1 were selected to define the cluster centroids. Average MAD values for this subset is threefold higher than the average of the remaining breakpoints, and twofold higher than the overall average (Fig. S1).

Samples from the METABRIC dataset have been assigned to the cluster whose centroid shows the highest correlation. Pearson correlation coefficient is a scale-independent measure and, thus, able to overcome the difficulty of comparing copy number values obtained with different platforms and on different scales. Nonetheless, we identified few samples (20 out of 594) with low correlation values with each centroids (pearson < 0.1). This subset is characterized by flat copy number spectra, such that all probes have identical values close to zero, with the exception of a few isolated probes likely to be either artifacts of the array or germline copy number variation (CNV). These flat samples have been assigned to the Copy Number Quiet subtype.

The quality of the clusters obtained with the centroids has been evaluated by different metrics: the clusters show similar proportions as those observed for the TCGA dataset (Fig. 1b), in-group proportion (IGP) values determined as in [27] are statistically significant (p < 0.001 for all clusters, Fig. 1b), and copy number pattern of alterations are remarkably similar to those observed in the TCGA dataset (Fig. S2).

Survival analysis

Survival analysis on the METABRIC dataset has been performed using the R package survival (http://cran.r-project.org/web/packages/survival/index.html) [26]. Patient follow-up has been limited to 15 years, and deaths related to other causes have been ignored in the analysis. The same analysis was performed for the Russnes et al. [19] dataset consisting of 77 Luminal A samples and combined three different cohorts: Ull cohort [19] (26 Luminal A samples), MicMa cohort [20] (9), and Chin cohort [18] (42).

Cox regression multivariate analysis was performed suing the

coxph

function from the

survival

R package. Multivariate analysis was used to assess dependencies between the classification induced by the Copy Number High (CNH) cluster and multiple covariates: tumor size, grade, stage, and age at diagnosis. This information was available for 468 Luminal A samples in the METABRIC dataset.

Subtype enrichment analysis

Given the tremendous heterogeneity displayed by Luminal A breast tumors, we analyzed the overall spectrum of genomic alterations across different subtypes looking for copy number subtype-specific patterns. Our approach relies on the general abstraction of gene alteration per sample, where each alteration belongs to one of the three categories: (a) gene is altered by mutations; (b) gene is primarily altered by CNA, and mRNA expression levels correlate with copy number changes; (c) “Wild-card” events (e.g., gene shows aberrant mRNA expression and/or methylation status independent of mutations and copy number).

These categories rely on two systematic approaches: for mutations we selectively analyzed the list of SMGs identified by the algorithm MuSiC [28], for copy number we analyzed frequently amplified and deleted region of interest (ROI) as identified by GISTIC [29]. We used the set of wide copy number gains and losses as the wild-card events. First, we selected the chromosome arms that were found to be recurrently gained or lost by GISTIC in [10]. Second, for each event we classified a sample as altered if segments accounting for at least 50 % of the whole chromosome arm length had values above (gain) or below (loss) selected thresholds. In this study, we used T = 0.15 as the absolute value for the threshold, where gains are defined by copy number level >T and losses by copy number levels <−T.

To systematically look for subtype-specific genomic events, we developed a method called Subtype-Enriched Alterations (SEA). Subtype enrichment is tested in two steps: (1) the distribution of alterations is compared to the expected distribution given in the number of samples that belong to each subtype by a goodness-of-fit test; (2) a hypergeometric p-value is derived for the subtype with the highest percentage of alterations when compared against all others. Enrichment p-values are then corrected for multiple testing [22]. Alterations in each category are tested separately and treated independently.

Differential mRNA expression analysis

To characterize the CNH subtype, beyond CNA and somatic mutations, we looked for genes differentially expressed between CNH tumors and the rest of Luminal A. We tested each gene by ANOVA, computed nominal p-values and corrected for multiple testing. High level amplification of the 8q region in the CNH group strongly influences this analysis with a high presence of genes located in this region. For this reason we used slightly less strict thresholds to select genes with nominal p-value < 0.05 and FDR-corrected q-value < 0.2.

Using this procedure we extracted two lists, one for genes that are up-regulated in CNH tumors, and one for genes that are down-regulated in the same set when compared to the other Luminal A tumors. Each list has been tested for functional enrichment using the DAVID Functional Annotation Tool [30].

Interestingly, despite a high number of genes from the highly amplified genomic region 8q21–24, only 3 out of 59 genes associated to “mitotic cell cycle” (GO: 0000278) are in this locus. Thus, mRNA up-regulation driven by this amplification does not seem to be related with mitosis regulation, spindle assembly, and chromatids segregation. We repeated the same functional enrichment test [30] after removing genes from 8q21–24, and separately for the up-regulated genes in this region to untangle potential different phenotypes associated to this copy number amplification and to other mRNA aberrations independent of it.

Up-regulated genes in CNH Luminal A that are not in the 8q21–24 region are similarly highly enriched for “mitosis cell cycle” (GO: 0000278), but also for related processes like “cell division” (GO: 0051301), “spindle” (GO: 000581), and “nuclear division” (GO: 0000280). These categories show elevated expression of both Aurora kinase A and B, as well as several genes in their pathways (e.g., PLK1, CDC25B, CCNA2, CDK1, INCENP, BIRC5, CDCA3/5/8, KIF2C). On the other hand, up-regulated genes in 8q21–24 were not found to be enriched for any functional annotation.

Mutual exclusivity analysis

All pairwise tests of mutual exclusivity were done using the switching permutation procedure described in [31]. This permutation strategy has the desirable property of preserving both number of alterations per gene and number of alterations per samples.

Mutual exclusivity modules (MEMo) were identified using the algorithm MEMo [31]. MEMo automatically identifies mutually exclusive alterations targeting frequently altered genes that are likely to belong to the same pathway. Genomic events were defined, as described in the previous section, following the “gene alteration per sample” abstraction and including focal copy number altered regions from GISTIC and recurrently mutated genes identified by MuSiC. Wild-card events included NF1 down-regulation (<1.5 standard deviations from the average) observed in 12 cases. The corresponding oncoprints for these modules are shown in Fig. S3 highlighting sample specific alterations in each module (samples are in the columns, altered genes in the rows). Module oncoprints were generated using the cBio Cancer Genomics Portal [21].

The landscape of Luminal A CNA

Luminal A tumors show great heterogeneity in terms of somatic mutations and CNA [10], indicating that additional substructure may be present within this large group. Dissecting the genetics of this tumor could be fundamental to inform therapeutics and predict clinical outcome.

We first explored the spectrum of copy number changes across Luminal A tumors, to identify novel and subset-specific alterations. We performed hierarchical clustering of Affymetrix 6.0 SNP copy number data from 209 Luminal A tumors from the TCGA dataset. Segments of uniform copy number value for each patient were compared to compute the set of unified breakpoints across the whole dataset, and the so determined set of minimal segments of change were used as features for the clustering procedure. Hierarchical clustering of copy number changes across the whole genome on the TCGA dataset revealed a complex structure of recurrent patterns of alterations. Based on clustering results and recurrent CNA, we were able to identify four major characteristic patterns and a mixed group (Fig. 2a; Table S1).

The first major pattern is characterized by 1q gain and 16q loss (1q/16q pattern, clusters a, b, and c). This pattern has been frequently observed in breast tumors and has been associated with the translocation der(1;16) [32, 33]. The 1q/16q pattern is dominant in cluster a, which features otherwise mostly diploid genomes. By contrast, cluster b is characterized by a broad deletion occurring on 6q, and cluster c has concurrent 11q13–14 focal amplification and 11q loss. High level amplification of the 11q13 and 11q14 loci is frequently observed in breast tumors and minimal regions of overlap target CCND1 and PAK1, respectively. Notably, these amplicons significantly co-occur with loss of the remaining part of the 11q arm across all tumors in the TCGA dataset (p = 6E−8, by one-tail-Fisher’s exact test).

Another group of Luminal A patients is characterized by a surprisingly quiet copy number spectrum (Copy Number Quiet pattern). These tumors (cluster d) have almost completely diploid genomes, with only few cases showing whole arm loss of 16q.

The third group includes clusters e, f, and g, and is strongly characterized by CNA of chromosome 8, with loss of 8p and gain of 8q (Chr8-associated pattern). Within this group, cluster f shows an interesting pattern of CNA, where 8p loss and 8q gain co-occur with 16p gain and 16q loss. These gains and losses affect the whole arms of the chromosomes, and in this group they do not co-occur with other CNAs. Cluster g, on the other hand, displays more copy number changes and is enriched for focal amplifications of 8p11.23–22 (FGFR1, ZNF703, and WSHC1L1), 8p11.21 (IKBKB), and 11q13–14 (Table S2).

The fourth group is characterized by the highest level of genomic instability among Luminal A tumors, including multiple focal CNAs (CNH pattern, cluster h). This group shows recurrent 20q gain, 5q loss, 8p loss, 8q gain, and is enriched for focal amplifications of the MYC oncogene on 8q24.21 (Table S2). Finally, the Mixed group is characterized by frequent whole-arm and whole-chromosome gains and losses, lacking the recognizable patterns seen in the other four groups (cluster i).

We validated our clusters using the 721 Luminal A samples from the METABRIC dataset [12]. We classified these tumors using centroids derived from the TCGA Luminal A dataset to identify clusters with the same alteration patterns (Fig. S2; Table S3). The clusters obtained from the METABRIC dataset occurred in similar proportions as the TCGA clusters and their quality was confirmed by the IGP [27] measure (Fig. 2b). Interestingly, the CN-Quiet and Chr.8-associated clusters show strong correspondences with two of the Luminal-enriched clusters from METABRIC and colleagues (IC4 and IC7, respectively), whereas components of the Mixed and CN-High groups are spread across multiple clusters (Fig. 2c).

The landscape of Luminal A somatic mutations

Distinct copy number patterns frequently come in tandem with equally variable landscapes of somatic mutations. Indeed, despite the lowest mutation rate among breast cancer subtypes, the Luminal A subtype shows the largest number of genes mutated with statistically significant recurrence (Fig. 3a; Table 2).

The most frequently mutated genes (>10 %) are PIK3CA, GATA3, MAP3K1, and TP53. Interestingly, all show significant associations with recurrent copy number patterns. PIK3CA and GATA3 mutations are mostly found in tumors with low CNA (CN-Quiet and 1q/16q), and in particular GATA3 mutations are enriched for the 1q/16q subgroup (p_GATA3 = 0.009). Notably, 9 out of 15 hotspot mutations for GATA3 are in this subgroup. MAP3K1 mutations are enriched in the Chr.8-associated subgroup (p_MAP3K1 = 4.22E−04). MAP3K1 mutations strongly co-occur with the 8p−/8q+/16p+/16q− pattern observed in cluster f (p_MAP3K1(f) = 1.9E−5), and thus, are largely mutually exclusive with focal amplification of 8p11.23–21 (Fig. 3b). All MAP3K1-mutated cases harbor at least one inactivating mutation, indicating loss of function, and most of them have at least two mutations, suggesting bi-allelic inactivation (Fig. 3c). Finally, the CNH subgroup shows an overall depletion for all recurrent Luminal A mutations, except for a strong presence of TP53 mutations (p_TP53 = 9.35E−6).

Luminal A tumors show also an interesting presence of hotspot mutations beyond PIK3CA and GATA3. AKT1 E17K activating mutations were observed in 3 out of 4 datasets and predominantly in Luminal A tumors (14 out of 20). Interestingly these mutations are perfectly mutually exclusive with those targeting PIK3CA. Additional hotspot mutations include those targeting KRAS (G12V/D), splicing factor SF3B1 (K666E, K700E), and c-Myc transcriptional repressor CTCF (R283C, H284P/Y/N).

Associations with clinical outcome reveal high-risk Luminal A subtype

The characterization of Luminal A tumors provided so far clearly identifies distinct subgroups within this tumor subtype, whose clinical relevance needs to be addressed now. The large dataset analyzed by Curtis et al. [12] has extensive clinical follow-up, enabling a reliable survival analysis. Kaplan–Meier analysis showed significantly different disease survival within the Luminal A subgroups (log-rank p = 0.015, Fig. S4). In particular, the CNH subgroup had significantly worse outcome (log-rank p = 4.6E−5, Fig. 4a) despite receiving similar treatments (Table S3). We validated the poor prognosis for the CNH tumors in an independent dataset consisting of 77 Luminal A tumors from three different cohorts [19] (log-rank p = 0.003, Fig. 4b; Table S4).

To assess dependencies between the CNH classification and other clinical covariates, we used Cox multivariate regression. We tested for dependencies for tumor grade, stage, size, and age at diagnosis. We found independent statistically significant association with outcome for tumor size (p = 6E−05), age at diagnosis (p = 0.004), and CNH classification (p = 0.01); no association was found between these covariates. The overall log-rank p-value for the combined covariates is p = 2E−08. The CNH classification showed, therefore, independent prognostic value (Table S5).

Finally, we compared scores derived from research-based versions of Oncotype DX [1], Mammaprint [34], and PAM50 Risk of Recurrence (ROR-S) [8] across Luminal A subgroups. We found that the CNH subgroup is consistently associated with a higher risk than the other Luminal A subgroups (Fig. S5), confirming the prediction of a worse prognosis.

As shown above, the CNH tumors are (1) characterized by high genomic instability, (2) depleted for typical Luminal A mutations (e.g., PIK3CA, 29 vs. 51 %, p_PIK3CA = 0.04), (3) highly enriched for TP53 mutations (48 vs. 7 % on average in the other Luminal A samples), and (4) tend to have MYC focal amplification, 20q gain and 5q loss (Fig. 4c). Interestingly, genes that are significantly over-expressed in CNH tumors compared to the rest of Luminal A cases are enriched for regulators of mitosis including Aurora kinases A and B, PLK1, Cyclin-A, Cyclin-E, CDK1, and Cdc25 (Fig. 4d; Table S6). A similar set of genes was previously found to be associated with 5q loss in Basal-like breast cancers [12]. Most of these genes have also been identified as “proliferation marker” genes [35]. The observed clinical outcome for CNH tumors, confirmed in two independent datasets, is therefore strongly associated to high genomic instability and proliferation as revealed by their molecular features.

Integrated pathway analysis reveals mechanisms of resistance to endocrine therapy

The great diversity of genomic lesions observed in Luminal A tumors maps to different cellular processes. To explore the role of these alterations in a pathway context, we used the MEMo algorithm [31], which identifies micro-pathways or modules whose components are frequently altered in a mutually exclusive manner. Statistically significant mutual exclusivity between recurrent alterations strengthens the hypothesis of functional relatedness and, more importantly, may reflect either functional redundancy, highlighting multiple ways to de-regulate the same pathway, or synthetic lethal interactions [31].

Modules extracted by MEMo in a Luminal A-only analysis highlight frequent alteration to the PI(3)K/Akt, MAP-kinase, and Ras/ERK signaling cascades (Figs. 5a, S3; Table S7). Alterations include PTEN inactivation, mutations of PIK3CA and AKT1, inactivation of MAP3K1 and MAP2K4, amplification of receptor tyrosine kinases (ERBB2 and IGF1R), and RAS activation either through activating mutations of KRAS or NF1 depletion by either DNA homozygous deletion or mRNA down-regulation (Fig. S6). Mutually exclusive inactivation of MAP3K1 and MAP2K4 was confirmed in the Ellis et al. [11] dataset, corroborating the hypothesis of reduced JNK signaling in Luminal A tumors and providing further insights into MAP3K1 mutations.

MEMo also identified a set of less frequent alterations targeting the N-Cor and SMRT complexes (Fig. 5b). While not statistically significant due to the small number of altered samples (28 out of 209 in total), alterations in these modules are almost completely mutually exclusive. Alterations include recurrent mutations targeting core components of the nuclear co-repressor complex (NCOR1, TBL1XR1, and GPS2) and the nuclear co-repressor 2 (NCOR2) or SMRT. Most of these mutations are either frame-shift or nonsense and thereby likely inactivating, they correlate with low mRNA expression, and frequently co-occur with hemizygous loss of the target gene. The modules also include homozygous deletions of ANKRD11, consistent with its ability to inhibit p160 steroid receptor co-activator recruitment [36, 37].

Co-repressor and co-activator complexes play a major role in regulating ER-α transcription and the inhibitory activity of Tamoxifen. Tamoxifen-bound ER has an increased affinity to co-repressors and specifically to N-Cor/SMRT. These co-repressors are required for the anti-proliferative effects of Tamoxifen (Fig. 4c). Repressing the NCor and SMRT complexes in human breast cancer cell lines turns Tamoxifen into an agonist of cell proliferation [38, 39], and lower or absent mRNA expression of NCor in patients correlates with shorter relapse [40, 41]. Here, for the first time, we identify distinct molecular mechanisms of inactivation of these complexes in patients. We directly link multiple genomic alterations, occurring in 13 % of the patients, to loss of co-repressors activity (Fig. 4d). These alterations may predict resistance to endocrine therapy.

Below is the link to the electronic supplementary material.