GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers

1Cancer forms through the stepwise acquisition of somatic genetic alterations, including point mutations, copy-number changes, and fusion events, that affect the function of critical genes regulating cellular growth and survival. The identification of oncogenes and tumor suppressor genes being targeted by these alterations has greatly accelerated progress in both the understanding of cancer pathogenesis and the identification of novel therapeutic vulnerabilities. Genes targeted by somatic copy-number alterations (SCNAs), in particular, play central roles in oncogenesis and cancer therapy. Dramatic improvements in both array and sequencing platforms have enabled increasingly high-resolution characterization of the SCNAs present in thousands of cancer genomes.

2However, the discovery of new cancer genes being targeted by SCNAs is complicated by two fundamental challenges. First, somatic alterations are acquired at random during each cell division, only some of which ('driver' alterations) promote cancer development. Selectively neutral or weakly deleterious 'passenger' alterations may nonetheless become fixed whenever a subclone carrying such alterations acquires selectively beneficial mutations that promote clonal dominance. Second, SCNAs may simultaneously affect up to thousands of genes, but the selective benefits of driver alterations are likely to be mediated by only one or a few of these genes. For these reasons, additional analysis and experimentation is required to distinguish the drivers from the passengers, and to identify the genes they are likely to target.

3A common approach to identifying drivers is to study large collections of cancer samples, on the notion that regions containing driver events should be altered at higher frequencies than regions containing only passengers. For example, we developed an algorithm, GISTIC (Genomic Identification of Significant Targets in Cancer), that identifies likely driver SCNAs by evaluating the frequency and amplitude of observed events. GISTIC has been applied to multiple cancer types, including glioblastoma, lung adenocarcinoma, melanoma, colorectal carcinoma, hepatocellular carcinoma, ovarian carcinoma, medulloblastoma, and lung and esophageal squamous carcinoma, and has helped identify several new targets of amplifications (including NKX2-1, CDK8, VEGFA, SOX2, and MCL1 and BCL2L1) and deletions (EHMT1). Several additional algorithms for identifying likely driver SCNAs have also been described (reviewed in).

4Yet, several critical challenges have not yet been adequately addressed by any of the existing copy-number analysis tools. For example, we and others have shown that the abundance of SCNAs in human cancers varies according to their size, with chromosome-arm length SCNAs occurring much more frequently than SCNAs of slightly larger or smaller size. Therefore, analysis methods need to model complex cancer genomes that contain a mixture of SCNA types occurring at distinct background rates. Existing copy-number methods have also used ad hoc heuristics to define the genomic regions likely to harbor true cancer gene targets. The inability of these methods to provide a priori statistical confidence has been a major limitation in interpreting copy-number analyses, an important problem as end-users typically use these results to prioritize candidate genes for time-consuming validation experiments.

5Here we describe several methodological improvements to address these challenges, and validate the performance of the revised algorithms in both real and simulated datasets. We have incorporated these changes into a revised GISTIC pipeline, termed GISTIC 2.0.

6Cancer copy-number analyses can be divided into five discrete steps (Figure 1): 1) accurately defining the copy-number profile of each cancer sample; 2) identifying the SCNAs that most likely gave rise to these overall profiles and estimating their background rates of formation; 3) scoring the SCNAs in each region according to their likelihood of occurring by chance; 4) defining the independent genomic regions undergoing statistically significant levels of SCNA; and 5) identifying the likely gene target(s) of each significantly altered region. Figure 1 depicts a schematic overview of this process, highlighting the specific methodological improvements we will address in the present manuscript.

7The first step, accurately defining the copy-number profile of each cancer sample, has been addressed by multiple previous studies and is not discussed in detail here. We assume that segmented copy-number profiles have been obtained for all samples and all germline copy-number variations (CNVs) have been removed, yielding profiles of somatic events. The following sections describe improvements to steps 2 to 5. We evaluate these improvements on a test set of 178 glioblastoma multiforme (GBM) cancer DNAs hybridized to the Affymetrix Single Nucleotide Polymorphism (SNP) 6.0 array as part of The Cancer Genome Atlas (TCGA) project (the 'TCGA GBM set'), and on simulated data. Full technical details for each step are described in the Supplementary Methods (Additional file 1).

