The cancer genome hyperbola

At the top of this hierarchical classification, we identified two main tumor classes of similar size, each characterized by distinct sets of SFEs (Fig. 2a). Unexpectedly, although the distinction between copy number alterations and mutations was not used as a feature in our classification, these characteristic events were predominantly somatic mutations in one class and copy number alterations in the other (Fig. 2b). To reflect this trend, we named these two classes the M class (primarily with mutations) and the C class (primarily with copy number alterations), respectively. Notably, TP53 mutations were an exception to this trend, as they were strongly enriched in the C class (q = 3 × 10⁻¹⁷⁶), consistent with early mutations in TP53 causing copy number genomic instability (Supplementary Fig. 1). This division into two main tumor classes indicates that recurrent copy number alterations and mutations are predominant in different subsets of tumors.

Open in a separate window

Figure 2

The first partition of the pan-cancer data set identifies two main classes primarily characterized by either recurrent mutations (M class) or recurrent copy number alterations (C class). (a) Each class is composed of multiple tumor types in different proportions. (b) SFEs were tested for significant enrichment (more frequent than expected in a random distribution) in each class (events along the x axis, log-scaled q values on the y axis). Highly enriched events are primarily mutations in the M class and copy number alterations in the C class. Mut, mutation; meth, methylation change; amp, amplification; del, deletion. (c) The distribution of SFEs in tumors indicates that the number of copy number alterations in a sample (x axis) is approximately anticorrelated with the number of somatic mutations in a sample (y axis). The number of samples for a given (x,y) position range from 0 (white) to 243 (dark blue). CNAs, copy number alterations. Tumor types abbreviated as in Table 1.

Closer inspection of the distribution of selected functional events showed a striking inverse relationship between copy number alterations and somatic mutations at the extremes of genomic instability, particularly in highly altered tumors (Fig. 2c). Such tumors had either a large number of somatic mutations or a large number of copy number alterations, never both. We refer to this trend as the cancer genome hyperbola.

Tumors in the C class and M class were positioned along the two axes of this hyperbola (Supplementary Fig. 2). Whereas individual tumor types (defined by tissue of origin) had varying proportions of copy number alterations and mutations (Supplementary Fig. 3), none had high numbers of both.

We verified this approximately inverse relationship by adding 907 tumor samples from 6 additional tumor types to the pan-cancer set of 3,299 samples (Supplementary Fig. 4). In this larger data set, we also identified two major classes, one primarily dominated by mutations and the other primarily dominated by copy number alterations (Supplementary Fig. 4), with a remarkably similar set of characteristic functional events (Supplementary Fig. 4).

Starting from this first major subdivision, we applied the network modularity algorithm recursively to the C class and M class tumors and to their subclasses. The result was hierarchical division into several levels of subclasses characterized by distinct patterns of functional alteration at each level of granularity (Fig. 3, Supplementary Fig. 5 and Supplementary Table 3). We found that sample assignment to each subclass was robust in that it varied little upon systematic subsampling (Supplementary Fig. 6).

Open in a separate window

Figure 3

Characteristic patterns of functional alterations and distinct oncogenic processes as determinants of oncogenic signature classes (OSCs). (a) The first partition of the tree-like stratification (starting with ‘all tumors’ on the left) identifies two main classes: the M class (green) and the C class (red). We identify 17 oncogenic signature subclasses for the M class (M1–M17) and 14 oncogenic signature subclasses for the C class (C1–C14) (one row per subclass). (b) Each subclass includes subsets of tumors from several cancer types (grayscale heatmap; gray intensity represents the fraction of samples in a particular tumor type (column) and a particular subclass (row)). (c) Tree classification is determined at each level by sets of characteristic functional events (color intensity represents the fraction of samples in a subclass (row) affected by a particular functional event (column)). For functional copy number alterations, we indicate, if present, known oncogenes and tumor suppressors in parentheses, for example, 8q24 (MYC). (d,e) Subclass characteristic events reflect particular cellular processes (color intensity represents the fraction of samples in a subclass (row) affected by alterations to a particular process (column)) (d) and altered pathways involved in each of the processes (e). RTK, receptor tyrosine kinase; DSB, double-strand break. Tumor types abbreviated as in Table 1.

