The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups

Abstract

The elucidation of breast cancer subgroups and their molecular drivers requires integrated views of the genome and transcriptome from representative numbers of patients. We present an integrated analysis of copy number and gene expression in a discovery and validation set of 997 and 995 primary breast tumours, respectively, with long-term clinical follow-up. Inherited variants (copy number variants and single nucleotide polymorphisms) and acquired somatic copy number aberrations (CNAs) were associated with expression in ~40% of genes, with the landscape dominated by cis- and trans-acting CNAs. By delineating expression outlier genes driven in cis by CNAs, we identified putative cancer genes, including deletions in PPP2R2A, MTAP and MAP2K4. Unsupervised analysis of paired DNA–RNA profiles revealed novel subgroups with distinct clinical outcomes, which reproduced in the validation cohort. These include a high-risk, oestrogen-receptor-positive 11q13/14 cis-acting subgroup and a favourable prognosis subgroup devoid of CNAs. Trans-acting aberration hotspots were found to modulate subgroup-specific gene networks, including a TCR deletion-mediated adaptive immune response in the ‘CNA-devoid’ subgroup and a basal-specific chromosome 5 deletion-associated mitotic network. Our results provide a novel molecular stratification of the breast cancer population, derived from the impact of somatic CNAs on the transcriptome.

Figures

Figure 1. Germline and somatic variants influence tumour expression architecture

a, Venn diagrams depict the relative contribution of SNPs, CNVs and CNAs to genome-wide, cis and trans tumour expression variation for significant expression associations (Šidák adjusted P-value ≤0.0001). b, Histograms illustrate the proportion of variance explained by the most significantly associated predictor for each predictor type, where several of the top associations are indicated.

Figure 2. Patterns of cis outlying expression refine putative breast cancer drivers

A genome-wide view of outlying expression coincident with extreme copy number events in the CNA landscape highlights putative driver genes, as indicated by the arrows and numbered regions. The frequency (absolute count) of cases exhibiting an outlying expression profile at regions across the genome is shown, as is the distribution across subgroups for several regions in the insets. High-level amplifications are indicated in red and homozygous deletions in blue. Red asterisks above the bar plots indicate significantly different observed distributions than expected based on the overall population frequency (χ2 test, P < 0.0001).

Figure 3. Trans-acting aberration hotspots modulate concerted molecular pathways

a, Manhattan plot illustrating cis and trans expression-associated copy number aberrations from the eQTL analysis (top panel). The matrix of significant predictor–expression associations (adjusted P-value ≤0.0001) exhibits strong off-diagonal patterns (middle panel), and the frequency of mRNAs associated with a particular copy number aberration further illuminates these trans-acting aberration hotspots (bottom panel). The directionality of the associations is indicated as follows: cis: positive, red; negative, pink; trans: positive, blue; negative, green. b, Enrichment map of immune response modules in the trans-associated TRA network, where letters in parentheses represent the source database as follows: b, NCI-PID BioCarta; c, cancer cell map; k, KEGG; n, NCI-PID curated pathways; p, PANTHER; r, Reactome.

Figure 4. The integrative subgroups have distinct copy number profiles

Genome-wide frequencies (F, proportion of cases) of somatic CNAs (y-axis, upper plot) and the subtype-specific association (−log10P-value) of aberrations (y-axis, bottom plot) based on a χ2 test of independence are shown for each of the 10 integrative clusters. Regions of copy number gain are indicated in red and regions of loss in blue in the frequency plot (upper plot). Subgroups were ordered by hierarchical clustering of their copy number profiles in the discovery cohort (n = 997). For the validation cohort (n = 995), samples were classified into each of the integrative clusters as described in the text. The number of cases in each subgroup (n) is indicated as is the in-group proportion (IGP) and associated P-value, as well as the distribution of PAM50 subtypes within each cluster.

Figure 5. The integrative subgroups have distinct clinical outcomes

a, Kaplan–Meier plot of disease-specific survival (truncated at 15 years) for the integrative subgroups in the discovery cohort. For each cluster, the number of samples at risk is indicated as well as the total number of deaths (in parentheses). b, 95% confidence intervals for the Cox proportional hazard ratios are illustrated for the discovery and validation cohort for selected values of key covariates, where each subgroup was compared against IntClust 3.