8Segmented copy number profiles represent the summed outcome of all the SCNAs that occurred during cancer development. Accurate modeling of the background rate of copy-number alteration requires analysis of the individual SCNAs. However, because SCNAs may overlap, it is impossible to directly infer the underlying events from the final segmented copy-number profile alone. Given certain assumptions about SCNA background rates, however, it is possible to estimate the likelihood of any given set of candidate SCNAs so as to select the most likely one.

9We have developed an algorithm ('Ziggurat Deconstruction' (ZD)) that deconstructs each segmented copy-number profile into its most likely set of underlying SCNAs (see Supplementary Methods in Additional file 1 and Supplementary Figure S1 in Additional file 2). ZD is an iterative optimization algorithm that alternatively estimates a background model for SCNA formation and then utilizes this model to determine the most likely deconstruction of each copy-number profile. Its output is a catalog of the individual SCNAs in each cancer sample, each with an assigned length and amplitude, that sum to generate the original segmented copy profile. We assume that most of these SCNAs are passengers, so that their distribution reflects, to a first approximation, the operation of the 'background' mutation process (see Supplementary Figure S2 in Additional file 3).

10A major advantage of the ZD method is its ability to separate arm-level and focal SCNAs explicitly by length. Prior studies have attempted to exclude arm-level SCNAs by setting high amplitude thresholds because, in contrast to focal SCNAs, few arm-level SCNAs reach high amplitude (Figure 2a). However, this approach suffers from at least two undesirable consequences: first, low- to moderate-amplitude focal copy-number events are eliminated from the analysis, reducing sensitivity to identify positively selected regions; and second, the amplitude threshold is left as a free parameter, allowing for potential over-fitting of the analysis to a desired result.

11We have previously shown that SCNA frequencies across cancers of diverse tissue origin are inversely proportional to SCNA lengths, with the striking exception of SCNAs exactly the length of a chromosome arm or whole chromosome (which are very frequent). This trend is preserved in the TCGA GBM samples (Figure 2b). This reproducible distribution provides a natural basis for classifying events as 'arm-level' and 'focal' based purely on length. Such length-based filtering of events allows for the computational reconstruction of 'arm-level' and 'focal' representations of the cancer genome (Figure 2c) and enables the inclusion of low- to moderate-amplitude focal copy-number events in the final analysis.

12To determine the benefits of this approach, we ran the original 'GISTIC 1.0' algorithm on the TCGA GBM set using three different thresholding approaches (Figure 3; Supplementary Table S1 in Additional file 4): 1) a low amplitude threshold (log2 ratio of ± 0.1) that only eliminates low-level artifactual segments; 2) a high amplitude threshold (log2 ratio of 0.848 and -0.737 for amplifications/deletions) used previously to eliminate arm-level events; and 3) the low amplitude threshold but also removing all SCNAs occupying more than 98% of a chromosome arm, leaving only the focal events.

13Filtering out arm-level events through use of either amplitude or length thresholds greatly increased the sensitivity of GISTIC for detecting focal amplifications and deletions (Figure 3; Supplementary Table S1 in Additional file 4). While entire chromosomes were scored as significant using only a low amplitude threshold, including gain of chromosome 7 and loss of chromosome 10 (Figure 3a), a number of recurrent focal alterations were missed, including amplifications surrounding CDK6, CCND2, and HMGA2. These alterations were detected using either the high amplitude (Figure 3b) or the focal length filters (Figure 3c).

14The benefits of length-based filtering result from the inclusion of low- to moderate-amplitude focal events. Amplification of PIK3CA and AKT1 and deletion of WWOX are detected using length-based filtering, but are not significant under the high amplitude filter (compare Figure 3b and 3c). Moreover, the length-based analysis identified significant SCNAs detected in neither of the amplitude-based analyses, including amplifications of MLLT10 and deletions of CDKN1B and NF1.

15No known GBM target gene was detected in either of the amplitude-based analyses that was not also detected by the length-based analysis. These results suggest that length-based filtering of arm-level events greatly improves the sensitivity of GISTIC to identify relevant regions of focal SCNA.

