The CNVrd2 package: measurement of copy number at complex loci using high-throughput sequencing data

Recent advances in high-throughout sequencing technologies have made it possible to accurately assign copy number (CN) at CN variable loci. However, current analytic methods often perform poorly in regions in which complex CN variation is observed. Here we report the development of a read depth-based approach, CNVrd2, for investigation of CN variation using high-throughput sequencing data. This methodology was developed using data from the 1000 Genomes Project from the CCL3L1 locus, and tested using data from the DEFB103A locus. In both cases, samples were selected for which paralog ratio test data were also available for comparison. The CNVrd2 method first uses observed read-count ratios to refine segmentation results in one population. Then a linear regression model is applied to adjust the results across multiple populations, in combination with a Bayesian normal mixture model to cluster segmentation scores into groups for individual CN counts. The performance of CNVrd2 was compared to that of two other read depth-based methods (CNVnator, cn.mops) at the CCL3L1 and DEFB103A loci. The highest concordance with the paralog ratio test method was observed for CNVrd2 (77.8/90.4% for CNVrd2, 36.7/4.8% for cn.mops and 7.2/1% for CNVnator at CCL3L1 and DEF103A). CNVrd2 is available as an R package as part of the Bioconductor project: http://www.bioconductor.org/packages/release/bioc/html/CNVrd2.html.


INTRODUCTION
Copy number variation (CNV) encompassing genes is a common phenomenon in the human genome, and has been shown to be associated with variation in phenotype Freeman et al., 2006;McKinney et al., 2008McKinney et al., , 2010Bentley et al., 2010). Accurate CN assignment is required for studies associating CNV loci with phenotype. However, accurately measuring gene CN by direct molecular methods over large numbers of samples is challenging, and is often complicated by the existence of paralogous gene pairs. As a consequence, phenotypic association data, largely obtained by quantitative polymerase chain reaction (Q-PCR) based methods, should be regarded with caution (He et al., 2009;Marques et al., 2010;Carpenter et al., 2011;McKinney and Merriman, 2012). In response to the increasing utility of high-throughput sequencing (HTS) technologies, we previously developed the read depth based method CNVrd (Nguyen et al., 2013). We demonstrated its ability to accurately assign CN using genome-wide HTS data duplications and deletions at the FCGR locus, where CN ranges from zero to four. To date, however, the utility of this approach for genotyping loci at which a greater range of CN exists is untested.
One complex CN variable locus is CCL3L1. Copy number variation at CCL3L1 has been associated with susceptibility to HIV infection (Liu et al., 2010), autoimmune disease (Burns et al., 2005;Mamtani et al., 2008;McKinney et al., 2008) and asthma (Lee et al., 2011). The median CN of CCL3L1 is 2 in European populations and >2 in other populations . Evaluation of the role of CCL3L1 CNV in common disease has been hampered by robustness of methodology, particularly that based on Q-PCR (He et al., 2009;Carpenter et al., 2011). Similarly, CNV within the beta-defensin locus on chromosome 8, which includes the DEFB4 and DEFB103A genes that vary in CN en bloc (a range of two to nine copies) (Groth et al., 2008), has been associated with various infectious and inflammatory phenotypes (Bentley et al., 2010;Mehlotra et al., 2012;Zhou et al., 2012). As is the case for CCL3L1, however, interpretation of these data, and further progress in studying the role of beta-defensin CN in medical phenotypes, is hampered by the difficulty in accurately measuring CN in large numbers of samples (Aldhous et al., 2010). Paralog ratio test (PRT) copy number data are available at both CCL3L1 and the beta-defensin locus (Armour et al., 2007;Walker et al., 2009), and thus provide an important resource for developing and validating new HTS read depth-based copy number assignment approaches. The PRT uses multiple probes to compare CN to an invariant paralog and is currently regarded as the gold-standard method for measuring CN at complex loci (McKinney and Merriman, 2012).
Here we describe the development of CNVrd2 as an extension of the CNVrd methodology, and demonstrate the application of this new approach by assigning copy number at the development locus CCL3L1, and test locus DEFB103A, both of which have previously had CN measured using the PRT.

THE WORKFLOW OF THE PIPELINE
The entire pipeline is represented in Figure 1. The pipeline is used to measure CN for a specific gene or locus. Four main steps are described below.