This classification highlights distinct mechanisms of oncogenesis as determinants of tumor subclasses, unexpected similarities between tumors originating in different tissues and new insights into alterations shared by multiple tumor types. Additionally, it provides a framework to explore therapeutic protocols on the basis of the genetic and epigenetic traits of tumors.

The M class

The M class of tumors included almost all the samples in kidney clear-cell carcinoma (KIRC), glioblastoma multiformae (GBM), acute myeloid leukemia (LAML), colorectal carcinoma (COADREAD) and uterine carcinoma (UCEC), with the exception of the serous-like subtype of UCEC. We identified 17 subclasses (M1–M17).

The first partition of the M class contained two main subclasses of mixed tumor type, which were characterized by distinct mutational events (Fig. 3a, Supplementary Fig. 5 and Supplementary Table 4). These subclasses had alterations in distinct oncogenic pathways, with alterations of phosphatidylinositol 3-kinase (PI3K)-AKT signaling characterizing the first main subclass (M1–M8) and with APC, TP53 and KRAS mutations most prominent in the second subclass (M9–M14).

Within the M class, we discovered recurrently mutated amino acids (hotspots) in the chromatin modifiers ARID1A and CTCF (Supplementary Fig. 7). ARID1A (Supplementary Fig. 7) is a member of the chromatin-remodeling complex SWI/SNF¹⁹ and, although truncating mutations in this gene have been reported in several tumor types²⁰, no recurrent hotspot had previously been identified.

CTCF encodes a chromatin-binding factor that acts as both a repressor and an activator of multiple genes, including known oncogenes and tumor suppressor genes (MYC, PLK, PIM1, CDKN2A and IGF2)²¹. CTCF achieves sequence-selective DNA binding by using different combinations of 11 zinc-finger domains (ZF1–ZF11)²². Mutations in CTCF were characteristic of subclass M5, which included several endometrioid tumors with microsatellite instability (MSI) and a small fraction of luminal A breast cancers (Supplementary Fig. 7). Mutations of CTCF affecting Arg448 have previously been reported^22,23 and occurred in multiple endometrioid tumors in subclass M5. Here we also identified seven mutations affecting residues upstream of ZF5 (Arg377 and Pro378), four mutations affecting ZF2 (His312 and Asn314), one of which targets one of the zinc-binding histidine residues (His312), and seven mutations affecting ZF1 (Gly261, Arg283 and His284), three of which affected the zinc-binding histidine residue His284 (Supplementary Fig. 7). Mutations observed in luminal A breast tumors specifically targeted ZF1, implying selective inactivation. We identified three splice-site mutations just upstream of exon 4, which encodes ZF1 and ZF2. One of these mutations caused an in-frame exon-skipping event (Supplementary Figs. 7 and 8). Even though the functional role of impaired CTCF activity in tumorigenesis is still unexplored, these mutations indicate that there is selection for specific zinc-finger loss and altered DNA-binding specificity that is not tumor type specific but broadly defines a subset of breast and endometrioid tumors.

Although most recurrent patterns of alteration characterize tissue-independent tumor subsets, subclasses M15–M17 were characterized by tumor type–specific mutational events (Supplementary Fig. 4); for example, EGFR amplification in GBM (M15), NPM1 mutation in LAML (M16) and VHL mutation in KIRC (M17). Our approach is therefore sensitive for reclassification both within and between tumor types.

The C class

The second major class was characterized primarily by TP53 mutations and multiple recurrent chromosomal gains and losses and is therefore called the C class. This class included almost all serous ovarian (OV) and breast (BRCA) carcinoma samples, as well as a large fraction of lung (LUSC) and head and neck (HNSC) squamous cell carcinomas and endometrioid tumors of the serous subtype (UCEC-serous).