16We set out to define a scoring framework for SCNAs that more accurately reflects the background rates of alteration. Ideally, we aim to score each region of the genome according to the probability with which the observed set of SCNAs would occur by chance alone. Scores using this framework have a clear interpretation: the higher the score assigned to a region, the less likely that the SCNAs in that region are observed entirely by chance, and the more likely that they underwent positive selection.

17The probability of observing a single SCNA of given length and amplitude can be approximated by the frequency of occurrence of events of similar length and amplitude across the entire dataset (as in Supplementary Figure S2 in Additional file 3). However, since cancer genomes do contain drivers, this procedure is likely to overestimate the probability of observing SCNAs under the null model. Specifically, driver events tend to be shorter in length and of higher amplitude than passengers and therefore constitute the majority of events in their length/amplitude neighborhood (Supplementary Figure S3 in Additional file 5).

18To avoid biasing our background model, we set out to fit the log-probability distribution of SCNAs to a functional form that would be insensitive to the presence of driver events in the data (Supplementary Methods in Additional file 1). We made use of a large collection of 3,131 cancer samples run on the Affymetrix 250K StyI SNP Array plus several hundred additional samples run on the Affymetrix SNP6.0 Array (data not shown). At the level of resolution provided by these arrays, the probability of observing a focal SCNA at a given locus under the background model is roughly independent of length. As a result, the functional form for the log-probability distribution is similar to the original GISTIC G-score definition (G = Frequency × Amplitude), with the notable exception being that the new score is proportional to the amplitude in copy-number space rather than log-copy-number space.

19Although this functional form was empirically derived from a large collection of samples run on two different array-based platforms, it does lead to increased sensitivity to differences in dynamic range across platforms as well as differential saturation characteristics of probes within the same array platform. To avoid this problem, we routinely cap the segmented copy-number data at a level representing the signal intensity above which most probes start to saturate (Supplementary Methods in Additional file 1). This ensures that we are using data that originate from the linear regime of the probes' response curves and therefore are more comparable across platforms.

20As with GISTIC 1.0, we obtain P-values for each marker by comparing the score at each locus to a background score distribution generated by random permutation of the marker locations in each sample (Supplementary Methods in Additional file 1). This procedure controls for sample-specific variations in the rate of copy-number alteration. We correct the resulting P-values for multiple-hypothesis testing using the Benjamini-Hochberg false discovery rate method.

21Some genes are affected by non-overlapping deletions, either on different alleles in one sample or across multiple samples. For such genes, a marker-based score does not weight the presence of all deletions affecting that gene, despite the fact that these events are likely to have similarly deleterious effects on gene function. We have developed a modified scoring and permutation procedure, termed GeneGISTIC, that scores genes rather than markers (Supplementary Methods in Additional file 1). In each sample, we assign each gene the minimal copy number of any marker contained within that gene, and then sum across all samples to compute the gene score. Because genes covering more markers are more likely to achieve a more extreme value by chance, the permutation procedure is adjusted to account for gene size; the score for a gene covering n markers is compared against a size-specific null distribution generated by computing minima overall running windows of size n in each sample and then randomly permuting these minimal values across the genome.

22To determine the effect of gene-based scoring of deletions, we compared the results of gene-based and marker-based scoring on the TCGA GBM set (holding all other parameters equal). As expected, GeneGISTIC ranks known tumor suppressor genes higher and is more sensitive for genes subject to non-overlapping deletions (Supplementary Table S2 in Additional file 6). For example, RB1 was ranked 5th out of 39 regions using gene-based scoring (q-value = 2.6e-10) but only 13th out of 38 using marker-based scoring (q-value = 0.0013), and CDKN1B was ranked 26th using gene-based scoring (q-value = 0.08) compared to 38th using marker-based scoring (q-value = 0.19). NF1 was focally deleted in 12 of the 178 GBM samples (6.7%), and these deletions were frequently non-overlapping (Supplementary Figure S4a in Additional file 7). As a result, NF1 was scored just over or just under the significance threshold using the marker-based score, depending on the parameters used. By contrast, NF1 was robustly identified using gene-based scoring across all parameter combinations (Supplementary Table S2 in Additional file 6 and data not shown).

