An optimization framework for unsupervised identification of rare copy number variation from SNP array data

A highly sensitive and configurable method for calling copy number variants from SNP array data is presented that can identify even rare CNVs


Background
Identifying DNA variants that contribute to disease is a central aim in human genetics research. Pinpointing these causal loci requires the ability to accurately assess DNA sequence variation on a genome-wide scale. In recent years, considerable progress has been made in identifying and cataloging single-nucleotide polymorphisms (SNPs) in many populations [1]. Commercial SNP microarray platforms can now genotype, with >99% accuracy, over one million SNPs in an individual in one assay [2,3].
The discovery of copy number variants (CNVs) as a significant source of variation has complicated the identification of genetic differences among humans. CNVs are defined as chromosomal segments at least 1,000 bases (1 kb) in length that vary in number of copies from human to human [4]. Since their discovery, several high-profile studies have been published associating copy number variation in the genome with a variety of common diseases. Recent examples include Alzheimer's disease [5], Crohn's disease [6], autism [7], and schizophrenia [8]. The significance of the gains (copy number greater than two) and losses (copy number less than two) that comprise these variants is increasingly evident, and cataloging them and assessing their frequencies has become an important goal.
SNP arrays contain hundreds of thousands of unique nucleotide probe sequences, each designed to hybridize to a target DNA sequence. When a DNA sample is properly prepared and applied to the array, specialized equipment can produce a measure of the intensity of hybridization between each probe and its target in the sample. The underlying principle is that the hybridization intensity depends upon the amount of target DNA in the sample, as well as the affinity between target and probe. Extensive processing and analysis of these raw intensity measures yield estimates of some characteristic of the target sequences in the sample -either target quantity [9,10], base composition [11,12], or both. In copy number inference, the objective is to identify chromosomal regions at which the number of copies per cell deviates from two; these include gains and losses.
There is now a large body of literature describing algorithms to infer copy number from SNP array data. All such algorithms address one or more of the three general steps: normalization, raw copy extraction, and CNV calling. Normalization is performed on the raw array intensity data in order to be able to compare these values fairly, thereby taking into account differences in overall array brightness and additional sources of nuisance variation. Raw copy number extraction entails converting the multiple measurements for each genomic site into a single raw measure of copy number. The word 'raw' here indicates that measurements from surrounding loci are not yet taken into account, and the measure is permitted to be non-integer. Since gains and losses occur in discrete segments often encompassing several such loci, true copy number is locally constant. Consequently, the final CNV calling step takes advantage of this fact, smoothing or segmenting the raw copy numbers into discrete segments of consistent copy number.
The Affymetrix SNP array was originally designed so that each SNP is interrogated by 24 to 40 unique probes. Of these, half are perfectly complementary to the sequence harboring the SNP site (perfect match probes), while half mismatch the sequence at the probe's middle nucleotide (mismatch probes). The mismatch probes were intended to capture background effects such as cross-hybridization. The perfect match/mismatch design was used for the 10,000, 100,000, and 500,000-SNP versions of the array. Most recently, Affymetrix has introduced the SNP Array 6.0, which interrogates nearly one million SNPs and differs fundamentally from previous versions. First, each SNP on the 6.0 array is interrogated only by six or eight perfect match probes -three or four replicates of the same probe sequence for each of the two alleles. Therefore, intensity data for each SNP consist of three or four repeated pairs of measurements. Second, the SNP probe sets are augmented with nearly one million CNV probes, which are meant to interrogate regions of the genome that do not harbor SNPs, but that may be polymorphic with regard to copy number. Each such CNV site is interrogated by only one probe.
For the Affymetrix platform, the community has largely settled upon quantile normalization [13] as a simple but effective normalization method. The next step, raw copy number extraction, typically entails fitting some model to raw probe intensity data [14][15][16][17]. Methods devoted to the final stepmaking CNV calls from raw copy number data -are numerous, and employ various strategies. Three commonly used strategies are hidden Markov models (HMMs) [17,18], circular binary segmentation [19,20], and adapted weight smooth-ing [21,22]. Although these methods appear to be quite different from one another in terms of the computational or statistical model they incorporate, at the core of each is an objective function whose optimum solution yields the method's copy number inference for a region. Each objective function is defined by the observed data (raw copy number) and is a function of inferred state (copy number call). The sequence of copy number calls (states) that optimizes the objective function gives the CNV call for each method.
In this paper, we present a general framework to call CNVs from raw copy number using optimization, based on an objective function that is composed of several explicitly formulated objective criteria. These criteria are carefully designed to quantify the desirability of a CNV assignment with respect to various biological insights and experimental considerations. Our general approach is to first apply a signal processing method to aggressively flag candidate gains and losses. The objective function is then optimized on each region and flanking sequence, yielding final CNV calls and boundaries. Note that the optimization process also filters out many candidate regions; that is, complete rejection of a candidate region is quite possible as it is part of the solution space for the corresponding optimization problem. This two-step procedure has the advantages of drastically reducing the computational time necessary to find the set of solutions, while identifying precise boundaries for each putative CNV. Indeed, for N markers and C CNV classes, the solution space of the optimal copy number assignment problem is of size O(C N ). Exhaustively searching for the optimal solution is quite infeasible unless N becomes very small. In our case, N ≈ 1.8 million, so we adapt a simulated annealing-based algorithm that efficiently searches the solution space at near-interactive rates.
We note here the distinction between CNVs and copy number polymorphisms (CNPs). CNPs are defined to be CNVs that are present and have identical boundaries (and are therefore likely identical-by-descent) in at least 1% of the human population [23]. Computationally, such higher-frequency polymorphisms present opportunities for detection that are not otherwise possible. A recent study [17] proposes separate methods to detect CNVs and CNPs, with the latter involving detecting correlations in raw copy numbers across samples. The current work is designed to address the problem of identifying rare and de novo CNVs, as it does not make use of multiple samples to convert raw copy number into CNV inferences.
A key feature of our method is that it is highly configurable, allowing researchers to define their own objective functions and tune parameters to emphasize the relative importance of different objective criteria. We demonstrate with a simple objective function involving a linear combination of variability, parsimony, and length, which performs surprisingly well. We evaluate the performance of our method on Affymetrix 6.0 array data from 270 HapMap individuals [1]. These sam-ples are increasingly well characterized with regard to CNVs and include 60 mother-father-child trios. Therefore, they serve as an excellent benchmark data set. We show via systematic in silico studies that the proposed method compares favorably with four methods that are currently publicly available. Furthermore, we experimentally validate, using laboratory techniques on genomic DNA, several CNVs newly discovered by our method. These results demonstrate the proposed method's potential to uncover human genetic variation that may be missed by other computational approaches.
The general framework described in this paper is implemented and freely available in a flexible, user-friendly R package ÇOKGEN*. ÇOKGEN works from the raw binary .CEL files produced by the Affymetrix protocol. It performs all of the steps in Figure 1, including quantile normalization, raw copy extraction, and CNV extraction (wherein the user may specify the desired objective function). Its graphical tools also allow the user to manually inspect the raw copy number data to gauge confidence in each putative aberration.