Overall hierarchical subdivision of the C class led to a first major partition into two groups, primarily determined by the absence (subclasses C1–C6) or presence (subclasses C7–C14) of gains and losses on chromosome 8 (Fig. 3, Supplementary Fig. 4 and Supplementary Table 5).

Subclasses C3 and C4, which included a large fraction of LUSC and HNSC tumors, provided an interesting example of cross-cancer similarity, in which genomic alterations are shared by subsets of tumors of different origin. Subclass C3 was characterized by mutation of TP53 (92%), amplification of 3q26 (64%) and deletion of CDKN2A (32%); in contrast, subclass C4 had recurrent focal amplification of 11q13 (82%) where CCND1 is located. Some of these genomic differences actually converged on the same pathway, as loss of CDKN2A (C3) and gain of CCND1 (C4) both impair Rb-mediated cell cycle control.

Amplification of the 3q26 locus spans multiple genes, including PIK3CA and TERC. To identify candidate functional targets of this copy number alteration, we analyzed the mRNA levels of all genes in the 3q26 peak in amplified and diploid samples across all tumor types. Combined differential expression analysis identified ZNF639 as the most upregulated gene in the region (Supplementary Fig. 9). The zinc-finger protein ZNF639, also known as ZASC1 (zinc-finger protein amplified in squamous cancer 1), has previously been associated with the pathogenesis of oral and esophageal squamous cell carcinomas^24,25. PIK3CA was also found to be upregulated when amplified, whereas no correlation between mRNA levels and copy number was found for TERC.

The second major set of subclasses, C7–C14, had the highest degree of copy number alteration and was strongly characterized by recurrent gains and losses on chromosome 8, including amplification of 8q24 where the MYC oncogene is located. Amplification of MYC and somatic mutations in TP53 were the most frequent events in this subclass.

Cell cycle regulation and the DNA damage response were additional pathways affected by copy number alterations in subclasses C7–C14. The G1/S checkpoint was compromised by CCNE1 amplification in subclasses C7 and C11 and was bypassed by E2F3 amplification in subclass C13. Subclass C13 also appeared to have defective cell cycle arrest in response to DNA damage owing to inactivation of BRCA1 and BRCA2, which is recurrent in basal breast and ovarian tumors^4,10. Finally, subclass C14 had recurrent amplification and overexpression of the regulator of mitosis AURKA (encoding Aurora kinase A). Notably, these alterations were not specific for a single tumor type but rather characterized distinct subsets of tumors across multiple cancer types.

In summary, we found inactivation of TP53, MYC-driven proliferation and dysregulated cell cycle checkpoints as the hallmarks of the C class of tumors, which is dominated by recurrent copy number changes (Supplementary Table 6).

Filtered calls

The M class of tumors included several MSI and POLE-mutant cases, both of which have been associated with an unexpectedly high mutation rate^5,11. These types of tumors, therefore, have a large number of mutations that are probably not functional. To limit the number of likely non-functional events, we restricted our set of mutations in this class to all truncating mutations and to only nonsynonymous, single-residue substitutions that occurred at specific residues (hotspots). Hotspot residues were defined as recurrently mutated amino acids (represented by at least three mutations) or amino acids directly adjacent to a recurrently mutated one.

Similarly, the C class had a large number of copy number events that frequently spanned large chromosomal regions. The non-focal nature and high numbers of these alterations in the C class are likely to generate false positive assignments. To reduce this effect in the C class, assignment of copy number events in the GISTIC peaks was conditional on concordant mRNA expression changes.

Bipartite network modularity for recurrent genomic alterations

The SFEs naturally identify a network of relationships between samples and alterations. This network is a binary graph G = [(S,A), E], where nodes are either samples (S) or alterations (A) and edges (E) only connect samples to alterations. The problem of clustering tumors according to recurrent alterations can therefore be formulated as a graph clustering problem.

