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

Model

RESIC is based on the principles of population genetics—the mathematical study of the dynamics of genetic variation within populations (6). Consider a population of N cells at risk of accumulating the genetic changes leading to cancer (Fig. 2A). Cells proliferate according to a stochastic process (7): At each time step, a cell is chosen proportional to fitness to produce a possibly mutated daughter cell. Subsequently another cell is chosen at random to die and is replaced by the newly produced cell to maintain homeostasis. A mutated cell can take over the population (i.e., reach fixation) or go extinct due to stochastic fluctuations (Fig. 2A, Inset). Depending on the order of appearance of particular mutations, the population of cells follows different evolutionary paths toward the fully mutated state (Fig. 2A). We developed a mathematical model describing the evolutionary dynamics of this system (see SI Text).

We assume that cancers originate from a single population of cells per person and study the evolutionary dynamics of individuals accumulating the mutations leading to cancer (Fig. 2B). We consider the dynamics of patients in steady state: There is a constant influx into the unmutated state, representing diagnosis of disease, and a constant outflux from the fully mutated state, accounting for deaths of patients or their cure. At steady state, the population is distributed across all possible states; this steady-state distribution can be compared to the numbers of clinical samples with the corresponding genotypes, where the total number of patients in a dataset is equal to the sum of patients in all states. This mapping is used to optimize a subset of parameters in the mathematical model (i.e., the fitness values of cell types) by minimizing the difference between the prediction and the observed frequencies in the dataset (see SI Text). Other parameters, such as cellular population size, mutation rate and influx rate, are estimated from experimental results (8, 9) and tested for robustness over several orders of magnitude (see SI Text). The output of RESIC is given as percent of the flux through the network via each particular path.

Results

The established sequence of genetic events arising during the multistep process of colorectal carcinogenesis (10) provides a unique example to test the ability of RESIC to recover this sequence from a cross-sectional dataset. We gathered genomic data from 70 advanced colorectal tumors (1, 2) (Table S1) and used RESIC to predict the temporal relationship between alterations in the APC, KRAS, and TP53 genes. RESIC predicts that the most likely sequence of events is homozygous inactivation of APC occurring before alterations of KRAS (Fig. 3A and Table S2). Similarly, we found that both APC alleles are likely mutated before TP53 is homozygously inactivated, and that at least one KRAS allele is likely mutated before inactivation of TP53 (Fig. 3A). Although the small number of samples prevents us from investigating the complete network of APC, KRAS, and TP53 mutations in a single computational analysis, the results of the separate smaller analyses can be combined into a sequence of events (see also Figs. S1 and S2) that coincides with the multistep model of colorectal cancer (10, 11). The results are very robust with regard to sampling stochasticity and variations of the population size, mutation rates, and influx values (Figs. S3–S6 and Table S3).

The temporal relationship between TP53 inactivation and mutations in the RAS pathway may not be the same in all cancer types. For example, genetic inactivation of neurofibromin-1 (NF1), a negative regulator of RAS activity (12), induces senescence of human astrocytes in the absence of TP53 mutations (13); furthermore, mice develop high-grade gliomas only if TP53 is inactivated prior to or simultaneously with NF1 (14). To investigate the temporal relationship between NF1 and TP53 inactivation in human primary glioblastoma (GBM), we applied RESIC to a genomic dataset containing sequence and copy number information regarding these two genes for 91 primary human GBM samples (3) (Table S1). We found that TP53 is likely inactivated before NF1 is lost (Fig. 3B). Because several GBMs in this particular sample set were collected after treatment with radiation and/or chemotherapy—therapies that might increase mutation rates (3), we reran RESIC using only GBMs that had not received prior therapy (n = 72). We again found that TP53 is likely mutated first (Fig. 3B). Our predictions are thus consistent with experimental data and suggest that TP53 inactivation may not only be an early genetic event in TP53-mutant secondary glioblastomas (15) but also in TP53-mutant primary glioblastomas.