Results and discussion
We applied our algorithm to Affymetrix 6. Figure 2 graphically displays all CNVs identified by our method. As expected, many common CNVs are located near the centromeres and telomeres, which are known to harbor variably repetitive elements.
The distribution of the CNVs among different ethnicities in the population is presented in Table 1. It is well known that Asian and Caucasian populations are genetically less diverse than African populations due to population bottlenecks. This is reflected in Figure 3, which shows a shifted frequency distribution in the YRI CNVs relative to the CEU and JPT/CHB CNVs.

Trio discordance as a copy number variant detection assessment tool
Although CNVs can arise in a de novo manner, it is believed that at least 99% of all CNVs in an individual's genome are inherited [23]. The 60 mother-father-child trios in the Hap-Map data set therefore provide an opportunity to assess the accuracy of CNV detection algorithms by measuring the rate of Mendelian concordance. A CNV in a trio child is said to be Mendelian concordant if it appears in at least one of the parents. Unless the CNV is de novo, any discordance is either the result of a false positive call in the child or a false negative call in one of the parents (in rare cases, discordance could also result from a parent harboring a duplication and a deletion at the same locus but on different chromosomal homologs). Dis-cordance rate, while useful, is imperfect as an assessment measure. In particular, it is possible for a CNV identification algorithm to have artificially low discordance rates by calling each CNV in a large number of samples. Even if the samples in which a gain or loss is called are randomly selected, frequently called CNVs will have a lower discordance rate, simply by chance. Therefore, while comparing the performance of algorithms according to trio discordance rate, we also account for the number of frequently called CNVs, as discussed in the next subsection.
In the current study, to decide whether two CNVs (of the same type -loss or gain), c 1 and c 2 , from two different samples correspond to the same event, we use the concept of minimum reciprocal overlap. We first define o(c 1 , c 2 ) as the number of markers existing in both c 1 and c 2 and l(c) as the number of markers in a CNV c. Minimum reciprocal overlap (MRO(c 1 , c 2 )) of c 1 and c 2 is defined as: This measure provides a standard way of determining the similarity in the chromosomal location of two CNVs, regardless of the scale of the events. For our discordance and sensitivity analysis, we use the MRO measure with a threshold of 0.5 to decide whether two CNVs identified in two different individuals correspond to the same event. That is, at least half of c 1 must be overlapping with c 2 and vice versa for c 1 and c 2 to be considered as the same CNV in different samples.