We addressed this problem using the notion of network modularity, originally introduced by Girvan and Newman³⁰. Given a partition of a graph G into distinct modules—subsets of nodes—the modularity associated with this partition is given by

where i is a module, e_ii is the fraction of edges with both ends in i, e_ij is the fraction of edges with one end in i and the other in j and r_i = Σ_je_ij is the expected fraction given the degree of the nodes in i.

This concept can be translated from a simple graph to a bipartite network. Recall that we defined our graph as composed of two sets of nodes: a set of samples S and a set of alterations A. Edges in our graph were defined as E = ((s,a) | s S, a A). Given a partition of the set of samples S, its modularity is the difference between the number of alterations shared by samples in the same module and the expected value of alterations. We defined the degree of each alteration a as d(a), equal to the number of samples connected to alteration a. Given a module m, d_m(a) is the degree of alteration a restricted to samples in module m. The e_ij term of the Girvan-Newman modularity can then be formulated as

where Z is a normalization factor. As with the Girvan-Newman modularity, given a partition of the set of samples S, the modularity measure defined above tells us how good this partition is in grouping together samples characterized by similar SFEs.

Modularity optimization by greedy partitioning

We adopted a greedy search procedure¹⁸ to optimize the modularity measure defined for our network of samples-to-alterations associations. This procedure starts by assigning each node or sample to a separate module and iteratively joining the pair of modules that produces the greatest increase in modularity. The approach is therefore similar to standard hierarchical agglomerative clustering. Although each step requires all (m²) pairs of modules to be scanned, the efficiency of this approach is derived from its requirement to compute, for each joined candidate pair, only the increase in modularity ΔQ. Note that ΔQ is given by

where e and r are intended to represent the corresponding sums over the set of modules and t is the iteration step.

In our network, upon joining modules m₁ and m₂ to form module m, we define e_m and r_m as

Therefore,

The algorithm stops when all nodes are grouped within the same module. The optimal partition is selected as the one with the highest modularity value among the ones generated through the optimization process. We used this optimization strategy to identify the optimal partition of our data set.

Hierarchical stratification of tumors by recursive modularity optimization

Community detection by network modularity optimization is limited by the size of modules³¹, and greedy partitioning tends to prefer incremental inclusions of single nodes in big modules rather than growing multiple modules simultaneously³². Moreover, the heterogeneous nature of our data set leads modularity optimization to be dominated by major differences between the main subclasses. The combination of these factors, although not affecting the optimal partition, limits the ability of this approach to capture the submodular structure of our data set at different levels of granularity. To address these issues, we recursively applied the greedy partition method. The algorithm proceeds as follows:

• Step 0: Determine the optimal partition P₀ of the whole data set.

• Step n: Determine the optimal partition $P_n^{(m)}$ for each module m contained in the partition $P_{n-1}^{(m*)}$, where m^* is the supermodule containing m.

At each step, the method subdivides a given module m by determining its optimal partition, i.e., the partition with maximal modularity. A module is no longer partitioned if (i) the modularity value of its optimal partition is below a limiting threshold Q_min and/or (ii) the module contains less than S_min samples. For this work, we set Q_min = 0.05 and S_min = 30. At each step, few small modules may be generated if the network contains isolated nodes or nodes with very few connections. In our data set, a minority of samples had few or no recurrent alterations. Small modules (with fewer than S_min samples), including samples with no or few uncharacteristic alterations, were ignored in the analysis. At each step, alterations were selected if they occurred in at least 1% of samples. The resulting set of partitions provides a hierarchical tree decomposition of the original data set, where the root is the whole set of samples and the leaves are modules that could no longer be partitioned.

In the analysis of oncogenic signature classes (OSCs) in Figure 3, we selected subclasses respecting conditions (i) and (ii) up to the third step of stratification. Exceptions included M1–M3 (fourth step), selected because of the marked and biologically relevant differences between these subclasses, and M15–M17 (second step), selected because each subclass was dominated by samples from a single tumor type.