23However, because this scoring method does not score regions of the genome that are not in annotated genes, it could underweight or completely miss deletions occurring in non-genic regions. For example, in our GBM samples, gene-based scoring did not identify a region just outside of PCHD9 on chr13q21.3 that scored as highly significant (q-value = 4.4e-9) using the standard marker-based score (Supplementary Figure S4b in Additional file 7). While many non-genic deletions may in fact represent technical artifacts or rare germline events, some may be functionally relevant.

24Individual SCNAs, and indeed significantly amplified or deleted regions of the genome, may extend over more than one oncogene or tumor suppressor gene. Other significant regions may contain no oncogenes or tumor suppressor genes, but achieve apparent significance due to their proximity to a target gene. Thus, an additional step is required after genome-wide scoring to identify independently significant regions.

25GISTIC 1.0 solves this problem through the use of an iterative 'peel-off' algorithm, which greedily assigns all SCNAs to the maximal peak on each chromosome, removes them from the data, and rescores until no remaining region crosses the significance threshold. This approach reduces the power to identify secondary peaks that are close to previously identified significant regions (Figure 4a). However, since it is possible for individual SCNAs to affect multiple driver regions, a less greedy approach might identify additional peaks without significantly increasing the false discovery rate.

26We have, therefore, modified the method to allow SCNAs to contribute to more than one peak ('arbitrated peel-off'). We first greedily assign the entirety of an SCNA's score to the most significant peak it covers. In subsequent steps, however, we allow scores of previously assigned segments to be redistributed before deciding whether a putative region is significant (Supplementary Methods in Additional file 1). Like the original algorithm, the process terminates when no region has an adjusted score that exceeds the significance threshold. A similar modification of GISTIC has recently been proposed.

27Arbitrated peel-off is more sensitive than the original algorithm (Figure 4a; Supplementary Table S3 in Additional file 8). We generated 10,000 simulated datasets each consisting of 300 samples, with each chromosome containing a primary driver event in 10% of the samples and a secondary driver event in 5% of the samples. We analyzed the sensitivity of standard and arbitrated peel-off to detect the secondary peak as we varied the percentage of secondary driver events that overlapped the primary driver peak between 0% and 100% (Supplementary Methods in Additional file 1). At 0% overlap, the two methods were nearly equally sensitive at identifying the secondary peak. However, arbitrated peel-off was vastly more sensitive than standard peel-off as we increased the rate of overlap between primary and secondary peaks from 5 to 50% (Figure 4b), recovering an average of 2.4 times (range 1.2 to 3.8) more secondary peaks. Over 80% of the novel peaks identified by arbitrated peel-off corresponded to an actual simulated driver peak, demonstrating that the increased sensitivity is accompanied by high specificity.

28The primary and secondary peaks tend to merge when the overlap is above 50%, obscuring any appreciable difference between the two methods (Supplementary Figure S5 in Additional file 9). Indeed, neither method was capable of independently identifying the secondary peak once the percent overlap rose above 80%. These simulations demonstrate both the superior sensitivity of arbitrated peel-off as well as the challenge of identifying neighboring drivers.

29The final step in the GISTIC pipeline is to determine the region that is most likely to contain the gene or genes being targeted for each independently significant region of SCNA (the 'peak region'). The standard approach is to focus on the minimal common region (MCR) of overlap (Figure 5a), the region that is altered in the greatest number of samples and therefore would be expected to be the most likely to contain the target genes. However, one or more passenger SCNAs adjacent to, but not overlapping, the target gene can result in an MCR that does not include the true target. This is a frequent occurrence, especially when the frequency of the driver event is low (< 5%; Figure 5b). An alternative method (utilized by the GISTIC 1.0) is to apply a heuristic 'leave-k-out' procedure to define the boundaries of each peak region (Figure 5a). This procedure assumes that up to k passenger SCNAs (typically, k = 1) may aberrantly define each boundary of the peak region. While the 'leave-k-out' procedure correctly identifies the target gene more often than the MCR (Figure 5b), it suffers from the potential for overfitting introduced by the free parameter 'k'. Moreover, the accuracy of 'leave-k-out' varies depending on the number of samples and the frequency of the event under question. For fixed k, the sensitivity of 'leave-k-out' increases for increasing driver frequency (Figure 5b) and decreases for increasing sample size (Figure 5c).