Performance of ÇOKGEN in comparison to existing software
We compared the performance of our algorithm with that of four other software packages. The DNA-Chip Analyzer (dChip) [24] is a Windows software package for Affymetrix platform and high-level analysis of gene expression microarrays and SNP microarrays [14,25]. Birdseye [17] is a rare CNV identification tool based on HMMs, and is part of the Birdsuite platform [17]. QuantiSNP [26] is an analytical tool for the analysis of copy number variation using whole genome SNP genotyping data. It was originally developed for Illumina arrays, but version 1.1 of this software supports Affymetrix 6.0 data files with additional data conversion steps. PennCNV [27] is the last software tool that we use for CNV detection for our comparative analyses. Although it is also designed to handle signal intensity data from Illumina arrays, it currently supports Affymetrix.
Comprehensive experimental results show that ÇOKGEN outperforms all of these four CNV identification tools in terms of general trio discordance. Overall, ÇOKGEN has a 30.8% discordance rate whereas Birdseye, dChip, QuantiSNP and PennCNV demonstrate discordance rates of 42.6%, 94%, 74% and 32.9%, respectively, on the same array data. It is important to note that dChip was originally optimized for Overview of the proposed CNV detection algorithm Figure 1 Overview of the proposed CNV detection algorithm. ÇOKGEN first extracts the intensity values from the Affymetrix .CEL files. It then obtains the raw copy numbers for each marker using regression with the help of the Affymetrix software's SNP genotype calls. The edge detection determines the candidate loss/gain regions from smoothed copy number signal, which is obtained by low-pass filtering the raw copy numbers. We determine the final class assignments using objective function optimization. The function is optimized using an iterative simulated annealing procedure, with initialization provided by the edge detection.   detecting somatic copy number aberrations in cancer cells from earlier versions of the Affymetrix platform, and Quan-tiSNP is designed for data obtained from the Illumina platform. Therefore, Birdseye, PennCNV, and ÇOKGEN's superior performance compared to dChip and QuantiSNP on Affymetrix 6.0 data is not surprising. For this reason, we restrict our assessment to ÇOKGEN, Birdseye and PennCNV in the remainder of this section.
As discussed in the previous section, the expected discordance rate of any algorithm approaches zero as it calls the CNV in more samples. At the extreme, if the algorithm identifies a CNV in all samples, the discordance rate will be zero. Therefore, a more precise assessment of accuracy can be achieved by stratifying discordance rate by call frequency. For this purpose, in Figure 4, we first examine how the discordance rate behaves across call frequency strata for ÇOKGEN, PennCNV, and Birdseye. As a reference, we also display the expected discordance of randomly called CNVs in this figure. As expected, the performance of all algorithms improves when more frequent CNVs are considered. Although the performance of PennCNV is similar to that of ÇOKGEN, our algorithm does attain a modest improvement in concordance over PennCNV at all strata. It is also clear in Figure 4 that ÇOKGEN outperforms Birdseye significantly at all strata. Furthermore, ÇOK-GEN performs consistently better than random CNV assignment at all strata, which shows its superior performance is not an artifact of the frequency of the CNVs it calls.
Another feature of Figure 4 is Birdseye's sharper decline in discordance rate as the frequency threshold increases. This is likely due to its higher average call frequency compared to ÇOKGEN. Figure 5a shows the empirical density for sample frequency of concordant CNVs. We find that 34% of the concordant CNVs identified by Birdseye have frequency larger than 60, whereas only 16% of the concordant CNVs identified by our algorithm and 14% of the CNVs identified by PennCNV have frequency larger than 60. Concordant CNVs with sample frequency larger than 90 make up 3% of those called by our algorithm and 4% of those called by PennCNV compared to 22% for Birdseye. This clearly shows that ÇOKGEN does not achieve its high concordance rate by overcalling a CNV in multiple samples. Figure 5b displays the density distribution of discordant CNVs as a function sample frequency for all algorithms. It is clear from the figure that most of the discordant CNVs for Birdseye are rare, whereas more frequent CNVs called by our algorithm turn out to be discordant. These two observations clearly show that ÇOKGEN's performance depends less on the sample frequency and demonstrate its ability to accurately detect rare events.

Sensitivity comparison across methods
Trio discordance is a reasonable hybrid measure of sensitivity (recall) and specificity (precision), but these two measures cannot be easily decoupled based only on discordance rate. A recent study [28] assembled a 'stringent dataset' comprising CNVs identified by at least two independent algorithms. The dataset contains a total of 808 autosomal CNV regions reported by the study to be harbored in at least one of the 270 HapMap individuals. Another study [23] identified 1,292 autosomal CNP regions in 270 HapMap samples. We use these two as 'gold standard' data sets to evaluate the sensitiv-   Discordance rate as a function of call frequency strata Figure 4 Discordance rate as a function of call frequency strata. The figure shows how the discordance rates behave as a function of the sample frequency threshold. Note that discordance rate is plotted cumulatively -that is, the value on the y-axis is the average discordance rate for CNVs with frequencies, at most, the corresponding value on the x-axis. The discordance value at the sample frequency threshold value t is calculated by finding the discordance rate across all CNVs with frequency at most t. ity of our method. We refer to sensitivity based on the data presented in [28] as sensitivity-Pinto and sensitivity based on the CNP data set presented in [23] as sensitivity-McCarroll.
In terms of sensitivity-Pinto, we observe that ÇOKGEN detects 696 of 808 (approximately 86.1%) CNVs from the study presented in [28]. PennCNV obtains the best result by a narrow margin, by identifying 716 of 808 (approximately 88.6%) CNVs. Birdseye achieves an 84.7% success rate, slightly less than that of our method. In terms of sensitivity-McCarroll, ÇOKGEN and PennCNV detect 20.7% and 25.5%, respectively. Birdseye detects 68.2%, which is the best sensitivity rate among all the methods compared for this data set; however, as mentioned in [23], Birdeye is one of the methods used for identifying the CNPs in this dataset. For this reason, this result is not surprising. PennCNV is slightly more sensitive than our method on this dataset, though this seems to be at the cost of a modest increase in trio discordance rate, as shown above.

Run time performance
To analyze the run time performances of ÇOKGEN, Pen-nCNV, and Birdseye, we compare ÇOKGEN with PennCNV on a Windows system, and time both ÇOKGEN and Birdseye on a Linux system (Birdseye is not available in a Windows version). Performances are measured from the time at which the CEL file is taken as an input to the time at which the list of CNVs is output. On a Windows system that has an Intel Core 2 Quad CPU with a clock speed of 2.4 GHz and 4 gigabytes of memory, we observe that ÇOKGEN processes 22 chromosomes of a single HapMap sample in an average of 343 seconds compared to an average of 271 seconds for the PennCNV package.
The Linux experiments are done on a dual Intel Xeon 3 Ghz Centos 5 × 86 64-bit machine with 4 gigabytes of memory. Since Birdsuite is designed to be run as a pipeline of consecutive steps, we are unable to run only Birdseye in isolation. Thus, we report the run time for the whole package rather than single steps, which may admittedly inflate the time that Birdseye would take to run alone. In this experiment, ÇOK-GEN processes 22 chromosomes of a single sample in an average of 702 seconds compared to 2,232 seconds for the whole Birdsuite pipeline.
In addition to computational efficiency, these experiments also highlight the user-friendliness of our package. Indeed, ÇOKGEN is wholly contained in a single, simple (composed of three commands) R package, making it completely platformindependent and available to Windows, Mac, or Linux/UNIX users. In contrast to the competing software, ÇOKGEN does not require the installation of additional tools such as Active Perl [29] or Affymetrix Power Tools [30].