We next applied RESIC to a leukemia dataset to test whether it also accurately predicts the order of mutations arising in “liquid” tumors. We investigated the order of mutations in the JAK2 and TET2 genes in a set of 57 secondary acute myelogenous leukemia (AML) samples transformed from a preexisting myeloproliferative neoplasm (MPN) (Fig. 4A). The classical MPNs are clonal disorders of hematopoiesis characterized by the presence of the JAK2V617F mutation in most patients as well as an increased likelihood of transformation to AML (16). Clonality studies suggest that acquired JAK2V617F mutations may not represent the earliest genetic event in MPN pathogenesis and are not required for transformation to AML (17, 18). Recent studies have identified mutations in the putative tumor suppressor gene TET2 in MPNs (17, 19–23), and clonal analysis suggested that TET2 mutations precede the acquisition of JAK2V617F mutations in MPN pathogenesis (19). However, these studies were done on a small number of patients. We therefore analyzed a larger set of patients with leukemic transformation from a preceding MPN, including samples from two different disease states (MPN and post-MPN AML) from 14 different patients (19, 24). When applying RESIC to the set of secondary AML samples for which JAK2 and TET2 mutational status was known, we found that JAK2 mutations likely preceded mutation in TET2 (Fig. 4B). Notably, 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 with concomitant JAK2/TET2 mutations (Fig. 4C) at the time of leukemic transformation. These data demonstrate, contrary to previous reports of smaller patient cohorts from a single time point, that TET2 mutations are more commonly acquired subsequent to JAK2 in MPN/AML pathogenesis (24), as predicted by RESIC.

After validation of RESIC using examples for which the temporal order of alterations is known, we analyzed a large integrated genomics dataset of primary glioblastoma samples (n = 594) (3, 4). These samples were used together with RESIC to identify genetic alterations occurring early during malignant transformation. Such alterations may result in “rewiring” of the signaling circuitry, confer a state of “addiction” to the new signal (25), and thus represent particularly promising targets for therapeutic intervention (26). We first identified areas of statistically significant gene copy number alterations using GISTIC (4). For each locus, we distinguished gain of a single copy from high-level amplification, and hemizygous from homozygous deletions (Table S1). We then identified alterations that are significantly positively correlated (Fig. 4D), because the determination of an order of oncogenic events is meaningful only for those events that co-occur sufficiently often. Although MDM2 and CDK4 on chromosome 12 and EGFR and MET on chromosome 7 were the most significantly correlated mutations, these associations are likely the result of large-region deletion or amplification events; because such changes cannot be attributed to independent mutational events, we excluded them from our analyses. Low- and high-level amplification of EGFR and homozygous loss of PTEN are the most significantly correlated events on separate chromosomes (p values < 10^-30 and 10^-24) (Fig. 4D). Because the dataset contains copy number information for 552 GBMs (excluding samples with homozygous loss of PTEN due to their relative infrequency) but sequence data for only 125 GBMs and because for both genes, copy number alterations and point mutations occur frequently, we first sought to investigate the robustness of the algorithm to the exclusion of sequence information. With only copy number data, we found that the two most frequent mutational paths through this network—initiating with EGFR low-level amplification and homozygous PTEN loss, respectively (Fig. S1)—have very similar frequencies. Therefore, RESIC identifies no clear order of events for this mutational network. Using both copy number and point mutation data, we found that EGFR biallelic alterations likely occur before PTEN loss (Fig. S1). Hence the only difference between the two analyses is loss of significance of the dominant path, suggesting that analyzing copy number information only may be a viable option.