30We developed a novel approach (termed 'RegBounder') to define the peak region boundaries in such a way that target genes would be included at a pre-determined confidence level, regardless of the event frequency or number of samples being studied (Figure 5a; Supplementary Methods in Additional file 1). RegBounder models the expected random fluctuation in G-scores within any given window size and uses this distribution to define a confidence region likely to contain the true driver at least γ% of the time, where γ is a desired confidence level. Unlike the MCR and 'leave-k-out' procedures, which are highly dependent on one or a few segment boundaries to define each region, RegBounder is designed to be relatively robust to random errors (either due to technical artifacts or passenger segments) in boundary assignment. When applied to real data, RegBounder captures known driver genes more effectively than 'leave-1-out' (and MCR) in regions with increased local noise (Figure 6a) and yet is capable of producing narrower boundaries than 'leave-1-out' in regions with little noise (Figure 6b).

31In simulated datasets, the performance of RegBounder was consistent across a wide range of driver SCNA frequencies (Figure 5b) and sample sizes (Figure 5c), and indeed controlled the probability of containing the driver. RegBounder captured the true driver gene in an average of 72%, 85%, and 95% of driver regions of varying frequency when run with a desired confidence level (γ) of 50, 75, and 95%, respectively. For no combination of sample-size, driver frequency, and γ did the average accuracy of RegBounder drop below γ.

32RegBounder also demonstrated a more optimal trade-off between peak region sensitivity (the likelihood of including the target gene) and specificity (the number of additional genes included) than the MCR or 'leave-k-out' approaches. The average size of the peak regions decreases with increasing driver frequency (Figure 7a) and sample size (Figure 7b) for all three approaches. However, RegBounder is more sensitive to these variables than the other methods, so that RegBounder peak regions (at 75% confidence) can range from an average of 90 times larger than the 'leave-k-out' peak regions (for datasets with few total driver events, in which the target gene locations are truly uncertain) to 37% smaller than the 'leave-k-out' procedure (for datasets with many total driver events). Thus, the increased confidence of RegBounder can even be achieved while producing narrower regions than the 'leave-k-out' procedure.

33RegBounder is also more consistent across datasets than the MCR and 'leave-k-out' methods. We randomly split the TCGA GBM set into two groups and compared the peak regions produced by RegBounder and the MCR and 'leave-k-out' procedures on each. Considering only those peaks that were identified by GISTIC in both datasets, only 23% of the MCRs and 31% of the 'leave-k-out' peak regions overlap between the two datasets, reflecting the low confidence with which these regions are assigned. By contrast, a majority (53%) of the RegBounder peak regions (at 75% confidence) overlapped, as expected (0.752 = 56%). This increased overlap came with only a modestly increased median size of the RegBounder peak regions (370 kb) compared to the leave-k-out (163 kb) or MCR (115 kb) peak regions.

34RegBounder regions are, on average, only 19% larger than the theoretically minimal peak region size for a wide range of driver frequencies (Figure 7c) and confidence levels (Supplementary Figure S6 in Additional file 10). These theoretically minimal peak region sizes were derived from the distribution of distances between the target gene and the MCR in our simulations (Supplementary Methods in Additional file 1). Our simulations reveal that RegBounder is capable of producing smaller peak regions than the 'leave-k-out' approach while simultaneously achieving greater target gene recall (compare Figures 5b and 7a; 'RegBounder 75%' versus 'leave-1-out', for driver frequencies > 5%). Thus, RegBounder is a robust algorithm for peak region boundary determination that demonstrates a more optimal trade-off between statistical confidence and peak resolution than previous heuristic approaches.

35The MATLAB source code for the GISTIC2.0 pipeline, along with a precompiled unix executable, will be available for download at. In addition, the entire pipeline can be accessed through the GenePattern analysis portal at.

36In addition to including all the methodological improvements described in this manuscript, the GISTIC2.0 source code has been designed to make efficient use of memory in storing segmented copy-number data (Supplementary Methods in Additional file 1). This improved memory efficiency should allow users with limited computational resources to run GISTIC2.0 on typical size datasets, and will be increasingly important for all users as the density of copy-number measuring platforms continues its rapid rise.