Experimental validation of copy number variants not previously reported
To gauge the ability of ÇOKGEN to uncover novel gains and losses, we compared the CNVs discovered by our method with those in version 6 (November 2008) of the Database of Genomic Variants [31]. We used multiplex ligation-dependent probe amplification (MLPA) [32] to verify some of the CNVs not reported in the Database of Genomic Variants but identified by ÇOKGEN (Table 2). In Figure 6, we also present the raw copy signal graphs generated by our software and the corresponding MLPA profiles for the first two CNVs given in Table 2.

The software package
Our software package, ÇOKGEN, is implemented in R and is able to output its results in two forms: tabular and graphical. The tabular output is a table of CNV entries with columns: sample ID, chromosome number, CNV start base position, CNV stop base position, and the CNV type. The graphical output allows the user to visualize the results of our CNV identification algorithm. The user can inspect the raw copy signal at any specified part of the genome along with the assigned, color-coded class values (examples are shown in Figures 6  and 7). Another aspect of the graphical output is the visuali- The frequency distribution of concordant and discordant CNVs for three calling algorithms

(b) (a)
zation of the signals of a family together, in which each member is represented by a different plotting symbol. This allows the user to see the CNV pattern for the whole family at the same locus of the genome and evaluate the algorithm's trio concordance visually. Besides its configurability in terms of tuning of parameters, ÇOKGEN also provides the user with the ability to specify their own objective criteria. With this functionality, users can construct their own objective functions that will best suit the characteristics and needs of their own experimental platform and application.

Conclusions
We present a method to detect germline CNVs from Affymetrix 6.0 SNP array data. Our approach, with its accompanying software, will be useful for researchers querying constitutional DNA for association of gains and losses with disease. Indeed, CNVs are emerging as important factors in a growing number of diseases, and the 6.0 array has the highest genome-wide resolution of current commercially available platforms. The current work shows that the problem of detecting CNVs from raw array data may be recast as an optimization problem with an explicit objective function. The objective function chosen here is quite simple and intuitive, but its effectiveness is clear. Our method is wholly contained in a freely available and flexible software package that efficiently processes raw probe-level .CEL files to produce lists of inferred gains and losses. The software allows the user to tune parameters for the desired specificity-sensitivity balance. With detailed experimental studies on the HapMap dataset, MLPA profiles and corresponding raw copy signals with class assignments for two CNVs not previously reported  50  60  70  80  90  100  110 120  130  140 150  50  60  70  80  90  100  110 120  130  140  we have demonstrated its sensitivity to detect both previously reported and novel CNVs, while keeping a low false positive rate, as demonstrated by high Mendelian consistency in trios.
The method described in this paper could also be adapted to other SNP arrays, including earlier versions of the Affymetrix platform, Illumina arrays, or array comparative genomic hybridization. Any platform that produces a measure of raw copy number at markers across the genome would be suitable. As SNP arrays continue to improve with regard to throughput and accuracy, our approach will be adaptable to handle the data as they become available. Raw copy numbers for sample NA12763 in a chromosome 12 region The optimization-based approach is the key to our method's flexibility. Although we have constructed our own default function to capture the criteria that we wish to emphasize, one may easily envision alternative criteria that other researchers would wish to incorporate. For example, since very long CNVs are quite rare in the human genome, researchers might wish to include a term in the objective function that takes into account the number of bases covered by a putative CNV region. Another possibility would be to incorporate allelic ratio intensity information at SNP markers, as is done in some HMM approaches [26,27]. We anticipate that users will design their own objective functions and apply them, using our software, to their own specific applications and data.
It should be emphasized that previously established approaches may actually also be considered special cases of functional optimization. For example, HMMs often used in the copy number setting [14,17,26] entail finding 'state paths' (marker-by-marker sequences of copy-number calls) that maximize a log-likelihood function. In HMM applications, however, the model parameters are often estimated simultaneously with the copy number states via a Viterbi algorithm [33], based on training samples. Precise parameter estimation relies on sufficient representation from each copy number state, which may be unrealistic for rare CNVs.
Another popular approach to inferring CNVs from raw copy number data is circular binary segmentation [19]. Rather than explicitly representing copy number state as a solution to an optimization problem, circular binary segmentation aims to find change points from one copy state to another. It does so by maximizing functions of marker indices. The optimum values of the function determine the boundaries of the CNV regions. A third example is the GLAD (Gain and Loss Analysis of DNA) algorithm [22], which has been adapted extensively using methods developed to analyze tumor DNA [15,34]. To find CNVs, GLAD explicitly models raw copy number as a function of position. The true underlying copy number is encoded in a position-dependent parameter. The CNV regions are inferred by maximizing a weighted likelihood function using an adaptive weights smoothing procedure [21]. Note that the objective functions in HMMs, binary segmentation and GLAD all make distributional assumptions about the raw copy number measurements. The function that we adopt in the current study makes no such assumptions, but could be modified to incorporate them. Furthermore, our CNV calling method is fully unsupervised in that it does not require any training samples in terms of known copy numbers. Lastly, rather than estimating and fixing parameters (thus fixing the performance of the algorithm), our method presents the opportunity to tune parameters, which makes it possible to adjust the performance of the algorithm to obtain the best results in a semi-automatic manner.
Three other studies have utilized various smoothing and edge detection algorithms: wavelet footprints [35], non-linear dif-fusion filtering [36], and kernel smoothing [37]. We also apply an edge detection scheme on low-pass filtered data to identify regions that potentially correspond to aberrations. Unlike other approaches, however, we apply edge detection rather aggressively to identify all candidate regions that may correspond to aberrations. This is because the raw copy number signal is extremely noisy due to the artifacts of microarray technology, as seen in Figure 7a. Furthermore, since the markers are distributed unevenly across the genome, the onedimensional signal represents a non-uniform sample of the actual copy number signal. Consequently, it is not straightforward to choose a smoothing and edge detection scheme that will be most appropriate for all experiments, samples, chromosomes, or even chromosomal segments. For example, in Figure 7b, the edge detection scheme identifies a single duplication as two separate duplications, since the markers at the middle of the region exhibit relatively low raw copy numbers, probably due to noise. This problem can be alleviated by smoothing the signal more aggressively to eliminate such artifacts, although this might result in falsely eliminating many aberrations that span relatively less numbers of markers. Motivated by these considerations, we use edge detection to identify all potential candidates and then use an optimization scheme with adjustable parameters to eliminate false calls among these candidates.
We also note that ÇOKGEN works on each sample individually and is therefore suited for rare CNV identification at the expense of losing some information to detect CNPs. The importance of rare CNVs is underscored by the recent deep sequencing of the entire genome of a single individual [38]. In that study, some 30% of the discovered CNVs had not been previously reported by any other study.
In addition to presenting a new software tool, the current work also casts Mendelian concordance, as an assessment tool, in a new light. While concordance rate is valuable as a metric to evaluate methods for calling germline variation, it is best viewed as a function of overall variant call rate. As we have shown, concordance rate can be artificially boosted simply by calling variants at a high rate. When evaluating the performance of future methods on family-based data sets, researchers may compare trio discordance results as a function of call frequency to the null expectation that we derive in the Materials and methods section.