PIPELINE WORKFLOW: IDENTIFICATION OF APPROXIMATE BOUNDARIES OF MEASURED REGIONS
Samples were segmented as previously described for the CNVrd procedure (Nguyen et al., 2013). The results of the segmentation process were used to obtain approximate boundaries of copynumber variable regions encompassing the gene being measured. CNVrd2 used different percentiles and standard deviations of segmentation results to visually obtain information of CNVs at loci. Then, standard deviations were used to obtain approximate boundaries as follows (Step C, Figure 1).
To identify approximate boundaries of putative polymorphic regions (Step C, Figure 1), the coordinates of regions having similar segmentation results for each sample were obtained from the segmentation process applied to the large 1-2 Mb regions across all samples. Next, all these coordinates were combined to generate sub-regions, and segmentation results generated by the original segmentation analysis for each sample were obtained. Finally, the standard deviations and various percentiles of segmentation results of each sub-region across the 2535 samples were calculated and plotted to visually confirm CNV. Windows exhibiting large standard deviations were used to identify approximate positions of putative regions being assigned CN.

PIPELINE WORKFLOW: OBTAINING SEGMENTATION SCORES FOR INDIVIDUAL SAMPLES AT THE MEASURED LOCUS
This step was similar to our previous pipeline (CNVrd) but with some modifications (Step D, Figure 1) added to obtain more reliable segmentation scores. If the entirety of the gene of interest was placed in a single segment then the segmentation score for the gene was the segmentation score for that region. If the gene was split into multiple segments, each having the same sign then the segmentation score for the gene was calculated as the lengthweighted average across the segments (NA18612 for example in Step D, Figure 1). Otherwise, if a segment at the boundary of the gene had a different sign from other sub-regions of the gene (HG00254 for example) and the boundary segment's length was less than a specific threshold (here the threshold was set to half the window size) then the z-score for the boundary segment was calculated, where the z-score is defined as the value of the standardized observed read-count ratio (Nguyen et al., 2013). If the z-score had the same sign as the segmentation scores from other sub-regions of the gene, then the segmentation score for the gene was again calculated as the length-weighted average across the other segments. For all other situations the gene of interest was assigned a segmentation score of zero, reflecting the population average at that locus.

PIPELINE WORKFLOW: ADJUSTMENT OF THE SEGMENTATION PROCESS FOR MULTIPLE POPULATIONS
The pipeline was applied to each of the five major populations (European, East Asian, West African, South Asian, and Americas) to obtain segmentation scores (SSP: the segmentation scores for a single major population). This identified the scores in each of the large populations, but did not allow comparisons between the populations. Therefore, we ran the segmentation process for all samples and obtained segmentation scores for samples in all populations [SSP (P) : the segmentation scores for pooled populations]. When we pooled samples and ran the pipeline, if the gene being analyzed was not segmented into a single region for a particular sample, then that sample was assigned a segmentation score of zero. For example, the median CCL3L1 CN for Europeanancestry sample sets was less than other populations, therefore the majority of segmentation results and z-scores (where neededsee above) were less than zero. If the CCL3L1 gene region of a European-ancestry sample was segmented into two sub-regions, both having negative signs in the all-population analysis, then this could either be recapitulated in the single-population analysis (implying they were within the same CN region) or they could be sub-regions that could span the population average. Thus, they could have one negative and one positive sign (implying they were within two distinct CN regions). In this case, however, the pipeline only recognized that sub-regions shared the same sign, and the two regions would therefore be considered as one region and would be merged. To improve this, for each of the five major populations, we fitted a linear regression line with SSP and SSP (P) as independent and dependent variables respectively (data for CCL3L1 are shown in Figure 2). This step aimed to reduce errors when performing segmentation within single major populations or pooled populations. The fitted mean values of the five single linear regression models were used as the final per sample segmentation scores for these five populations. Segmentation scores were transformed and standardized, and a normal mixture model was used to cluster the scores into CN genotype groups.

PIPELINE WORKFLOW: ASSIGNING COPY-NUMBER COUNTS
A normal mixture model was used to cluster final segmentation scores into different groups. If clear clusters were observed then the Expectation Maximization (EM) algorithm (Dempster et al., 1977) was used in the clustering process. If data were complicated (e.g., loci with high copy-number variants such as CCL3L1 and DEFB103A) then a Bayesian clustering approach was used (Step F, Figure 1; details of the analysis process for CCL3L1 and DEFB103A are described below).