Validation of the modularity optimization method

We tested our approach on two well-characterized data sets frequently used as benchmarks for network modularity detection. The first network is known as the Southern Women Event Participation network³³. It represents women’s attendance of social events in the Deep South, using data collected by Davis and colleagues in the 1930s to study social stratification. For this network, our approach was able to identify the two-module structure of the network (Supplementary Fig. 12) that coincides with the solution proposed by Guimera and colleagues³⁴ and, except for one woman, with the subjective solution proposed by the ethnographers that conducted the study.

The second test network is derived from data on corporate interlocks in Scotland in the twentieth century³⁵. The largest connected component of this network is composed of 131 directors and 86 firms. Our approach identified a nine-module solution with modularity value Q = 0.65 (Supplementary Table 8). The same component was analyzed by Barber³⁶ using an approach based on the eigenspectrum of the adjacency matrix of the network. In this work, the best solution obtained with the standard approach had Q = 0.566. A solution with Q = 0.66 was found using a modified version of their method that performed a search to optimize the number of modules rather than letting this number emerge from the modularity optimization procedure.

Enrichment analysis of genomic alterations

Each node of the tree, except for the leaves, represents a partition of a set of samples into separate modules or clusters. At each step, we identified the determinants (particular SFEs) of the partition by testing for statistically significant enrichment of each SFE in each class. For each SFE, we first tested for significant deviation from the expected distribution of its occurrences using a χ² test. Second, we selected the cluster with the highest fraction of samples altered by the particular SFE and tested for statistical enrichment by Fisher’s exact test. All P values were corrected for the false discovery rate (q value). SFEs listed in Figure 3 (middle) were selected as the most significantly enriched in each subclass at each branching of the tree decomposition (q < 0.001) or as the most frequent in each subclass.

Robustness of the subclasses

The robustness of the subclasses was assessed by removal of different percentages of samples and reclassification of the reduced data sets. During each run, hierarchical stratification obtained with the reduced data set was mapped to the original one by mapping each module from the reduced classification to the module from the original classification that maximizes the overlap (Jaccard) coefficient (J) associated with the two sets³⁷. Given a module m₁ and a module m₂, the J coefficient of m₁ and m₂ is defined as

with J = 1 if the two sets are identical and J = 0 if they are completely disjoint. Each mapping was scored with the average J obtained by the mapped modules. For each classification derived from a reduced data set, we derived a corresponding randomized version with the same hierarchical structure but permuted class memberships. Robust solutions were those with high average J values, averaged over repeated removal runs.

Reduced data sets were generated randomly by removing 5%, 20% and 50% of the samples (15 instances for each reduction), with the set with 50% fewer samples only used to evaluate the robustness of the M and C classes. We evaluated the robustness of subclasses separately at different levels of the hierarchical stratification and, for each evaluation, by estimating the expected J value using the randomized classifications.

We gratefully acknowledge contributions from the TCGA Research Network and its TCGA Pan-Cancer Analysis Working Group (contributing consortium members are listed in the Supplementary Note). The TCGA Pan-Cancer Analysis Working Group is coordinated by J.M. Stuart, C.S. and I. Shmulevich. We also thank E. Oricchio, X. Jing, S. Domcke, R. Sinha, J.J. Gao, G.B. Mills, J.J. Hsieh, B.S. Taylor, D.B. Solit, G. Rätsch, D.S. Marks and D. Bemis for discussions and/or critical reading of the manuscript. This work was supported by US National Cancer Institute funding of the TCGA Genome Data Analysis Center (U24 CA143840), US National Institutes of Health funding of Pathway Commons (U41 HG006623), by a Stand Up To Cancer Dream Team Translational Research Grant, a Program of the Entertainment Industry Foundation (SU2C-AACR-DT0209) and by US National Institutes of Health funding of the National Resource for Network Biology (P41 GM103504).