Materials and methods
Our method takes as input raw .CEL files and produces a table of inferred genome-wide gains and losses. The software package, ÇOKGEN, provides a configurable platform for CNV identification, allowing users to: adjust the parameters of our default formulation to tune the behavior of the method to the target application (for example, aggressive versus conservative in calling CNVs); and to specify their own target objective functions. ÇOKGEN also produces 'zoomable' plots of raw copy number at the chromosome and sub-chromosome level for manual inspection of identified copy numbers. An overview of the methods implemented in ÇOKGEN is given in Figure 1. Details of each step are provided in this section.

Intensity extraction and normalization of raw data
The raw probe intensities for each array are encoded in the binary .CEL files output by the Affymetrix instrument, one file for each array. As a first step, we use the R package affxparser [39] to extract the intensities for each array locus from the .CEL files. Next, we quantile normalize [13] the intensities across all arrays in the experiment. This enables fair comparison of intensities, taking into account systematic non-biological differences such as overall array brightness.

Raw copy number for SNP markers
The genomic loci interrogated on the Affymetrix 6.0 array fall into two categories -SNP markers and copy number (CN) markers.  [23]. Therefore, if we define total probe intensity at marker i as: we can safely assume in general that most of the middle two quartiles, across all samples, of PI i are from individuals with two copies of the chromosomal segment that contains marker i. In other words, the individuals that fall into these quartiles for the corresponding marker are likely to carry diploid genotypes AA, AB, or BB. Consequently, we fit the model based on these samples' genotypes.
Given the 12 parameter estimates for a SNP marker i, we generate raw estimates of A and B copy numbers for all samples by re-applying the model to each sample's six probe intensities. That is, for a sample with probe intensity values A i1 , A i2 , A i3 , B i1 , B i2 , and B i3 , the raw A and B allele copy estimates are A i and B i where: Finally, using these estimates, we calculate the raw copy number R i at marker i as the distance from the origin in the (A, B) plane: 1  2  3  1  2  3 ,˘,˘,˘,˘,(

Raw copy numbers for CN markers
Approximately half of the marker loci represented on the 6.0 array do not correspond to SNPs, but rather CN markers. Since these markers are each measured by only one probe, they must be treated separately. As above, we consider the samples within the middle two quartiles of (normalized) total probe intensity for the marker to be representative of individuals with copy number two. Therefore, the scaling factor for CN marker i is the least-squares estimate of the parameter β i from the model: fit to the middle two quartiles of the normalized probe intensities PI i . Again, e i is the error term. The raw copy number for a sample with CN probe intensity PI i is then calculated as: Using these two separate procedures for SNP and CN markers yields raw copy numbers R i for all markers i from 1 to ~1.8 million, ordered along the genome according to hg18 (build 36 of the human genome) coordinates. All 270 HapMap samples are used to parameterize the regression model for raw copy number estimation of both SNP and CN markers. Figure  7a gives an example of raw copy numbers for a 394-marker region.

Algorithm for copy number variant detection
Key to our approach is the observation that CNV identification can be formulated explicitly as an optimization problem without any requirement of reference models or training data. Based on general knowledge of the microarray technology and basic biological insights on copy number variation, we specify various quantitative measures that gauge the suitability of copy number assignments based on observed array intensities. We then formulate an objective function that captures the trade-off between these measures, so that the minima of this function represent optimal CNV assignments. This function is characterized by user-defined parameters, allowing the user to tune the performance of algorithms based on the requirements of the specific application (for example, minimizing false positives due to the cost of experimental verification versus minimizing false negatives to capture existing variation comprehensively). We next introduce the objective criteria that are included in the default objective function implemented in ÇOKGEN and the motivation behind these criteria. Researchers may wish to design an objective function of their choice, and indeed our software takes the objective function as an argument precisely to accommodate this. We describe the function as applied to a chromosome with M markers since each chromosome is processed separately.