IDENTIFICATION OF WINDOW SIZE PARAMETERS FOR CNVrd2 USING TEST LOCI
One of the new innovations in CNVrd2 was the merging of segmentation results inside regions being measured if these regions www.frontiersin.org August 2014 | Volume 5 | Article 248 | 5 FIGURE 2 | Linear relationship between the segmentation scores called for single populations and for all populations at CCL3L1. SSP, the segmentation scores of a single large population; SSP( P ), the segmentation score of pooled populations.
were divided into different parts (Step D, Figure 1). We merged these parts if their signal values were above (below) a threshold (namely "testThreshold2Merge" in the CNVrd2 package) for duplicated events (deleted events). Another parameter was the length of a boundary region, the default for which was set to "half window size." Real data sets were used for investigation of the parameters of CNVrd2. In order to make reliable data sets we intersected the CNV results of two sets of data (Conrad et al., 2010;Campbell et al., 2011). From these intersecting data we chose eight loci whose copy-number assignments were ≥45% identical and whose lengths were in different ranges (∼3, ∼5, ∼8, ∼20, ∼24, ∼45, ∼60, ∼100 kb). The coordinates and CN status of the 8 loci are described in Table 3. We only retained samples having identical results between the methods. We downloaded BAM files of the 1 Mb regions around these loci and let CNVrd2 automatically assign CN for the loci (CN groups of these loci were known in advance from the micro-array based results, Table 3).

COMPARISON WITH CNVrd
We compared the new pipeline with our previous pipeline at the eight loci mentioned above (Table 3) and the eleven loci previously measured in Nguyen et al. (2013). For the eight loci, both pipelines were run and the percent concordances were calculated.

APPLICATION OF CNVrd2 TO TWO COMPLEX LOCI: CCL3L1 AND DEFB103A
Windows of 500 and 1000 bp were used to analyze CCL3L1 and DEFB103A, respectively. Steps described above were used to obtain putative boundaries of CCL3L1/DEFB103A genecontaining regions and final segmentation scores for the two loci. For each locus, one population having clear clusters of segmentation scores was used to obtain prior information for all populations: the European and South Asian population for CCL3L1 and DEFB103A, respectively. This approach was implemented using the EM algorithm (Dempster et al., 1977). The results for the European and South Asian populations (means  Conrad et al. and Campbell et al. (Conrad et al., 2010;Campbell et al., 2011) and the 1000 Genomes Project. b In parentheses are the number and percentage of samples with the specified CN, as obtained from Conrad et al. (2010) and Campbell et al. (2011). and variances of groups and the distances between the means of the groups) were used to provide prior information for a Bayesian normal mixture model that was applied to the other populations. This was done to improve model stability and fit for populations exhibiting high levels of CN polymorphism. The Markov chain Monte Carlo sampling method implemented in the rjags package (Plummer, 2011) was used to obtain posterior estimates for each normal mixture model parameter-this allows control of parameters using prior information. An adaptive phase of 100 iterations was run, followed by a burn-in period of 1000 iterations. Next, we ran 20,000 iterations and calculated the means, standard deviations, and proportions of the mixture components from 20% of the iterations. Convergence was assessed by using Heidelberger and Welch's stationarity and halfwidth test implemented in the coda package (Plummer et al., 2006).

COMPARISON WITH OTHER METHODS: CCL3L1
Copy number assignments at CCL3L1 by cn.Mops (Klambauer et al., 2012), CNVnator (Abyzov et al., 2011) and Sudmant et al. (2010) were compared to assignments of CNVrd2 on the 180 samples measured by the modified PRT. To run CNVnator, we downloaded all of the chromosome 17 data for the 180 samples. CNVnator was applied as previously described (Nguyen et al., 2013).  Figure 4. At CCL3L1 a total of 395,078,047 reads across all samples were aligned to a 2 Mb region around the gene (Chr17:33670000-35670000). These reads had lengths ranging from 36 to 160 bp (median of 91.5 bp and the highest frequency, 44.9%, was 100 bp), and mapping qualities from 0 to 70 (median of 33.5 and the highest frequency, 73.4%, was 60). The majority of reads (73.5% and 72.5%) aligning to the CCL3L1-containing region (chr17:34617501-3465201) and CCL3L1 gene (Chr17:34623842-34625730) had a mapping quality of 0, presumably reflecting multiple alignments to the paralogs within the locus.

