A mathematical framework to determine the temporal sequence of somatic genetic events in cancer

Abstract

Human cancer is caused by the accumulation of genetic alterations in cells. Of special importance are changes that occur early during malignant transformation because they may result in oncogene addiction and represent promising targets for therapeutic intervention. Here we describe a computational approach, called Retracing the Evolutionary Steps in Cancer (RESIC), to deduce the temporal sequence of genetic events during tumorigenesis from cross-sectional genomic data of tumors at their fully transformed stage. When applied to a dataset of 70 advanced colorectal cancers, our algorithm accurately predicts the sequence of APC, KRAS, and TP53 mutations previously defined by analyzing tumors at different stages of colon cancer formation. We further validate the method with glioblastoma and leukemia sample data and then apply it to complex integrated genomics databases, finding that high-level EGFR amplification appears to be a late event in primary glioblastomas. RESIC represents the first evolutionary mathematical approach to identify the temporal sequence of mutations driving tumorigenesis and may be useful to guide the validation of candidate genes emerging from cancer genome surveys.

Fig. 1.

Schematic diagram of RESIC. For cancer types with clinicopathologically defined stages (e.g., colorectal cancer), the temporal sequence in which genetic alterations arise during tumorigenesis can be identified by genotyping samples from patients at different stages of disease progression. For cancer types that are diagnosed de novo without detectable precursor lesions (e.g., primary GBM), the order of alterations cannot be identified with a similar approach. We present an evolutionary computational algorithm (RESIC) to identify the temporal sequence of events arising during tumorigenesis utilizing genomic data from a large number of samples (one per patient) of a particular histological type. In step 1, we use an algorithm such as GISTIC (4) to identify recurrent genetic aberrations in the genomics dataset. In step 2, these aberrations are ranked according to their pairwise association (statistically significant correlation, e.g., Fisher’s exact test). In step 3, the most likely sequence of these associated events is identified using RESIC. The results generated by RESIC are used to reconstruct the order in which alterations arise during development of a particular cancer type (step 4). Our methodology is applicable to large-scale datasets and can be used to identify the temporal sequence of many genetic alterations.

Fig. 2.

Evolutionary dynamics of genetic alterations leading to cancer. (A) Transition between mutational states and schematic representation of different evolutionary trajectories toward cancer. Initially, the population consists of N cells with genotype i and fitness (i.e., growth rate) r_i (detail). During each time step, a cell is chosen at random proportional to fitness to divide, and its daughter cell replaces another randomly chosen cell. During each cell division, a mutation arises with probability u. The mutated daughter has genotype j and fitness r_j. If r_j > r_i, the mutated daughter cell is advantageous as compared to the mother cell; if r_j < r_i, it is disadvantageous, and if r_j = r_i, it is selectively neutral. The probability that a mutated cell takes over the population is given by its fixation probability, ρ(r_i,r_j) = [1 - 1/(r_j/r_i)]/[1 - 1/(r_j/r_i)^N]. If r_j = r_i, then ρ(r_i,r_j) = 1/N. The transition rate between states i and j is given by Nuρ(r_i,r_j) in small populations. A population of wild-type cells may accumulate mutations in different orders; an example path from the unmutated population to a state with three mutations is highlighted in green. (B) Population dynamics. The dynamics of patients accumulating mutations is represented in this network where nodes (i.e., mutational states) are populated according to the transition rates from one mutational state to the next. In the example shown here, cells can accumulate two mutations. The number of patients harboring cells with no mutations are denoted by X₀, whereas those harboring mutations are denoted by X₁, X₂, and X₃. There is a constant influx of cases into the initial node. Cells in these patients accumulate mutations and populate the mutational states. The outflow from the fully mutated state eventually drives the system into steady state. An optimization algorithm is used to identify the transition rates for which the number of patients in each node at steady state coincides with the observations in a cross-sectional genomics dataset. The optimized parameter values of the evolutionary process serve to identify the most likely trajectory through the network.

Fig. 3.

Validation of RESIC utilizing colorectal cancer and glioblastoma data. We tested the predictions of RESIC in colorectal cancer and glioblastoma because the order of some events leading to those tumor types has been identified (10, 14). In the schematics of the networks, nodes represent the numbers of patients with a particular genotype, whereas black arrows represent transitions between mutational states. (A) The order of APC, KRAS, and TP53 in colorectal cancer. APC is shown in red, KRAS in green, and TP53 in blue. All mutation rates are 10^-7 per allele per cell division. We apply a pseudocount of 1 to the entire system to prevent states with zero observations. Schematics of the networks are shown at Left. We display the numbers and frequencies of patients in each mutational node in the network and show the most frequent paths through the network in the histograms at Right. Detailed results are listed in Table S2. (First Row) The APC–KRAS network. RESIC predicts that biallelic inactivation of APC likely occurs before any KRAS alteration. (Second Row) The APC–TP53 network. RESIC predicts that biallelic inactivation of APC likely occurs before TP53 inactivation. (Third Row) The KRAS–TP53 network. RESIC predicts that an alteration of KRAS likely occurs first. (B) The order of NF1 and TP53 in glioblastoma. NF1 is shown in orange. We study the mutational network of heterozygous alterations only since all NF1 and most TP53 mutations in the dataset are heterozygous (Table S1). A schematic of the network is shown at Left. Detailed results are listed in Table S2. We show the number of samples with each genotype observed in the complete set of 91 The Cancer Genome Atlas samples (black) and in the restricted set of 72 untreated samples (blue) (3). For both the unrestricted and the restricted sets, RESIC predicts that a TP53 point mutation likely occurs before NF1 is altered (Right).

Fig. 4.

Application of RESIC to secondary AML and primary GBM. (A) We analyzed a dataset of 57 patients with AML, including samples from two different disease states (MPN and post-MPN AML) from 14 different patients (19, 24). The data of AML patients were analyzed with RESIC (see B), whereas data from both MPN and AML patients were sequenced for JAK2 and TET2 (see C). (B) When applying RESIC to a set of secondary AML samples for which JAK2 and TET2 mutational status was known for both AML and MPN disease states, we found that JAK2 mutations likely precede TET2 mutations in this sample set. (C) Analysis of 14 patients for which samples were available from the MPN and AML disease states showed that TET2 mutations were present in the AML, but not preceding MPN sample, in 5 patients. (D) Genetic alterations in primary glioblastoma. The statistical significance of correlations between genetic alterations was calculated with Fisher’s exact test. Color codes range from 10^-20 (red) to 10⁰ (white). Significance after Bonferroni correction is marked at ∼3 × 10^-7 (orange). Note that lesions colocalized on the same chromosome have stronger correlations likely caused by large-region amplifications or deletions and thus cannot be considered as independent genetic events. (E) Prediction of RESIC for the PTEN-p16-EGFR network in primary glioblastoma. We show the frequencies of the initiating and final mutational events of this network. RESIC predicts that p16 deletion or EGFR low-level amplification are the most common initiating events with frequencies of about 35–39% each, whereas high-level amplification of EGFR is the most frequent last event of this network with a frequency of 56.4%.