Variability in raw copy numbers within each copy class should be minimized
The R i for markers in each gain or loss region should be separable from normal regions. Therefore, CNV identification lends itself to a clustering-like problem -one of partitioning the R i s into three classes so as to minimize the internal variability of each class. For a given CNV assignment S, we define the set of markers assigned to class c on a chromosome with M markers as: and denotes the mean raw copy number for class c. Then, the total intra-class variability induced by this assignment is given by: Consequently, a desirable S is expected to minimize σ(S) (subject to other constraints). Note that this formulation does not make any assumption about the expected raw copy numbers of the markers and is therefore robust to any systematic bias that might be encountered in measurement and normalization of R i .

Parsimony principle: observed variability should be explained via minimum number of anomalies
In general, there are relatively few regions of gain or loss in an individual's genome, relative to normal regions. Therefore, the CNV calls should be as contiguous as possible. Motivated by this observation, we formulate the parsimony principle as a criterion that seeks to minimize the total number of copy number state changes induced by a CNV assignment on the chromosome. Formally, for given CNV assignment S, we define total cut as the number of pairs of adjacent markers that are assigned different copy numbers: Here I(.) denotes the indicator function (that is, it is equal to 1 if the statement being evaluated is true, and 0 otherwise).

Filtering out noise by eliminating smaller regions
Longer CNVs indicate higher confidence as it can be statistically argued that shorter sequences of markers with deviant raw copy numbers are more likely to be observed due to noise. Thus, we explicitly consider CNV length as an additional objective criterion. To do so, we first define a CNV region, r, as a maximal set of contiguous markers all assigned to the same copy number state in {C + , C -}, and ζ(S) denotes the set of all CNV regions. Furthermore, we denote the number of markers in the CNV region r by l(r). We then define as an objective criterion that penalizes shorter CNVs (e denotes the natural logarithmic base).

Filtering out noise by eliminating possible false positives
Candidate CNVs with a median raw copy number much larger or much smaller than two indicate higher confidence since a CNV region with median raw copy number close to two is less likely to be valid. For this reason, we require that the median raw copy number of a called loss be below a certain threshold (T loss ) and the median for a called gain be above a certain threshold (T gain ). We define ζ + (S) and ζ -(S) as the set of all CNV gain and loss regions, induced by assignment S, respectively. Furthermore, median(r) denotes the median raw copy number value of the markers in the region r. We now incorporate into the objective function to minimize the effects of the noisy signal. Here, T gain and T loss are user-defined parameters that basically define the upper and lower limits for the raw copy number of markers in the set Π(C 0 ) (that is, the set of markers assigned to the normal class). As T gain is increased and T loss is decreased, candidate regions are penalized more harshly. In our experiments, we use 2.25 and 1.75 for T gain and T loss , respectively, since these values provide reasonable performance.

Putting the pieces together: a single objective function for copy number variant identification
We use a linear combination of the criteria above as an objective function. Namely, we define the optimal copy number assignment as the mapping: such that the function is minimized at S = S*. Here the tunable coefficients k σ , k χ , k λ , k δ adjust the relative importance of the objective criteria with respect to each other. In our experiments, for k λ and k δ , we choose large values such as 10 5 and 10 6 , respectively, to prohibitively eliminate candidate regions that are likely to be false positive during the course of the algorithm (as opposed to filtering them out in a post-processing phase).
The parameters k σ and k χ are used to adjust the apparent trade-off between the 'parsimony' and the 'variability' components of the objective function. Variability favors the genetic diversity on the genome by permitting many CNVs. On the other hand, according to the parsimony criterion, the variability in the raw copy estimates of markers should be explained via as few CNVs as possible, hence minimizing the number of evolutionary events that have had to occur. Without loss of generality, we require that k σ +k χ = 1 to highlight the trade-off between these two criteria. To systematically evaluate the effect of these two parameters on performance and determine the best k σ and k χ values based on our benchmarking data, we conduct a series of computational experiments. The results of these experiments are presented in Figure 8. Here, sensitivity is a measure of the performance of the algorithm in capturing previously reported CNVs. As seen in the figure, sensitivity-Pinto rapidly improves as more weight is assigned to variability and nears saturation around k σ = 0.35 and k χ = 0.65. On the other hand, sensitivity-McCarroll keeps improving until it settles around k σ = 0.6 and k χ = 0.4. In contrast, discordance rapidly declines as we increase the contribution of variability to the objective function, achieves a minimum around k σ = 0.35 and k χ = 0.65, and grows until it settles around 0.8 for k σ = 0.8 and k χ = 0.2. As The trio discordance and sensitivity as a function of k σ Figure 8 The trio discordance and sensitivity as a function of k σ . The figure shows how trio discordance and sensitivity are affected as we alter the relative weights of variability and parsimony in the objective function. Based on these observations, we set k σ = 0.35 and k χ = 0.65 as our defaults.