To perform a computational analysis on a large mutational network, we next determined that p16 homozygous deletions frequently co-occur with EGFR and PTEN alterations (p value < 10^-8) in the dataset containing copy number alterations only (Fig. 4D). When studying the mutational network of EGFR, PTEN, and p16 (Figs. S2–4), RESIC predicts that the most common early alterations are EGFR low-level amplification and p16 deletion, with similar likelihood (Fig. 4E). Although there is no single most frequent path through the network (Fig. S3), the frequency of paths concluding with high-level amplification of EGFR is highest; the second most frequent final event is homozygous p16 deletion (Fig. 4E). These data suggest that glial progenitor cells may tolerate full EGFR activation only after inactivation of p16 or PTEN. This result agrees with the fact that EGFR overexpression is insufficient for tumorigenesis in mouse models of glioblastoma (27, 28), providing support for the temporal order of events predicted by RESIC.

Discussion

We have presented a rational methodology, RESIC, to identify the temporal order of oncogenic events during tumorigenesis. This computational pipeline uses cross-sectional genomic data as input and provides the most likely temporal order of genomic alterations as output. We validated the predictions of RESIC with events for which the temporal order is known and then applied it to a large integrated genomic dataset of primary glioblastoma samples. RESIC is based on three assumptions: (i) For each network under consideration, all cancers initiate without any of the mutations in the network. Hence all patients enter the network through the unmutated node, but might already have accumulated other mutations causing malignancy. This assumption ensures that the model is consistent with the use of cross-sectional data and is important for situations in which more than one set of mutually exclusive mutations lead to cancer; in such cases, tumor samples can be subdivided into classes of different tumor subtypes, and each subtype must be analyzed individually. (ii) The likelihood of diagnosis is uniform across all states. This assumption can be partially relaxed if more biological data are included in the model, which is equivalent to changes in the number of observed samples with each combination of mutations; variations in these numbers have negligible impact on the results (Fig. S5). (iii) The order of mutations arising during tumorigenesis can be inferred from the order of mutations arising after diagnosis. This assumption links the results obtained from a database containing diagnosed cancer samples to the processes of tumorigenesis before diagnosis. Interactions between tumor cells and the microenvironment or immune system can be considered as a modulation of mutation rates and fitness values but are excluded from the current implementation of RESIC for clarity. The results of RESIC are very robust to changes in the number of samples per genotype due to variability in sampling or the rate of diagnosis (Figs. S5 and S6), the population size of cells at risk of accumulating mutations, the influx value into the unmutated state, and the mutation rate (Table S3).

The frequencies of evolutionary paths through mutational networks identified by RESIC are never 100%. This finding may suggest that there is not a unique order in which these alterations occur—a hypothesis proposing that the genetic model of colorectal cancer (10) describes the order of mutations of only a subset of colorectal tumors (29). Alternatively, it may be attributed to the alteration of signaling pathways through alterations of mutually exclusive but functionally equivalent genes; it could also be the result of distinct mutational classes within the same cancer subtype. The latter issue is addressed by studying highly correlated mutations because strong correlation may reflect functional associations.

We have restricted our analysis to only a subset of mutations implicated in tumorigenesis because RESIC requires correlated alterations as input; the question of the sequence of alterations is meaningful only if those alterations co-occur sufficiently frequently. Even the largest currently available dataset, the The Cancer Genome Atlas data, combined with separate glioblastoma studies (n = 594), contains only a small number of associated lesions (Fig. 4D). These correlated lesions also reflect individual molecular subtypes of glioblastoma. Based upon our analyses, we expect that of the order of 100 samples are sufficient to analyze mutational networks containing several loci. However, RESIC can be applied to any tumor type by analyzing significantly correlated genetic alterations in separate computational analyses (see Fig. 3). Furthermore, RESIC can be applied to individual subtypes of cancers after these subtypes have been identified utilizing gene expression or other data.

Based on the examples presented here, we anticipate that our algorithm will provide the research community with a tool for the identification of tumor-initiating events using the emerging cross-sectional cancer genome datasets and will help with the generation of hypotheses about carcinogenesis that can be tested using modern mouse models of human cancer (30). The identification of the order of genetic alterations in specific cancer (sub)types may lead to important insight into cancer biology and should inspire studies aimed at elucidating how specific genes cooperate (or in a different order may prevent) tumorigenetic processes.

Supplementary Material