Alignment results are presented in
A total of 386,195,488 reads across all samples were aligned to the 2 Mb region on Chr 8 (Chr8:6500000-8500000) around

FIGURE 4 | Read lengths and mapping qualities (top), mapping qualities (middle) and average read depth (bottom).
Data for the 2 MB CCL3L1 region are on the left and the 2 Mb DEFB103A region on the right.
the defensin genes (Figure 4). These reads had lengths ranging from 36 to 160 bp (median of 91.5 bp and the highest frequency, 45%, was 100 bp), and mapping qualities from 0 to 70 (median of 33.5 and the highest frequency, 40.2%, was 60). The majority of reads (64.1 and 99%, respectively) aligning to the DEFB103A-containing gene (chr8:7641001-7742001) and DEFB103A region (Chr8:7738726-7740105) had a mapping quality of 0. The average coverage of 2535 samples on the 2 Mb region was between 2.8 and 30.9× with a median of 6.9× (Figure 4).

CNVrd2: INFERRING RELIABLE VALUES FOR PARAMETERS
We used eight loci derived from the data sets of two microarray approaches (Conrad et al., 2010;Campbell et al., 2011) to obtain parameters of CNVrd2. CNVrd2 had good performance when the window size was 200 or 500 bp and a testThreshold2Merge value of 0.35 was used (Figure 5). A small window size (e.g., 100 bp) was not stable at all loci, possibly because there were multiple low-coverage samples in the 1000 Genomes Project data (e.g., 464/2535 = 18.2% samples having coverage in the CCL3L1 genecontaining 2 Mb region < 5×, Figure 4). This could lead to some windows having very low read counts, and thus being segmented incorrectly. For example, some windows in normal regions (i.e., CN = 2) would be segmented into deleted regions or some windows in duplicated regions would be segmented into normal regions. Similarly, large windows (e.g., 5000 bp) also generated unreliable results (Figure 5).

CNVrd2: COMPARISON WITH CNVrd AT EIGHT TEST LOCI
Our previous work was focused on loci having low copy-number ranges. Even though this present study is concentrated on loci with CN range from zero to >5 CN, we also compared with CNVrd (Nguyen et al., 2013) to assess the performance of the modified pipeline, especially in the steps requiring obtaining and adjusting segmentation scores for each sample. Loci used in this comparison had copy number ranging from 0 to 5, and their boundaries were known (Table 3 and Nguyen et al., 2013). Therefore, a simple normal mixture model was used in CNVrd2 to cluster segmentation scores into different groups. CNVrd2 had more stable results than those of CNVrd at almost all window sizes (Figures 5, 6). Merging sub-regions and adjusting boundaries (Step D and E, Figure 1) produced more reliable results for CNVrd2 at two complex regions (chr17:34436344-34481815 and chr1:110222301-110242933). At these two regions, CNVrd showed very low concordance with microarray-based approaches (Figure 6).

