Background
Somatic copy number alterations (CNAs) are common genetic events in the development and progression of various human cancers, and significantly contribute to tumorigenesis 12. The coverage of CNAs in tumors varies from a few hundred to several million nucleotide bases, consisting of both deletions and amplifications with highly complex patterns 34. Recent advances in oligonucleotide-based single nucleotide polymorphism (SNP) arrays have made it possible to detect regional amplifications and deletions with high resolution on a genome-wide scale 56. A critical challenge in the genome-wide analysis of CNAs is to distinguish between the “driver” mutations that allow the tumor to initiate, grow, and persist, and the “passenger” mutations that represent random somatic events accumulated during tumorigenesis 137. Identification of these “driver” alterations can provide important insights into the cellular defects that cause cancer and suggest potential diagnostic, prognostic, and targeted therapeutic strategies 178.
By studying a sufficiently large collection of cancer samples, Significant Copy Number Aberrations (SCAs), defined as significantly recurrent CNAs that affect the same region in multiple tumors, are widely considered as informative surrogates of “driver” mutations that may help pinpoint novel cancer-causing genes 39. Past studies have detected many SCAs in a wide range of cancer types, with an impressive coverage of many known oncogenes and cancer suppressor genes 127. Several methods for finding regions of SCAs using CNAs data have been described in the literature, where the task of distinguishing between sporadic CNAs and SCAs is largely a statistical significance testing. Two reviews with qualitative comparison of different methods have been published 1011. Despite the use of different algorithms, a common theme in these methods is that they often adopt a four-step strategy: (1) detect CNAs and separate deletions and amplifications; (2) design and calculate ensemble test statistics associated with a genomic locus; (3) construct and/or estimate the probability distribution of test statistics under the null hypothesis; (4) perform multiple testing on a pool of genomic loci.
Significance testing for aberrant copy number (STAC) starts by converting the normalized log-ratios into a binary matrix, with zeros indicating no change and ones indicting losses and gains 12. STAC then proposes two statistics (footprint and frequency) to define regions of SCAs while adjusting for multiple comparisons, where the null hypothesis is that the detected CNAs from single-sample analysis are the realizations of random CNA placements whose probability distribution is generated by permutations on CNA segments 13. Genomic Identification of Significant Targets in Cancer (GISTIC) works on the real-valued step function of log-ratios that allows GISTIC to exploit both the type (amplification/deletion) and amplitude of CNAs 13. Using a semi-parametric permutation assuming independence between probes, GISTIC calculates a score that is based on both the amplitude and frequency of CNAs at each probe position and subsequently identify regions of SCAs, where amplification and deletion CNAs are handled separately, and armed-level and focal CNAs are further analyzed independently 14. Aimed to correlate information from neighboring probes with the amplitude and frequency of CNAs at each probe position, Kernel Convolution – a Statistical Method for Aberrant Regions detection (KC-SMART) uses varying-width kernel functions to calculate the testing statistics from the original log-ratios across multiple samples, producing the kernel smoothed estimate (KSE) at each locus by locally weighted regression 15. SCAs are selected based on a permutation-generated null distribution and Bonferroni correction. To substantially reduce computational burden in analyzing high-resolution and large-population data, correlation matrix diagonal segmentation (CMDS) identifies SCAs based on a between-chromosomal-site correlation analysis directly using the raw intensity ratios across all samples 16. CMDS uses a correlation statistics to detect SCAs with a standard normal null distribution whose parameters are estimated directly from the data and adjusts for multiple comparisons by false discovery rate.
Existing methods have several limitations. When working with unprocessed raw intensity ratios 131516, most methods are oblivious to noise clutter that can significantly confound estimation of the null distribution about true yet sporadic CNAs 917. Furthermore, these methods cannot distinguish between contributions of amplifications and deletions to the calculated overall test statistics that may affect the power to detect SCAs. While some effort has been made to incorporate correlation among neighboring probes into the test statistics, most methods assign a score to, and test the significance at, each individual probe locus 1415. In addition, while it is widely accepted that CNAs signals at adjacent probes are highly correlated 9131415, the assumption of probe independence is often adopted in constructing and learning the null distribution, probably for mathematical convenience 316. Moreover, existing permutation experiments using multiple samples cannot distinguish between the contributions of sporadic CNAs (obeying null distribution) and actual SCAs (deviating from null distribution) to the estimation of null distributions, resulting in theoretically conservative estimations especially when the number of true SCAs participating in the permutation is large.
We now report Significant Aberration in Cancer (SAIC), a carefully motivated method for accurately identifying SCAs using CNAs data from multiple samples. To distinguish between different biological roles of CNAs types and between noise and sporadic CNAs, we use discretized CNAs data and separately analyze copy number amplifications and deletions. By exploiting the intrinsic correlation among consecutive probes, we calculate and assign a score (test statistics) to each CNA unit instead of each single probe, based on both the amplitude and frequency of CNAs within the unit. To accurately estimate the null distribution governing sporadic CNAs, we perform random positional permutations on CNA units that preserve correlations inherent to the copy number data. More importantly, to minimize the unwanted participation of true SCAs in determining the null distribution 314, we iteratively detect SCAs and estimate an unbiased null distribution by an SCA-exclusive permutation scheme.
We tested SAIC on extensive simulation data sets, observing significantly improved performance with larger areas under the Receiver Operating Characteristics (ROC) curves and higher sensitivities at acceptable low false discovery rates, as compared to four popular peer methods (GISTIC, STAC, KC-SMART, and CMDS). We then applied SAIC to four real benchmark data sets, successfully identified the majority (84%) of previously reported SCAs harboring regions associated with well-known tumor-causing genes, and more importantly, detected some novel SCAs partially validated by the presence of known cancer-related genes.
Methods
Data format and definitions
Preprocessed log-ratio data are stored in a numeric N ×M matrix X. Each entry xnm represents DNA copy number (in log2-ratio) for sample n at probe m, where each row Xn corresponds to copy number for nth sample at M probes. Copy number amplifications and deletions are analyzed separately. We use the indicator function to divide matrix X into two matrices X = Xamplification +Xdeletion, where
X
amplification
=
I
x
nm
≥
θ
amplification
·
x
nm
,
X
deletion
=
I
x
nm
≤
θ
deletion
·
x
nm
,
with θamplification and θdeletion being the pre-specified thresholds. For brevity, we focus all subsequent discussion on Xamplification and make comments on Xdeletion when necessary.
Definition 1
Any copy number probe m whose associated copy number is amplified or deleted in at least one of N samples is called a CNA probe.
To exploit correlations inherent in copy number data, we first merge consecutive CNA probes into CNA regions, leaving the gaps consisting of only non CNA probes, see Figure 1. Within each CNA region, the Pearson correlation coefficient ρij between CNA probes i and j is then calculated for
i
≠
j
∈
M
:
ρ
ij
=
∑
n
=
1
N
(
x
ni
−
x
¯
i
)
(
x
nj
−
x
¯
j
)
(
N
−
1
)
s
i
s
j
,
where
x
¯
i
,
x
¯
j
, si and sj are the estimated means and standard deviations of copy numbers at probes i and j across N samples, respectively. If ρij is less than a pre-specified threshold θρ, a breakpoint occurs between probes i and j.
<p>Figure 1</p>An illustration on how CNA units are defined
An illustration on how CNA units are defined. Left: Consecutive CNA probes are merged into two intervals, with the first interval containing probes 1–10 and the second interval containing probes 14–16. Right: Each of the two intervals is split into CNA units according to the correction coefficients between CNA probes defined by Eq. (2), e.g., the first interval is split into three independent CNA units.
![]()
Definition 2
A sequence of consecutive CNA probes with no breakpoints is defined as a CNA unit, denoted by u (k, L) with k being the starting probe index and L being the length of the CNA unit.
Intuitively, a CNA unit consists of a sequence of highly correlated consecutive CNA probes. Figure 1 illustrates the concepts of CNA region and CNA unit, where two CNA regions contain 10 and 3 CNA probes, respectively, and the first CNA region is further split into three CNA units due to two breakpoints within the CNA region.
Summary statistics and significance assessment
Units that exhibit high or low average copy number are of interest, so it is natural to examine summary statistics for each unit. SAIC identifies significant aberration units through two steps. First, the method calculates a statistic (U score) that incorporates both the frequencies of occurrence and the amplitudes of the CNA probes within the unit, leading to the unit summary statistics given by
U
k
,
L
=
1
L
N
∑
n
=
1
N
∑
l
=
k
k
+
L
−
1
x
nl
.
Second, the method assesses the statistical significance of each CNA unit by comparing the observed statistic to the U scores that would be expected by chance.
Sporadic CNA units often occur throughout the genome, so a null distribution for Uk, L under the hypothesis that no SCAs are present, can be estimated by randomly permuting the overall pattern of presumed all-sporadic CNA units across the genome 391215. Though various permutation schemes can be adopted, due to different rates of CNA and different percentages of normal tissue contamination in tumor samples 18, permutation of CNA units across rows/samples should be avoided. As aforementioned, permutation should be performed on CNA units (instead of single CNA probes) that preserve correlations inherent to the copy number data, even if the CNA units are sporadic 3915. Another subtle but conveniently ignored issue is the different background rates of CNA units with varying lengths 1. Short CNA units occur at a frequency inversely related to their lengths and long CNA units occur approximately 30 times more frequently than would be expected by the inverse-length distribution. This observation is seen across all cancer types, is applicable to both copy gains and losses, and is supported by the calculated genome-average background rates for CNAs as a function of length 1. These considerations motivate our carefully designed SAIC permutation scheme.
Let L
denote the integer set containing the lengths of all the observed CNA units in X, K
denote the integer set containing the starting probe indices of all the observed CNA units in X, and X(t) be the random positional permutation of X for t = 1,2,…,T, with T being the total number of permutations. We now describe our method for analyzing CNA units for evidence of significant alteration in cancer, where we account for the difference in background rates between CNA units of different lengths by considering them adaptively.
Algorithm 1
Assessing the statistical significance of Uk, L
(1) Perform T random within-row positional permutations X(1), X(2), …, X(T) of the data matrix X on CNA units;
(2) Compute the value of summary statistic
U
k
,
L
X
t
for each permuted data set t = 1,2,…,T, and for each starting probe
k
=
1
,
2
,
…
,
M
−
L
+
1
and each length
L
∈
L
;
(3) Calculate and assign a P-value to each observed CNA unit u (k, L) for
k
∈
K
based on the extreme right-hand tail probability given by 919
P
U
k
,
L
X
=
1
+
∑
t
=
1
T
I
max
k
'
U
k
'
,
L
X
t
≥
U
k
,
L
X
T
+
1
,
where
I
·
is the indicator function.
The empirical P-values on Xdeletion are calculated by the extreme left-hand tail probabilities and reversing the inequality in Eq. (4). Both definitions produce P-values that are easy to interpret, and the “max” operation automatically adjusted P-values for multiple comparisons across CNA units thus controls the family-wise error rate 9.
In algorithm 1, it is important to note that when we generate a randomly permuted dataset based on the observed data, we do not re-define the CNA units but re-use the already-defined CNA units. Specifically, in each permutation, we randomly place the already-defined CNA units over the whole genome or each chromosome within each sample, and calculate the summary U score for each length of CNA units. Thus, independent of the unit length, the observed CNA units will always be retained (implicitly) in the permuted dataset. Moreover, when the number of permutations is sufficiently large, the p-values of observed CNA units can be accurately estimated. More precisely, to assess the p-value associated with an observed CNA unit of length L, we calculate the U scores for any consecutive L probes (probes do not need to reside within the same unit) across the genome, and compare the maximum score with the score of the observed CNA unit.
Iterative estimation of unbiased null distribution
One important issue concerning Algorithm 1 is the presence of true SCAs (departing from null distribution) in cancer genomes that presumably contribute high copy number deviations to the estimation of overall null distribution (governing only sporadic CNAs), potentially reducing power to detect less-extreme SCAs due to theoretical conservativeness 914. Loss of power is particularly critical in real-world applications where the number of true SCAs in cancer genomes may be large. Thus, to minimize the unwanted participation of true SCAs in determining the null distribution, we iteratively detect SCAs and estimate an unbiased null distribution by applying an SCA-exclusive permutation scheme. SAIC assesses the ‘new’ SCAs conditional on having found the ‘existing’ SCAs, successively correcting for true SCAs in order to better dissect and detect SCAs. Specifically, the CNA units associated with the ‘existing’ SCAs are masked as zeros after each iteration, resulting in a new data set X-SCAs in which already-detected SCAs becomes null.
Algorithm 2
Assessing iteratively the statistical significance of Uk, L
(1) Perform Algorithm 1;
(2) Check whether ‘new’ SCAs are detected. If ‘yes’, continue; if “no”, stop and re-calculate the P-values for all SCAs using truth converging null distribution;
(3) Mask the CNA units associated with newly detected SCAs as zeros and let
X
=
X
-
SCAs
, then go to step (1).
It has been shown experimentally that additional power to detect SCAs can be gained by removing the effect of newly detected SCAs after each iteration 9. However, an iterative SCA-exclusive permutation scheme raises another subtle yet critical issue concerning the convergence of null distribution learning and potential bias due to the expected false positive SCAs under the truth-converging null distribution. Fortunately, based on the careful design of Algorithm 2, the following theorem shows that, if we apply a significance level
α
'
=
α
/
1
+
α
where α is the targeted false positive rate (FPR), an unbiased estimation and detection results can be readily obtained using Algorithm 2 (see formal proof in Appendix A).
Theorem 1
Suppose that Algorithm 2 is used to iteratively detect SCAs and estimate truth converging null distribution. Let αbe the targeted FPR and
α
'
=
α
/
1
+
α
be the significance level used to detect SCAs. Then an unbiased truth converging null distribution can be obtained together with a theoretical FPR α.
SAIC algorithm and data preprocessing
Figure 2 shows the flowchart describing the entire SAIC algorithm. Our algorithm begins with two data preprocessing steps 18. First, the extracted raw copy number signals from CEL files are normalized using benchmark methods such as dChip (DNA-Chip Analyzer) 2021. Second, the normalized copy number signals are segmented into CNA regions using existing single-sample analysis methods such as CBS (Circular Binary Segmentation) 2223. The preprocessed log2-transformed ratios are subsequently analyzed by the novel algorithm described here.
<p>Figure 2</p>Schematic flowchart of combined SAIC algorithms 1 and 2
Schematic flowchart of combined SAIC algorithms 1 and 2.
![]()
Results
In the absence of definitive ground truth about the recurrent CNAs in the cancer genomes, the validation of a new method for detecting SCAs is always problematic 913161824. We first validate SAIC on multiple realistic simulation data sets and then proceed to evaluate the method using real CNA data sets. All data sets were analyzed according to the algorithm described in Figure 2. We tested SAIC and the four peer methods (GISTIC, STAC, KC-SMART, CMDS) on realistic simulation data sets. Comparative performance was based on the ground truth in terms of detection power 18 and the Receiver Operating Characteristics (ROC) curves 16. When applied to real CNA data, we compared and discussed biological plausibility of the implicated SCAs, and examined relative SCAs coverage between SAIC and GISTIC on benchmark data sets using Venn diagrams. To assure a meaningful and differential comparison, we emphasized experiment suitability when choosing algorithm parameter settings. For example, the algorithm parameter settings cannot be too “simple” (if there are only a few arm-level SCAs, all methods may perform equally well) or too “complex” (if there are many weak focal SCAs, no method will perform consistently well) 14.
Simulation studies
Multiple simulation data sets with definitive ground truth and various design or parameter settings were generated based on the modified benchmark models proposed in 9161824 and as used to assess various performance characteristics 91618. We first assessed the family-wise type 1 error rate (FWER) whose accuracy is crucial for methods that detect SCAs based on their P-values. If the FWER is either too conservative or too liberal, the P-value loses its intended meaning and does not reflect the actual false positive rate. Thus, we cannot control how many false positives are detected by setting a P-value based threshold 25. A large number of simulated null data sets (under the null hypothesis that no recurrent CNAs are present) were generated based on the realistic model proposed in 9 and subsequently analyzed with SAIC; results are presented in Table 1. Algorithm 2 was repeated 10,000 times, and the observed FWER was estimated by the proportion of at least one Uk, L (X) in X that was significant at α = 0.05 level 9. Values of the observed FWER in Table 1 (0.0497) suggest that SAIC is almost perfect when compared with slightly conservative values (0.0452) by similar method 9.
<p>Table 1</p>Empirical type 1 error rate for simulated data sets under the null hypothesis
Null simulation model
Empirical FWER atα= 0.05level
|
Copy number data
0.0488
|
Clumped copy number data (25%)
0.0500
|
Clumped copy number data (50%)
0.0493
|
Clumped copy number data (75%)
0.0505
|
We then assessed the detection power of SAIC as compared to GISTIC. Based on the simulation model proposed in 18, we generated 100 simulation data sets under each combinatorial parameter setting, resulting in a total of 1,900 simulation data sets, where each data set consists of N = 40 ~ 80 samples and each sample contains M = 5,000 probes. To replicate the effect of inevitable normal cell contamination 18, the copy numbers at every probes are simulated by a mixture of normal and tumor genomes, where the normal cell fraction λ is randomly drawn from a normal distribution
N
μ
λ
,
σ
λ
with μλ and σλ being the mean and standard deviation of normal cell fraction in the sample. Each sample contains two sporadic CNA regions, one deletion and one amplification randomly drawn from integer sets {0, 1} and {3, 4,…,8}, respectively. Each data set contains two recurrent CNA regions that are contributed from a fraction of samples according to a specified frequency ω, one deletion and one amplification similarly designed as aforementioned. The length of both sporadic and recurrent CNA regions is randomly assigned from 150 to 250 probes, realistically reflecting the estimated background rate of focal CNAs in a typical cancer sample genome 1. To equally assess the power in detecting deletion or amplification SCAs, we calculate the detection power of SAIC or GISTIC as the rate of successfully detecting inserted, deleted or amplified SCAs across 100 data sets. Table 2 summarizes the comparative detection power of SAIC and GISTIC for a total of 19 parameter settings across 1,900 data sets. These comparative experimental results consistently show that SAIC outperforms GISTIC with significantly increased detection power in 18 out of 19 simulations.
<p>Table 2</p>Power to detect SCAs by SAIC and GISTIC in simulation studies
N = 60,ω = 0.2,μλ = 0.6,σλ=
0.15
0.2
0.25
0.3
0.35
|
GISTIC
89%
86%
79%
74%
72%
|
SAIC
96%
94%
86%
86%
82%
|
N = 60, ω = 0.2, σλ = 0.25, μλ =
0.4
0.5
0.6
0.7
0.8
|
GISTIC
83%
81%
82%
72%
79%
|
SAIC
93%
91%
87%
79%
74%
|
ω = 0.2, σλ = 0.25, μλ = 0.6, N=
40
50
60
70
80
|
GISTIC
58%
73%
79%
86%
89%
|
SAIC
65%
83%
87%
93%
94%
|
N = 60, σλ = 0.25, μλ = 0.6, ω =
0.1
0.15
0.2
0.25
|
GISTIC
30%
58%
80%
92%
|
SAIC
37%
72%
87%
97%
|
We further assessed the overall performance of SAIC, measured by both sensitivity and specificity via ROC curves, as compared with the four peer methods (GISTIC, STAC, KC-SMART, CMDS). Based on the modified benchmark model proposed in 24, we generated 100 simulation data sets under each combinatorial parameter setting, where each data set consists of N = 50 samples and each sample contains M = 5,000 probes. The log-ratios at every probe are simulated by a mixture of normal and tumor genomes, with the normal cell fraction λ being randomly drawn from a uniform distribution
U
0.2
,
0.8
. Zero-mean Gaussian noise is then added to each sample with three levels of standard deviation σ randomly drawn from uniform distributions
U
0.2
,
0.4
U
0.4
,
0.6
, and
U
0.6
,
0.8
. To make the simulations more realistic, for each simulated sample genome, we insert 2 to 10 randomly located background CNA regions with the lengths varying from 10 to 50 probes. There are three ‘amplification’ (L = 30, 20, 10) and one ‘deletion’ (L = 20) ground truth SCAs embedded in each of the simulation data sets with a baseline frequency ω = 0.1. The copy numbers associated with amplification SCAs are 3, 4 and 5, and deletion SCAs are 0 and 1. In our simulation software, we use two parameters βL and βω to modify the length and frequency of these SCAs. Other parameter settings include θρ = 0.75, θamplification = 0.1 and θdeletion = −0.1 (default setting by GISTIC and CBS) for defining CNAs probes and units. Based on the estimated true positive rate (TPR) and corresponding FPR at different significance levels, Figure 3 presents ROC curves of SAIC and peer methods derived from the simulation studies. These comparative experimental results consistently show that SAIC outperforms the peer methods in terms of larger areas Az under the ROC curves or increased sensitivity at low FPR. More simulation studies are given in Additional file 1, where we report the power in detecting the boundaries of SCAs by these methods, and once again, showing outperformance of SAIC as compared to the peer methods 314.
<p>Additional file 1</p>Table S1. Comparative detection rates of ground truth SCA boundaries by STAC, GISTIC, KC-SMART, CMDS, and SAIC for simulation data sets under various model parameter settings. The results are calculated based on 100 replications for each of the parameter settings and using p-value (or q-value) cutoff threshold <0.05.
Click here for file
<p>Figure 3</p>Comparative performance of SAIC and four peer methods (STAC, GISTIC, KC-SMART, CMDS) on realistic simulation data sets, quantified by the partial ROC curves (north-west) (TPR: true positive rate; FPR: false positive rate)
Comparative performance of SAIC and four peer methods (STAC, GISTIC, KC-SMART, CMDS) on realistic simulation data sets, quantified by the partial ROC curves (north-west) (TPR: true positive rate; FPR: false positive rate). The results are the averages calculated based on 100 replications under each of various parameter settings.
![]()
Application to four real cancer copy number data sets
We applied SAIC to four real cancer copy number data sets and identified many SCAs that encompass established or potentially novel cancer ‘driver’ genes. The data sets are from ovarian cancer 2627, prostate cancer 818, lung adenocarcinoma 17, and glioblastoma 13. Due to their distinct biological functions in cancer development, SAIC analyzes separately chromosomes 1–22 and chromosome X/Y. To account for the different background CNA rates across chromosomes, we identify SCAs by performing SAIC on individual chromosomes. Other parameter settings include T = 1000 and α = 0.05 (theoretical significance level or FPR/FWER). To provide a somewhat independent verification, we compared the SCAs detected by SAIC with what reproduced by GISTIC on lung adenocarcinoma and glioblastoma data sets that have been previously reported 37.
Results on the ovarian cancer data set
Our in-house ovarian cancer data set consists of N = 63 tumor samples 262728. Copy number signals were acquired using the Affymetrix Human Mapping 250 K Sty SNP Array platform 1. Each sample contains a total of 238,230 probes across the whole genome. Other algorithm parameter settings include θρ = 0.95, θamplification = 0.263 (2.4 copies) and θdeletion = −0.322 (1.6 copies) 14. The genome-wide landscapes (via -log10P) of recurrent or sporadic CNAs observed in the data sets are given in Figure 4, where amplifications and deletions are separately shown (left and right sides). SAIC detected several SCAs (both amplification and deletion), many of which are biologically plausible and include known oncogenes (e.g., KRAS, CCNE1 and CCND2) and tumor suppressor genes (e.g., CDKN2A and CDKN2B) 26272930. Full lists of the genes covered by these SCAs are given in Additional file 2 (ST 2). SAIC also identified many other cancer driver genes within individual chromosomes (ST 3), such as SKIL, CDK4, PIK3CA, PTEN, FGD4, FGFR1.
<p>Additional file 2</p>Table S2 and Table S3. Details about the implicated SCAs and full list of genes covered by these SCAs, derived from the ovarian cancer data set.
Click here for file
<p>Figure 4</p>Genome-wide landscapes of recurrent or sporadic CNAs derived from 63 ovarian cancer samples
Genome-wide landscapes of recurrent or sporadic CNAs derived from 63 ovarian cancer samples. Amplifications and deletions are displayed on the left and right sides, separately, where dashed lines correspond to the significance level α = 0.05 for calling SCAs.
![]()
Results on the metastatic prostate cancer dataset
Our in-house prostate cancer data set consists of N = 55 clustered metastatic tumor samples, obtained from 13 prostate cancer patients. Copy number signals were acquired using Affymetrix Genome-Wide Human SNP Array 6.0 818. Each sample contains a total of 1,868,857 probes across the whole genome. To discount the potential bias due to imbalanced subject-cluster sampling 8, we chose to analyze the N = 13 representative samples and to detect global recurrent CNAs by SAIC. Other algorithm parameter settings include θρ = 0.95, θamplification = 0.263 and θdeletion = −0.322, the same as used in analyzing ovarian cancer data. The genome-wide landscape of recurrent or sporadic CNAs observed in metastatic prostate cancer data is given in Figure 5, where SAIC detected 15 amplification SCAs (318 genes) and 21 deletion SCAs (756genes). Full list of the genes covered by these SCAs are given in Additional file 3 (ST 4). Many of these genes are cancer related (e.g., EGFR, BRCA2, TP53, ATBF1, MYC and RB1). In individual chromosome analysis of the data set, SAIC identified many other SCAs involved with cancer driver genes, such as PTEN (ST 5).
<p>Additional file 3</p>Table S4 and Table S5. Details about the implicated SCAs and full list of genes covered by these SCAs, derived from the prostate cancer data set.
Click here for file
<p>Figure 5</p>Genome-wide landscapes of recurrent or sporadic CNAs derived from 13 metastatic prostate cancer samples
Genome-wide landscapes of recurrent or sporadic CNAs derived from 13 metastatic prostate cancer samples. Amplifications and deletions are displayed on the left and right sides, separately, where dashed lines correspond to the significance level α = 0.05 for calling SCAs.
![]()
Results on the lung adenocarcinoma and glioblastoma datasets
The lung adenocarcinoma data set consists of N = 371 tumor samples, publicly available at http://www.broad.mit.edu/cancer/pub/tsp7. Copy number signals were acquired using Affymetrix 250K Sty SNP Array, where each sample contains a total of 216,327 probes across the whole genome 7. To assure the general comparability of the results produced by SAIC and GISTIC, we adopted similar algorithm parameter settings used by GISTIC for detecting focal SCAs: θamplification = 0.848 and θdeletion = −1.15, in addition to θρ = 0.9. The genome-wide landscape of recurrent or sporadic CNAs observed in lung adenocarcinoma data is given in Figure 6, where SAIC detected 23 amplification SCAs and 26 deletion SCAs (after combining some of 98 recurrent CNAs within the same cytobands). Full list of the genes covered by these SCAs is given in Additional file 4 (ST 6). The Venn diagram in Figure 7 reveals the numbers of common and distinctive SCAs detected by SAIC and GISTIC. It can be seen that SAIC successfully detected most (87% amplification and 75% deletion regions) of the SCAs that have been detected by GISTIC, while also revealing many additional SCAs (10 amplification and 23 deletion regions) 7. In addition, the result from within-chromosome analysis of the data set is listed in Additional file 4 (ST 7).
<p>Additional file 4</p>Table 6 and Suplementary Table 7. Details about the implicated SCAs and full list of genes covered by these SCAs, derived from the lung adenocarcinoma data set.
Click here for file
<p>Figure 6</p>Genome-wide landscapes of recurrent or sporadic CNAs derived from 371 lung adenocarcinoma samples
Genome-wide landscapes of recurrent or sporadic CNAs derived from 371 lung adenocarcinoma samples. Amplifications and deletions are displayed on the left and right sides, separately, where dashed lines correspond to the significance level α = 0.05 for calling SCAs.
![]()
<p>Figure 7</p>Venn diagram on the numbers of common and distinct focal SCAs detected by SAIC and GISTIC in the lung adenocarcinoma samples
Venn diagram on the numbers of common and distinct focal SCAs detected by SAIC and GISTIC in the lung adenocarcinoma samples.
![]()
The glioblastoma data set consists of N = 141 tumor samples, publicly available at http://www.broad.mit.edu/cancer/pub/GISTIC, where each sample contains a total of 115,593 probes across the whole genome 3. Once again, we adopted the similar algorithm parameter settings used by GISTIC for detecting focal SCAs. The genome-wide landscape of recurrent or sporadic CNAs observed in glioblastoma data is given in Figure 8, where SAIC detected 15 amplification SCAs and 30 deletion SCAs (after combining some of 67 recurrent CNAs within the same cytobands). Full list of the genes covered by these SCAs are given in Additional file 5 (ST 8). The Venn diagram in Figure 9 reveals the numbers of common and distinctive SCAs detected by SAIC and GISTIC. It can be seen that SAIC successfully detected most (88% amplification and 75% deletion regions) of the SCAs that have been detected by GISTIC, while it also revealed many additional SCAs (8 amplification and 27 deletion regions) 3. In addition, the result from within-chromosome analysis of the data set is listed in Additional file 5 (ST 9).
<p>Additional file 5</p>Table S8 and Table S9. Details about the implicated SCAs and full list of genes covered by these SCAs, derived from the glioblastoma data set.
Click here for file
<p>Figure 8</p>Genome-wide landscapes of recurrent or sporadic CNAs derived from 141 glioblastoma samples
Genome-wide landscapes of recurrent or sporadic CNAs derived from 141 glioblastoma samples. Amplifications and deletions are displayed on the left and right sides, separately, where dashed lines correspond to the significance level α = 0.05 for calling SCAs.
![]()
<p>Figure 9</p>Venn diagram on the numbers of common and distinct focal SCAs detected by SAIC and GISTIC in the glioblastoma samples
Venn diagram on the numbers of common and distinct focal SCAs detected by SAIC and GISTIC in the glioblastoma samples.
![]()
The common SCAs regions (e.g., 7p11.2, 12p12.1, 9p21.3, etc.) are highly consistent with previous reports, and largely encompass well-known oncogenes or tumor suppressor genes. For example, EGFR (epidermal growth factor receptor) is an oncogene within 7p11.2 whose mutations or amplifications have been shown to contribute to uncontrolled cell division (a predisposition for cancer) 31. Many additional SCAs regions (e.g., 8p23.2, 21q22.2) contain or adjacent to disease-related genes (e.g., CSMD1 and TMPRSS3) and may warrant further study.
Discussion
SAIC is similar to many peer methods in that it assesses statistical significance of SCAs using a permutation-based null distribution 912141516. However, in contrast to the existing procedures, the CNA units used by SAIC preserve the essential correlation structures of serial probes whose estimated average correlation coefficient can be as high as 0.985 32. Moreover, by automatically adjusting P-values for multiple comparisons 3334 and iteratively re-estimating the null distribution exclusive of detected SCAs 9, SAIC can preserve the intrinsic false positive rate, without compromising detection power to resort to sometimes overly conservative schemes 3141516. Theoretic analysis and extensive experimental results show that SAIC preserves both type 1 error and detection power, see Tables 12. Furthermore, the novel concept of CNA unit and associated scoring and permutation scheme neatly parallels many considerations in the revised GISTIC2.0 14, for example, serial probes covering driver events should be more highly correlated than probes covering only passengers and thus more likely to identify the target genes. The flexible length-adaptive significance assessment of CNA units via Eq. (4) automatically accounts for distinct background rates according to their lengths and thus more likely to detect independent SCAs.
As for the θamplification and θdeletion parameters in the SAIC algorithm, there is no general guideline about how to select their values 14, since different types of cancers usually have different rates and magnitudes of background CNAs 142635. In addition, various degrees of normal cell contamination 18 and intratumor heterogeneity 3536 occur in many samples and these further complicate the selection of parameter values. Practically, lower thresholds were used to define broad (arm-level) CNAs while higher thresholds were used to define focal CNAs 314. A newly proposed strategy is to apply joint magnitude-length thresholds 14 and to correct normal cell contamination using BACOM 18. Since our main objective here is to identify focal CNAs, we have largely adopted the same strategy used in 314, i.e., we used relatively higher thresholds to define focal CNAs for subsequent analyses. Specifically, based on the observation that the magnitude of CNAs in ovarian and prostate cancers is relatively low, we used relatively lower and commonly used thresholds (2.0 ± 0.4), i.e., 2.4 copies for amplification and 1.6 copies for deletion. In contrast, on the datasets of lung adenocarcinoma and glioblastoma, we applied relatively higher thresholds (2.0 + 1.6, 2.0–1.1), i.e., 3.6 copies for amplification and 0.9 copies for deletion, that are similar to the thresholds used by GISTIC algorithms 314.
Similar situation occurs to the selection of θρ in defining CNA units 9. Lower values of θρ often produce longer CNA units while higher values of θρ often produce shorter CNA units. It has been reported that the average successive probe correlation of the segmented data can be as high as 0.985 932. In our experience in analyzing real cancer datasets, a value of
θ
ρ
taking between 0.7 and 0.95 would be a suitable choice.
It is important to note that the general conclusion on the relative performance of our SAIC and peer methods, at least based on the extensive simulation studies, remains largely true. We have used the same parameter values in all methods so that a fair comparison on their relative performances can be assured. Based on our analysis of real datasets using current parameter settings, it appears that SAIC performs well when compared to peer methods. In addition, the results of extensive simulation studies, performed under a variety of probe correlation schemes, show that SAIC preserves well the expected type 1 error, even when the probes follow non-stationary correlation structures similar to those found in real data 9.
SAIC currently can perform either genome-wide (except X/Y chromosome due to its distinct biological role) or chromosome-based CNA unit permutations. In the application of SAIC to real cancer data sets, we performed genome-wide, autosome-based, and X/Y-chromosome-based permutations. The combined results from using different permutation schemes contain more SCAs that may involve novel cancer driver genes. By exploiting the novel concepts of CNA probe, CNA unit, and multiscale permutation, experimental results show that SAIC can accurately detect the boundaries of SCAs with different lengths, see Additional file 1.
We have also performed simulation studies (data not shown) that indicate that detection power of SAIC can be further improved by correcting for normal tissue contamination using a recently developed BACOM method 18. However, the current version of BACOM requires paired tumor-normal sampling, availability of two-channel signals, and existence of deletion CNAs. Thus, we leave the combination of SAIC and BACOM as an extension for future research.
Conclusions
We have presented a novel approach to accurately detect significant recurrent CNAs in cancer genomes which is both statistically-principled and which, as illustrated by real examples, can be very effective at revealing SCAs within data. The concepts of CNA unit and iterative permutation are relatively simple to interpret, yet still convey considerable novel mathematical insights into data structure and bias correction.
It is worth noting that there are three novel features associated with SAIC. First, we define CNA unit to capture the intrinsic correlation structure in copy number data. Second, we perform iterative SCA-exclusive permutation to produce an unbiased null distribution. Third, we apply SAIC to real cancer copy number datasets and detect most previously reported SCAs covering well-known cancer genes.
Two important pending issues with the present algorithm are the expected significant impact of intratumor heterogeneity and normal cell contamination 183536. We are currently investigating applications of BACOM based normal cell correction 18 and hierarchical bi-clustering that optimize critical steps such as the selection of various thresholds and identification of subtype-specific copy number alterations.
Appendix A
Proof of theorem 1. Let α’ be the significance level used in each iteration to detect SCAs in Algorithm 2. Under the truth converging null distribution, we have
Pr
SCA
(
r
)
= 'yes'
|
SCA
(
r
−
1
)
= 'yes'
=
α
'
,
for iterations
r
=
1
,
2
,
…
,
∞
since SAIC assesses the ‘new’ SCAs at the rth iteration conditional on having found the ‘existing’ SCAs at the (r-1)th iteration.Considering
Pr
SCA
(
2
)
= 'yes'
=
Pr
SCA
(
2
)
= 'yes',SCA
(
1
)
= 'yes'
=
Pr
SCA
(
2
)
= 'yes'
|
SCA
(
1
)
= 'yes'
Pr
SCA
(
1
)
= 'yes'
=
α
'
·
α
'
=
α
'
2
.
Therefore for the rth iteration,
Pr
SCA
(
r
)
= 'yes'
=
Pr
SCA
(
r
)
= 'yes',SCA
(
r
−
1
)
= 'yes',
…
,SCA
(
1
)
= 'yes'
=
Pr
SCA
(
r
)
= 'yes'
|
SCA
(
r
−
1
)
= 'yes',SCA
(
r
−
2
)
= 'yes',
…,
SCA
(
1
)
= 'yes'
·
Pr
SCA
(
r
−
1
)
= 'yes'
|
SCA
(
r
−
2
)
= 'yes'
,
SCA
(
r
−
3
)
= 'yes',
…
,SCA
(
1
)
= 'yes
·
…
·
Pr
SCA
(
2
)
= 'yes'
|
SCA
(
1
)
= 'yes'
Pr
SCA
(
1
)
= 'yes'
=
Pr
SCA
(
r
)
= 'yes'
|
SCA
(
r
−
1
)
= 'yes'
·
Pr
SCA
(
r
−
1
)
= 'yes'
|
SCA
(
r
−
2
)
= 'yes'
·
…
·
Pr
SCA
(
2
)
= 'yes'
|
SCA
(
1
)
= 'yes'
Pr
SCA
(
1
)
= 'yes'
=
α
'
·
α
'
·
α
'
⋯
α
'
=
α
'
r
.
The rationale behind the above derivation is that
SCA
(
r
−
1
)
= 'yes'
already implies
SCA
(
r
−
2
)
= 'yes',
…,
SCA
(
1
)
= 'yes'
. In other words, we have
Pr
SCA
(
r
)
= 'yes'
=
Pr
SCA
(
r
)
= 'yes',SCA
(
r
−
1
)
= 'yes',
…,
SCA
(
1
)
= 'yes'
and
Pr
SCA
(
r
)
= 'yes'
|
SCA
(
r
−
1
)
= 'yes',SCA
(
r
−
2
)
= 'yes',
…,
SCA
(
1
)
= 'yes'
=
Pr
SCA
(
r
)
= 'yes'
|
SCA
(
r
−
1
)
= 'yes'
.
Let α be the targeted FPR, we have
α
=
∑
r
=
1
∞
Pr
SCA
(
r
)
= 'yes'
=
α
'
+
α
'
2
+
…
+
α
'
r
+
…
=
α
'
1
−
α
'
,
α
'
<
1
.
Accordingly, we have
α
'
=
α
/
1
+
α
.
Additional files
Competing interests
The authors declared that they have no competing interests.
Authors’ contributions
XY, GY and YW participated in the design of concepts and methods. XY and GY developed the permutation strategy and CNA simulation algorithm. XY implemented the C++ code. RRW implemented the R code of GISTIC. GY, XY and XH analyzed and evaluated the algorithm. XH and YW constructed and proved Theorem 1. YW, XY and GY drafted the manuscript. IMS and EPH interpreted the results on real cancer data. JZ, RC and EPH help edited the manuscript. YW, RC and ZZ conceived of the study, participated in its design and coordination, and helped edited the paper. All authors read and approved the final manuscript.