Two phase copy number variant identification
Since the solution space of the optimal copy number assignment problem is exponential, we require a good initial solution and a heuristic algorithm that iteratively improves the solution. For this purpose, we use a two-phase algorithm: we first determine a set of candidate gain and deletion regions via a filtering and aggressive edge detection procedure that we consider as an initial CNV assignment, S (0) ; and then we employ an iterative improvement based algorithm to adjust the boundaries of duplications and deletions accurately, and eliminate false positives.
In order to identify the boundaries of CNV regions, it is necessary to smooth the raw copy number signal since it is highly noisy. We use a simple discrete low-pass filter with filter kernel [1/3; 1/3; 1/3], that is, the first filtered copy number estimate is given by: Applying the filter for a second time, we obtain: Consequently, introducing an adjustable repetition parameter W, we obtain as a smooth version of the copy number intensity for a user defined value of W. Here, larger W provides smoother signals, thereby eliminating false positives, at the cost of missing true CNVs that span a smaller number of markers. For the ÇOKGEN's default value, we chose W = 20, for which we obtain reasonable results. Figure  7b demonstrates how the raw copy number R i in Figure 7a is converted into a smooth signal using the low pass filter.

Identification of candidate copy number variation regions via edge detection
Based on the observation that gains and losses manifest themselves as (respectively up or down) concavities in raw copy number of the low-pass filtered data, an edge detection scheme, which we describe below, is a useful tool for the identification of initial CNV assignment S (0) . Thus, after low-pass filtering, we apply our edge detection algorithm on the smoothed signal, first identifying high gradient markers that may correspond to transitions between regions with different copy numbers. For this purpose, we interpolate the discrete signal to obtain a real-valued function on the continuous interval . This task is performed using the built-in splinefun function of the R language, which performs cubic spline interpolation of given data points. Next, we generate two sets of high-gradient markers, denoted D max and D min , for which the function attains maximum increase and maximum decrease, respectively. Specifically, we define: where denotes the derivative of at marker i. These markers are the approximate inflection points of the signal . Now let Q ij denote the indices corresponding to the set of contiguous markers on the genome starting from marker i and ending at marker j, where i ≤ j. Given the user defined thresholding parameter T gain (see above), we designate Q ij as a candidate gain region (that is, ∀ k ∈ Q ij , S (0) (k) = C + ) if it satisfies the following conditions: i ∈ D max and j ∈ D min ; there exists at least one marker p, i ≤ p ≤ j, such that ; max(Q ij ∩ D max ) < min(Q ij ∩ D min ); and Q ij is a maximal set of contiguous markers satisfying the first three conditions.
The first condition ensures that the region starts with a marker with locally maximal positive gradient and ends with a marker with locally maximal negative gradient in terms of the raw copy number values. The second condition guarantees that the region contains markers with copy number estimates that might indeed correspond to a gain. The third condition specifies that the region does not contain any interior concavities, that is, all maximum positive gradient markers in Q ij appear before any maximum negative gradient marker in the region. Finally, the fourth condition ensures that Q ij can be enlarged at neither the right nor the left borders. Examples of regions that violate each of the conditions are shown in Figure 9. The designation of Q ij as a candidate loss region is done in a completely analogous manner.
All markers m that are not included in a candidate loss or gain region are preliminarily designated as 'normal', that is, S (0) (m) is set to C 0 . As a special case, if a candidate gain/loss region identified by edge detection is very close to another candidate region of its type, then we merge these two candi- date regions into a single region, since they are likely to correspond to the same aberration.
This procedure gives us an initial CNV identification assignment S (0) . This solution is quite aggressive in the sense that many truly normal (copy number two) markers are likely to be placed in the gain or loss classes. To eliminate these false positives and obtain S*, we use an optimization-based algorithm to tune the boundaries of candidate gain and deletion regions as discussed in the next section.

Fine tuning of the region boundaries using optimization with simulated annealing
This phase of the algorithm begins with initial class assignments, S (0) , and iteratively improves them with regard to the value of the objective function f by making moves in a way to quickly reach an optimum and avoid being trapped into undesirable local optima. For a given copy number assignment S, we define a 'move' as the extension or contraction of a CNV region's boundaries by changing the copy number states assigned to a contiguous group of markers (either inside or outside the region) bordering the region. In short, at each iteration of the algorithm, a random number of contiguous markers is selected from the right or left boundary of a candidate region Q ij ∈ ζ(S) and the corresponding move is defined as the assignment of these markers to either the class of neighboring markers (if the selected markers belong to Q ij ) or to Q ij 's class (if the selected markers are outside of Q ij ). The concept of a move is illustrated in Figure 10. As seen in the figure, we restrict possible moves to those that can enlarge or shrink a candidate aberrant region, but can never create a candidate region from scratch or divide a candidate region into two candidate regions. The size of the valid moves set that can shrink a candidate region Q ij of size n is 2n -1. This set contains all moves that change the class value of contiguous markers in either the left or right boundary of Q ij to the class of neighboring markers outside Q ij . The size of the valid moves that can enlarge a candidate region Q ij is limited to 2ψ (that is, at most, ψ markers from the left and ψ from the right are converted to the class of region Q ij ) where ψ is a userdefined parameter that determines how aggressively a candidate region is to be enlarged. In our experiments, we set ψ = 5, which limits the permissible expansion of a CNV region. We set such a threshold since we want our algorithm to Illustration of criteria for selection of candidate CNV regions Figure 9 Illustration of criteria for selection of candidate CNV regions. In this figure, the red lines tangent to the signal curve indicate the marker points that are elements of the sets D max and D min . (a) Q ij violates the first condition, which requires that the start and end markers be elements of D max and D min , respectively. (b) Q ij violates the second condition, which dictates that at least one marker must exceed the T gain threshold value. (c) Q ij does not satisfy the third condition since three markers of D min appear before two markers of D max in region Q ij . (d) Q ij violates the fourth condition, which requires that Q ij is a maximal set satisfying all the first three conditions.
expand the candidate region gradually. Thus, the total number of valid moves at each stage of our algorithm for a candidate region of size n is 2n + 9.
We quantify the quality of such a potential move in terms of the difference between the value of the objective function before and after the move, commonly referred to as the 'gain' of a move. In the context of the optimal copy number assignment problem, the gain associated with move ν is defined as: where S (t+1) denotes the copy number assignment if move ν is made and S (t) is the current copy number assignment. We use a stochastic algorithm based on simulated annealing [40] to determine the move. Simulated annealing is an iterative improvement heuristic that proceeds by repeated moves to improve the quality of the solution. Key to its efficiency is the stochastic nature of the selection of moves. At each step, the algorithm first randomly chooses a candidate gain or loss region, Q ij , from the set ζ(S) and then chooses a move v from the set of all moves that are validly defined on Q ij . If the gain γ(v) associated with the candidate move is positive, then the move is made. If the gain is not positive, the move is still made with a certain probability, which is proportional to the gain and declines as a function of time in the course of the algorithm. Therefore, simulated annealing starts its course with aggressive moves to jump out of undesirable local optima, and becomes more conservative as the algorithm proceeds, smoothly converging to a locally optimum solution. The procedure is repeated until either there is no valid positive gain move left to be done on the current solution or a user-defined number of negative gain moves, τ, have already been done consecutively (we use τ = 5 as our default). The mapping obtained at the end of the procedure is reported as S*.
We note that our algorithm allows us to consider the candidate regions in ζ(S) independently (as opposed to the entire chromosome) because the candidate regions with potential aberrations are sparse, and we therefore work on local subproblems associated with each candidate region separately. This results in significant improvements in computational efficiency. Since the distribution of raw copy numbers in the neighborhood of Q ij ∈ ζ(S) provides a good sampling of raw copy numbers in the entire chromosome, the quality of the solution to this local problem does not deviate significantly from the global problem. Indeed, in our experimental evaluations, we observe that there is no significant difference between solving the class assignment globally (applying the above algorithm to a whole chromosome) or locally (as we describe above) in terms of their specificity and sensitivity in predicting copy number variations.