37We describe a number of analytical improvements to the standard copy-number analysis workflow that increase the sensitivity and specificity with which driver genes may be localized. We also demonstrate the utility of each of these changes using both simulated and real cancer copy-number datasets. While these changes have been specifically implemented in GISTIC 2.0, the challenges we describe apply broadly to the general task of identifying significantly aberrant regions of SCNA in cancer, and we anticipate that the approaches we have described can be adapted to other copy-number analysis workflows.

38The procedure we outline enables data-driven estimation of the background rates of SCNA and how these rates vary with features of the SCNA, such as length or amplitude. The specific trends we have observed are likely to depend on the resolution and characteristics of the measuring platform used to generate our datasets (the Affymetrix 250K StyI and SNP6.0 arrays). As more cancer samples are characterized using higher-resolution array and sequencing platforms, new trends are likely to emerge. Further improvements would account for such trends, possibly taking into account additional features that may determine SCNA background rates, such as the presence of known fragile sites of the genome or the surrounding sequence context. Indeed, we and others have recently shown that somatic deletions frequently occur in genes with large genomic footprints, suggesting the existence of a contextual bias in the rate of somatic deletion that is presently unaccounted for in our background mutation model. Our probabilistic scoring framework allows such trends to be accounted for once the background model has been specified.

39For the significant SCNAs, the background rate estimates also enable the delineation of regions likely to contain the target genes at predetermined confidence. RegBounder, the algorithm we devised to assign these boundaries, is more robust than either MCR- or 'leave-k-out'-based methods. RegBounder achieves this higher sensitivity by producing wider peak regions when the number of informative segments at a driver locus is small, but we find that RegBounder performs well compared to the theoretically optimal performance. However, RegBounder's underlying assumptions may not always be satisfied, including the assumption that each peak region contains a single dominant target gene and the expectation that copy-number breakpoints are independently distributed around the driver locus. To the extent that these assumptions are violated, RegBounder's performance may be worse than our simulations suggest.

40While the arbitrated peel-off approach described in this manuscript reflects a more sensitive way of identifying independently targeted regions of amplification and deletion than our prior approach, it is still an imperfect attempt to decipher the complexity of cancer copy-number alterations. One major limitation stems from the fact that array-based measurements map SCNAs onto a linear reference genome. However, many SCNAs are preceded by rearrangement events that juxtapose genomic regions separated by great physical distance in the germline (even different chromosomes). This level of detailed structural information is impossible to infer from probe-level copy-number estimates but can be obtained by sequencing paired-end libraries. Indeed, we anticipate that copy-number information derived from shotgun sequencing of cancer samples will become more common as sequencing costs continue to plummet. Tools for estimating and segmenting copy-number values from sequencing coverage data already exist, and these segmented copy-number profiles can, with only slight modification, be run through the GISTIC 2.0 workflow. Fully exploiting the level of detailed information provided by these technologies will, however, require a significant extension of the background mutation model to include the probability of random genomic rearrangements, as well as the ability to perform significance analysis, segment peel-off, and peak finding across non-contiguous regions of the reference genome. The data provided by these sequencing efforts should lead to new insights into the cellular and molecular processes underlying SCNA generation in different cancer types, and will allow for the development of vastly more detailed and accurate models of the background mutation rate of such events during tumor development.

41Full methods are available in the Supplementary Materials (Additional file 1).

42CNV: copy number variation; GBM: glioblastoma multiforme; GISTIC: Genomic Identification of Significant Targets in Cancer; MCR: minimal common region; SCNA: somatic copy number alteration; SNP: single nucleotide polymorphism; TCGA: The Cancer Genome Atlas; ZD: Ziggurat Deconstruction.

43RB and GG developed and coded the original GISTIC algorithm. CHM, SES, RB, and GG developed and coded the algorithmic modifications contained in GISTIC 2.0. CHM, MM, RB, and GG conceived and designed the present study. CHM, SES, and BH debugged and packaged the GISTIC 2.0 software release. CHM, MM, RB, and GG wrote the manuscript. All authors read and approved the final manuscript.