CNVrd2: COMPARISON WITH THE ELEVEN LOCI USED IN CNVrd
We also applied CNVrd2 to eleven CNV regions (including FCGR3A/3B) to which CN was assigned using CNVrd in our previous work (Nguyen et al., 2013

APPLICATION OF CNVrd2 TO TWO COMPLEX LOCI; CCL3L1 AND DEFB103A
Information on read length, mapping quality, read depth and read count ratios is presented in Subjects and Methods ("HTS and microarray data used") and in Figure 4. Based on the results of the eight test loci (above, Figures 5, 6), the BAM files were processed using 500 bp windows for CCL3L1 and 1000 bp windows for DEFB103A, and a testThreshold2Merge = 0.35. 1000 bp windows were chosen for DEFB103A because there was higher correlation between segmentation scores of DEFB103A and DEFB103B than for 500 bp windows (r = 0.80 and 0.78, respectively). Owing to the small length of each gene (1888 bp for CCL3L1 and 1379 bp for DEFB103A) larger window sizes were not investigated as this could result in not detecting breakpoints nearby or inside the genes.

POLYMORPHISM OF THE CCL3L1 LOCUS
Using a 2 Mb region it took approximately 3 h to assign copy number at CCL3L1 on 2535 samples using a 4-core computer with 8 Gb memory. Standard deviations and different percentiles of the segmentation results of sub-regions on the 2-Mb region were calculated to visually detect this region. Using a standard deviation threshold of 0.5, we identified CCL3L1 within a large polymorphic CN variable region of approximately 329 kb (chr17:34436001-34815000, which included a 50 kb gap) (Figure 7). This region included small CNV blocks of clearly defined increased standard deviation (SD) which suggested that there would be multiple recombination events inside the CNV region. CCL3 was outside this region, CCL4 was at the boundary, while CCL3L1 was in a second block with CCL3L3 and CCL4L1. The CCL3L1 gene was located in a sub-region which had SDs fluctuating around 0.84 and 0.85. CCL3L3 and CCL4L1 were also in a sub-region with SDs fluctuating around 0.83 and 0.86, but these sub-regions were separated from the CCL3L1 genecontaining sub-regions by a decrease (0.8) in SDs between two blocks (Figure 7). We merged these sub-regions and the boundaries of the CCL3L1-containing region were chr17:34617501-34652500. Using CNVrd2, we obtained segmentation scores for both the CCL3L1 gene-containing region and the larger CCL3L1 region. The segmentation scores of the two regions were strongly correlated (r = 0.98, Figure 8).

ASSIGNING CCL3L1 COPY NUMBER
The European CCL3L1 segmentation results from the smaller 35kb gene-containing region were assigned by CNVrd2 to five groups (Figure 9). All clusters were clearly delineated, although the highest group exhibited some scatter. The results obtained from the EM algorithm applied to the European samples were used as prior information to obtain CCL3L1 CN estimates for the 2535 samples from the 1000 Genomes Project. The final segmentation scores across all samples ranged from −1.66 to 4.53 (Figure 9), with only one score being larger than 3.42. This value was considered to be an outlier, and the associated sample (HG00620 in the CHS population) was assigned to the largest copy number group with a probability of 1. Thus, the segmentation score range used to determine the number of groups was between −1.66 and 3.42. The distances between the groups in the European population ranged between 0.39 and 0.62. As a result, 10 groups were chosen to encompass the range of CN values for the full collection of samples. The Markov Chain Monte Carlo chain was convergent after 10,000 iterations and the segmentation scores were clustered into 10 groups corresponding to CN from 0 to ≥9 (Figure 9).
The CN assignments for European-ancestry populations ranged between 0 and 5, with the highest frequency (43.6%, 220/505) being 2 copies ( Table 1). Similar to the Europeanancestry populations, CN = 2 was the most common (42.1%; 208/494 samples) in the South Asian sample set. The East Asian sample set had CN between 0 and ≥9, with the most common being CN = 3 (36.7%, 189/515). CN = 3 was also the most frequent in the Americas sample set, with the highest CN = 8. The African sample sets had no individuals having 0 copies and just 10 samples (1.5%) having 1 copy, with 3-5 copies being the most common ( Table 1).

VALIDATION OF CCL3L1 CN ASSIGNMENTS WITH PARALOG RATIO TEST (PRT) DATA
The CNVrd2 results on 180 of the samples were compared to data where CCL3L1 CN had previously been determined using PRT-45 European samples measured by Carpenter et al. (2011) and 135 African and Asian samples measured by Janyakhantikul et al. (2010). There was 77.8% (140/180) identity between the two methods ( Figure 3A). The majority of discordant results were for high CN samples and differences were of one CN.

COMPARISON OF CNVrd2 WITH OTHER READ DEPTH-BASED METHODS AT CCL3L1
The packages CNVnator (Abyzov et al., 2011) and cn.Mops (Klambauer et al., 2012) were compared to CNVrd2 (refer also to Subjects and Methods). In addition, we also used the 159 genome data set of Sudmant et al. (2010). Copy number for samples in this data set were measured using a read depth-based method, but the authors used the mrsFAST aligner (Hach et al., 2010) which obtains all possible positions of a read. The concordance rates for CNVnator, cn.Mops, and the approach of Sudmant et al. (2010) with CCL3L1 PRT-based results were 7.2% (13/180, r = 0.95), 36.7% (66/180, r = 0.65), and 0% (0/111, r = 0.93) respectively (Figure 3). The highest CNs called by CNVnator and cn.Mops were only 5 and 4 respectively. CN assignments of Sudmant et al. (2010) were significantly higher than the modified PRT results, although assignments of the two methods were highly correlated (r = 0.93) (Figure 3).
To investigate the discordant results of the other packages with the PRT-based method at CCL3L1, observed read-count ratios of European sample sets and all 2535 samples were calculated for both the CCL3L1 gene region and CCL3L1-containing region (Figure 10). The two medians of these ratios were considerably less than 1 (0.7 and 0.8, respectively), while the expected median was 1. For the larger CCL3L1-containing region, 74.7%  (1893/2535) of samples had ratio <1, 76.4% (1938/2535) for the 35kb CCL3L1 region. Reads of samples having zero copies which were aligned to CCL3L1 gene region were very close to zero (Table 4), and thus the density of reads at CCL3L1 was not as high as expected. Approaches based solely on read depth with no cross-sample standardization would therefore have difficulty in accurately assigning CCL3L1 CN, especially for high-CN samples, which would be assigned to lower-CN groups.

POLYMORPHISM OF THE DEFB103A LOCUS
The DEFB103A region lies within a polymorphic inversion with complex structure (Sugawara et al., 2003). We used quantile values and standard deviations across sub-regions to visually identify the boundaries of the polymorphism encompassing the gene. Using a standard deviation threshold of 0.5, DEFB103A was identified within a large CN polymorphic region (consistent with Groth et al., 2008) of 1078 kb (chr8:70200000-7474000 7532001-8098000 including a 58 kb gap) (Figure 7). This CNV region also included small polymorphic blocks which suggested that there would be multiple recombination events inside the CNV region. DEFB103A was in the same block as DEFB104A and DEFB107A. The standard deviations of this block fluctuated around 0.70 and 0.74. We merged these sub-regions and the boundaries of the DEFB103A gene-containing region were chr8:7641001-7742000. The segmentation scores of this region and the larger DEFB103A gene region were strongly correlated (r = 0.95, Figure 8). We used this region to calculate CN for DEFB103A. The DEFB4A gene was outside this block owing to a relatively high peak (SD was 0.81 whilst the maximum SD of the DEFB103A gene-containing region was only 0.74) between DEFB4A and the DEFB103A block (Figure 7). We also calculated segmentation scores of DEFB103B (paralog of DEFB103A) to compare segmentation scores of this gene and those of the DEFB103A gene-containing region. The segmentation scores of the two regions were not as highly correlated (r = 0.83), but they made clear clusters of distinct groups (Figure 11).

ASSIGNMENT OF DEFB103A CN
The segmentation scores of the gene had clear clusters (Figure 9) with some scattered values >3, therefore these values were allocated into the largest (>3) group. The South Asian population had the clearest clusters therefore we used this population to obtain prior information and ran the EM algorithm upon the segmentation score with six groups. The distances between the groups were 0.50-0.69 therefore 0.57 was set as the prior distance. Eight groups were defined for all segmentation scores and CNVrd2 was used to assign DEFB103A copy-number using similar parameters as for CCL3L1 in the Bayesian clustering process ( Table 2). Because the observed read count ratio range was reasonable (0.53-1.93) (Figure 11), this suggested that the lowest CN was greater than zero, with the scattering of values on the left (Figure 9) assigned to the smallest group (CN = 2) during the clustering process. A CN of 4 was the most common in all populations ( Table 2; 26.3-47.9%), consistent with previous reports (Armour et al., 2007;Hardwick et al., 2011).

DISCUSSION
The CNVrd2 package, depicted in Figure 1, is an improved version of CNVrd previously used to assign CN at the FCGR locus  (Nguyen et al., 2013). CNVrd2 was designed to correct issues that can occur during the segmentation process which can result in CNVrd assigning a segmentation score of zero, resulting in some samples being assigned to unreasonable CN groups. In addition, in CNVrd2, a linear regression model was used to adjust the differentiation in CN between populations. Clustering data into different groups is challenging at complex loci with high ranges of CN. Therefore, a Bayesian clustering strategy was integrated into CNVrd2 in order to obtain more reliable results at CCL3L1 and DEFB103A. Finally, plots of percentiles and standard deviations of segmentation results were also added to CNVrd2. These can be used as visual tools to identify copy-number variable regions encompassing loci being measured (Figure 7). In evaluating our results, it should be noted that CNVrd2 was developed on the FCGR and CCL3L1 loci. Therefore, the concordance comparisons are biased toward CNVrd2. However DEFB103A can be regarded as a genuine test locus, on which CNVrd2 performed well (Figure 3). There is a paucity of loci for which published PRT gold-standard data are available on complex CN loci.
CNVrd2 obtains higher concordance with PRT-based methods than other CN-assigning approaches. However the concordance was not 100% (Figure 3; 77.8% at CCL3L1 and 90.4% at DEFB103A), which would limit its use in case-control analysis of complex CN loci in common disease. We used BWA alignments from the 1000 Genomes Project, where the majority of mapping qualities were zero owing to reads having multiple hits within segmental duplications. Further improvement of CNVrd2 could come from the use of alignments better optimized to detect CNV, such as mrsFAST (Hach et al., 2010) which allows diverged paralogs (1-2% difference) to be differentiated. It is important also to stress that CNVrd2 was developed to measure CN at specified loci, while other comparison packages are more automated approaches that both detect CNV and assign CN. CNVrd2 is automated to produce individual segmentation scores at a locus and assign CN using the Bayesian Information Criterion (Schwarz, 1978). However it is better to manually undertake the clustering process as we did in this study. The reason is that clear clusters were not seen at complex loci CCL3L1 and DEFB103A. Note that CNVrd2 is enhanced by some prior knowledge of CN.
From the observed read count ratio analysis (Figure 10), it can be seen that if only the read-count information from one sample is used to assign CN (as is done by CNVnator) the majority of samples at the CCL3L1 locus would likely be called as deletions because of the inherently low number of reads in that region, possibly resulting from CCL3L1 reads mapping to paralogs nearby. The majority of populations had observed read-count ratios of less than 1. This could also be the reason that the cn.Mops package did not call high CN for multiple samples (Figure 3). Complete discordance was seen between the Sudmant et al. (2010) CCL3L1 results obtained using multiple hits of a read and all other methods (Figure 3). The results of Sudmant et al. (2010) had been validated by PCR-based methods with good correlation. It is important to note that PCR-based methods are unreliable in measuring complex gene CN (Walker et al., 2009;McKinney and Merriman, 2012)-any differences in amplification efficiency that alters the ratio between the test and reference genes will alter the apparent CN (Armour et al., 2007). One possible reason for the Sudmant et al. (2010) discordance could be the use of different clustering strategies, but it could also be due to using different alignment methods when there are >1 read-mapping possibilities: choosing a random position (CNVrd2) vs. using all possible positions for a read (Sudmant et al., 2010). We checked the only sample (NA11831) called CN = 1 by Sudmant et al. (2010) and zero by CNVrd2, cn.Mops, CNVnator, and PRT-based methods. There was only one read aligned to CCL3L1 (Table 4), while the average coverage for this sample on the 2 Mb region was adequate (5.8×). Similarly, checking all samples called as 0 copies by CNVrd2, we also saw that the number of reads aligned to the CCL3L1 gene region were extremely low ( Table 4). This suggests that the sample NA11831 was correctly assigned CN = 0 for CCL3L1 by the PRT methods and the BWA aligner based next-generation sequencing methods (CNVrd2, CNVnator, cn.Mops).
At the beta-defensin locus, we identified a CNV region which had nearly equal standard deviations of sub-regions. This CNV region encompassed DEFB103A. The segmentation scores of this and the DEFB103B region (Figure 7; on the polymorphic section before the gap) were moderately correlated (r = 0.83) but they made clear clusters of different groups (Figure 11). This suggested that individual CNs of the two regions would be nearly identical. However, the identified region did not include the DEFB4A gene because we saw a high peak between the region and DEFB4A gene (Figure 7). We retested with different windows (500 and 2000 bp) and still observed this peak (data not shown). An enlarged region that included the peak and DEFB4A was strongly correlated with the smaller region (r = 0.95) but had multiple samples with a score of zero, because in CNVrd2 if a region was segmented into different sub-regions having different signs then the segmentation score of the region was assigned to zero. This suggested that there was not full duplication of the enlarged region (Figure 11). This situation could occur because of an error in the alignment process, or because there were heterogenous breakpoints within the enlarged region.
At the complex CCL3L1 and DEFB103A loci CNVrd2 assigned CN with greater accuracy than the other read depth-based CN assignment methods CNVnator and cn.Mops. The CNVrd2 package is implemented in R (as of version 3.0.1) and is available from Bioconductor (as of version 2.14) at http://www. bioconductor.org/packages/release/bioc/html/CNVrd2.html, where installation and usage instructions can be found.