Data
For the application of our method, we used Affymetrix 6.0 array data from a total of 270 HapMap individuals. In the data set, there are 30 mother-father-child trios from the Yoruba people of Ibadan, Nigeria, 45 unrelated individuals from Tokyo, Japan, 45 unrelated individuals from Beijing, China and 30 Caucasian trios that were collected in 1980 from US residents with northern and western European ancestry by the Centre d'Etude du Polymorphisme Humain (CEPH).

Multiplex ligation-dependent probe amplification method
Each MLPA probe was designed synthetically to match sequences within the region of interest avoiding all SNPs in the area. Control probes were used from previously published work [41]. Oligonucleotides were synthesized by IDT (Coralville, IA, USA) with 5'-phosphorylation of each downstream Illustration of moves in the proposed iterative-improvement based optimization algorithm Figure 10 Illustration of moves in the proposed iterative-improvement based optimization algorithm. The CNV assignment for a hypothetical region after t iterations is shown in the initial figure. Suppose Q ij is selected for the (t+1) st iteration. (a) All the markers in Q ij are assigned to C 0 , which completely eliminates Q ij as a gain region. This is a valid move. (b) Some markers that initially have a C 0 class on the right border of Q ij are assigned to the C + class which merges them with Q ij by a valid move. (c) Some markers that initially exist in Q ij are assigned to the C 0 class, which contracts Q ij . This also represents a valid move. (d) An invalid move that divides the Q ij into two sub-gain regions. (e) Another invalid move that introduces a completely new gain region that is not identified in the previous solution.

(e)
probe and tagged with common PCR primer sequences [32]. Probes were hybridized with 100 ng of DNA sample using MLPA reagents (part number EK1, MRC-Holland BV, Amsterdam, The Netherlands) in accordance with the recommended protocol. Samples were diluted 20-fold and analyzed on a 3130xl Genetic Analyzer from Applied Biosystems (Foster City, CA, USA) with GeneMapper software. Control DNA used were male and female genomics DNA pools (Promega, Madison, WI, USA). Peak height ratios were normalized to the mean of the entire data set, with subsequent elimination of outlier samples from the calculation of the mean.

Other methods
For analysis using dChip [24], we downloaded the version with a build date of 21 August 2008. We used the HMM as the Inferred copy method option with 50% of the samples trimmed. For Birdseye [17], we used version 1.5.1 of the Birdsuite package, which can be downloaded from [42]. The default parameters for that package were used. The latest version of PennCNV software, available as of 18 November 2008, was downloaded from [43] for analysis using PennCNV. We followed the steps described at [44] for the PennCNV-Affy protocol and used the default parameters for analysis. For QuantiSNP, we downloaded version 1.1 from [45], followed the steps described at QuantiSNP in the Affymetrix tutorial document located at [46] and used the default parameters.
We also note that we combined copy number 0 and 1 into one category -loss -and copy number greater than 2 into one category -gain -for the results obtained by all packages, in order to compare their results with ÇOKGEN's results fairly.

Computation of expected discordance rate
Suppose that CNV calls are random in Φ parent-child trios. When randomly assigning a CNV to n of the 3Φ individuals, the expected discordance, denoted by ED(n), can be calculated as: where ED(n | k) denotes the expected discordance rate when k of the n CNVs are assigned to children and P(k) denotes the probability of assigning the CNV to k children. We calculate P(k) as: ED(n | k) is equivalent to the probability of having a child assigned a CNV being discordant given that k children have the CNV. This probability can be calculated by dividing the number of ways to assign this CNV to parents other than the discordant child's parents by the number of all possible ways to assign this CNV to parents. Thus, ED(n | k) is: Consequently: .