On the association analysis of CNV data: a fast and robust familybased association method
 Meiling Liu†^{1, 2},
 Sanghoon Moon†^{3},
 Longfei Wang†^{4},
 Sulgi Kim^{5},
 YeonJung Kim^{3},
 Mi Yeong Hwang^{3},
 Young Jin Kim^{3},
 Robert C. Elston^{6},
 BongJo Kim^{3}Email author and
 Sungho Won^{4, 7, 8}Email author
DOI: 10.1186/s128590171622z
© The Author(s). 2017
Received: 15 June 2016
Accepted: 31 March 2017
Published: 18 April 2017
Abstract
Background
Copy number variation (CNV) is known to play an important role in the genetics of complex diseases and several methods have been proposed to detect association of CNV with phenotypes of interest. Statistical methods for CNV association analysis can be categorized into two different strategies. First, the copy number is estimated by maximum likelihood and association of the expected copy number with the phenotype is tested. Second, the observed probe intensity measurements can be directly used to detect association of CNV with the phenotypes of interest.
Results
For each strategy we provide a statistic that can be applied to extended families. The computational efficiency of the proposed methods enables genomewide association analysis and we show with simulation studies that the proposed methods outperform other existing approaches. In particular, we found that the first strategy is always more efficient than the second strategy no matter whether copy numbers for each individual are well identified or not. With the proposed methods, we performed genomewide CNV association analyses of hematological trait, hematocrit, on 521 Korean family samples.
Conclusions
We found that statistical analysis with the expected copy number is more powerful than the statistic with the probe intensity measurements regardless of the accuracy of the estimation of copy numbers.
Keywords
CNV Association analysis Score test HematocritBackground
Copy number variants (CNVs) are widely distributed throughout the human genome [1, 2] and have been considered as important genetic factors for human diseases [3, 4]. Several different methods, such as array comparative genomic hybridization (aCGH) and next generation sequencing, have been suggested to identify CNVs. Thanks to the recent improvement of sequencing technology, sequencing cost decreases very fast and becomes much cheaper. Furthermore, aCGH cannot detect aberrations such as mosaicism that do not result in copy number changes. However, in spite of this advantage of sequencing, aCGH is still cheaper and many array data have already been produced. Thus, it may be a cost effective choice at least for a while. In this report, we focus on CNV analysis with aCGH data–though the proposed method can be readily extended to other types of CNV data.
For aCGH data, gene copy numbers are not directly observed and have to be estimated with their intensity measures for association analyses. True unknown copy numbers will be called as unobserved copy numbers in the remainder of this report. CNV association requires estimation of copy numbers, and several algorithms, such as PennCNV [5], QuantiSNP [6], dChip [7] and GTC [8], have been developed to detect unobserved copy numbers. Then statistical methods such as linear regression and chisquare tests have been utilized to detect CNV association with estimated copy numbers [9]. Barns et al. [10] calculated the posterior probability for each possible copy number, and likelihoods weighted by these posterior probabilities were used to build a likelihood ratio test. As an alternative to CNV analysis using the estimated copy number, the probe intensity measurements can be used to detect the CNV association [11]. The probe intensity is assumed to be proportional to the unobserved copy number, and its distributions can be compared between affected and unaffected individuals. If copy numbers are correctly estimated, the analysis using the expected copy numbers seems to be an efficient choice. However, estimates of copy numbers are often uncertain and this effect has not been carefully considered in the statistical analysis [4]. In this report, we considered both approaches and compared them with simulation studies for a large variety of parameter settings.
For association analysis, phenotypic correlations between individuals have the effect of sample size reduction, and thus independent populationbased samples have often been preferred to maximize the statistical efficiency. However, familybased association analyses have been useful for certain scenarios because family members are genetically homogeneous [12, 13]. For instance, FBATstatistics based on the socalled withinfamily component [14] are robust in the presence of population substructure and they are often preferred, in particular for candidate gene studies. Withinfamily component indicates the distribution of nonfounders’ genotype when their parental genotypes are conditioned. The distribution of founders’ genotype is called betweenfamily component and the statistical power of FBATstatistics has been improved by combining FBATstatistics with the betweenfamily component in a robust way [11, 15, 16]. This twostage analysis can achieve efficiency comparable to that of independent samples. However, due to the assumption that between and within family components are equally informative, this method can suffer from statistical power of loss if the numbers of founders and nonfounders are different.
In this report we propose two statistics, T _{1} and T _{2}, for CNV association analysis using familybased samples; for T _{1}, the phenotypes are regressed on the expected copy number, and for T _{2} they are regressed directly on the probe intensity measurements. A random effect is included to model the phenotypic covariance between family members, and the variance components for the phenotype are estimated with a restricted likelihood. Our results show that statistical analysis with the expected copy number is usually more efficient than the statistic with probe intensity measurements. We applied the proposed methods to detect CNV association with a hematologyrelated trait, hematocrit, in Korean familybased samples.
Methods
Notations and the disease model
We assume that K intensity measurements are observed at a particular CNV region for each individual, there are n families, and n _{ i } individuals in family i. For simplicity, we consider only trio families, but the methods can be extended to large extended families. We assume that j = 1, 2 indicates the parents in each family. We let x _{ ijk } indicate the observed intensity measurement on probe k for individual j in family i. X _{ ij } indicates the column vector, (x _{ ij1},…,x _{ ijK })^{ T }, for individual j in family i. We let λ _{ ij } be the unobserved copy number for individual j in family i, and denote a set of possible realizations of λ _{ ij } and their corresponding frequencies respectively as C and Θ. We denote the phenotype for individual j in family i by y _{ ij }, and let Z _{ ij } be a vector of measured environmental factors, including an intercept as the first element. The intensity matrix, X _{ i }, and phenotype vector, Y _{ i }, for family i are respectively defined as \( {\mathbf{X}}_i={\left({\mathbf{X}}_{i1}^T,\dots, {\mathbf{X}}_{i{ n}_i}^T\right)}^T \) and \( {\mathbf{Y}}_i={\left({y}_{i1},\dots, {y}_{i{ n}_i}\right)}^T \). We include a random effect, b _{ i }, to allow for the phenotypic correlation between family members. λ _{ i } and ε _{ i } indicate respectively an unobserved copy number vector and a measurement error vector for family members in family i. If we let N = Σ_{ i } n _{ i }, an N × K design matrix X and an N × 1 vector Y are respectively obtained by stacking all X _{ i } and Y _{ i } vertically. λ, b and ε are N × 1 vectors and are obtained by stacking all λ _{ i }, b _{ i } and ε _{ i } vertically.
Signal model
The parameters for \( {\boldsymbol{\Sigma}}_{\lambda_{ij}} \) will be denoted by Σ, and this proposed model will be called the signal model in the remainder of this manuscript.
Phenotype model
Copy number model
Parameter estimation with the EM algorithm
To derive a score test for CNV association analysis, β in the phenotype model was assumed to be 0, and the variance component parameters were estimated with the restricted maximum likelihood (REML) method. The copy number vector λ and the random effect vector b are considered as missing variables for the EM algorithm, and the conditional expectation of a complete data loglikelihood was maximized to estimate all the parameters. Individuals were separated with Kmeans clustering [19], and the empirical mean and covariance matrix were used as the initial values for the signal model.
Identifying the number of clusters
If sil _{ ij } is close to one, it indicates that the corresponding individual is wellclustered, whereas if sil _{ ij } is close to −1, it means that the individual is badly clustered. If sil _{ ij } is close to zero, there may exist a better cluster for the corresponding individual. Therefore, we first estimated the copy numbers for each individual for different choices of M. Then we calculated silhouette scores for all the individuals, and the value of M that maximized the mean silhouette score was considered as the optimal choice.
Statistical inference
The detailed derivation of T _{2} is shown in Additional file 1: Text 2. In particular, we can utilize a transformed value for X in T _{2}. For instance, the mean intensity measurement over all probes or the first principal component (PC) score can be utilized, and then T _{2} follows the chisquare distributions with a single degree of freedom. Implementation of the methods is assembled in an R package PedCNV, which is available from CRAN.
Simulation studies
Data generation
Here the dissimilarity between probe intensity measurements in different clusters for probe k is proportional to the value of s _{ k } and we considered three scenarios by using three different choices of s _{ k }: badly separated clusters (BSC), moderately separated clusters (MSC) and well separated clusters (WSC). The different means for the probe intensities were provided by \( {p}_{\lambda_{ij}} \) generated from N(0, 0.9(λ _{ ij } + 1.5)^{2}). The variance of each probe intensity measurement was provided by q _{ v } generated from Γ(0.025,0.0016^{2}). The parameter settings in the signal model are described in Additional file 1: Table S1.
Analysis of a hematological trait
Subjects
Basic characteristics of study participants and hematological trait
Variables  Discovery (family)  Replication (cohort) 

Sample size (n)  521  4694 
Age (years)  38.2 ± 18.3  54.0 ± 9.0 
Male (%)  45.7%  47.1% 
Hematocrit (%)  41.3 ± 4.3  41.1 ± 4.5 
CNV discovery
All samples were genotyped with NimbleGen HD2 3 × 720 K aCGH which contains more than 720,000 probes. Around 360,000 probes were designed based on previously reported CNVs, and the other probes were spaced uniformly throughout the whole genome as a backbone. Sample NA10851obtained from the HapMap lymphoblastoid cell line (LCL) DNA was used as a reference, and NimbleScan version 2.5 was used to process the array image files (.tif) according to the manufacturer’s protocol. Extracted signal intensity was transformed to log2 ratio with hg18/NCBI build 36. Subsequently, we set the log2 ratio thresholds less than −0.25 for a deletion and greater than 0.25 for a duplication, with more than 10 consecutive probes required to assign a CNV.
CNV selection
We used a reciprocal overlap threshold > 50% to find CNVs with similar boundaries for association analysis. According to this threshold, clusters of overlapping CNVs at the sample level are merged into one CNV. Overlapping CNVs with very different sizes and sequentially connected CNVs were excluded from further study. Moreover, we selected CNV clusters which are wellseparated and have multiclass CNVs in order to assign individuals to copynumber classes with high confidence [24]. In total 500 CNVs were utilized for association analyses.
CNV association
PedCNV was applied to an association study with a hematological trait: hematocrit (Hct). The association of CNVs with Hct was analyzed using T _{1} and T _{2}, with age, age^{2} and sex included as covariates. The resulting statistics were adjusted by using genomic control to allow for population substructure.
CNV validation by PCR experiment
To confirm CNV genotypes, a PCR using the AccuPrime Taq DNA Polymerase High Fidelity (invitrogen, CA, USA) was performed on 10–16 individuals selected from each cluster (Additional file 1: Figure S1 (A)). The primers were designed to give rise to amplicons with different lengths to detect both the deleted (690 bp) and normal (1519 bp) alleles (Additional file 1: Table S2). Genomic locations for designed primers based on human genome assembly hg18 were converted to those based on hg19 by liftOver of the UCSC genome browser. PCR was carried out on a GeneAmp PCR system 9700 (Applied Biosystems, Calif., USA) with the following PCR conditions: 5 min at 95 °C, followed by 33 cycles of 30s at 95 °C, 30s at 60 °C, 2 min at 68 °C, and final extension at 68 °C for 7 min. The resulting PCR products were visualized by electrophoresis separation on a 1.5% agarose gel with SafePinky DNA gel staining solution (Genedepot, TX, USA). Moreover, to confirm exact breakpoints of the CNVs, PCR products were sequenced using an ABI 3730 DNA analyzer (Applied Biosystems, CA, USA).
Replication study
We have previously implemented KGVDB, which includes 3601 multiclass CNVs and their tagging SNPs, from 4694 communitybased cohort samples, as a part of the Korean Genome Epidemiology Study (KoGES) [25]. We used these unrelated individual samples to pursue replication of the identified CNV from the discovery association study. Table 1 shows a summary of the participants’ characteristics. In short, all the 4694 samples were also genotyped with NimbleGen HD2 3 × 720 K aCGH. The NA10851 sample was again used as a reference. NimbleScan version 2.5 was used to extract signal intensity. Subsequently, quality control, such as normalization and waviness correction, was conducted using the R package (http://cran.rproject.org) and WaveNorm [26]. For CNV detection, the Genome Alteration Detection Analysis algorithm (GADA) was used with T = 10, alpha = 0.2 and MinSegLen = 10. Moreover, an average log2 ratio of ±0.25 was set as a cutoff value [25]. Among the detected CNVs, we selected those CNVs having a similar boundary with any CNV significant in the discovery association study. Additional file 1: Figure S2 shows the overall process of the replication study.
CNV validation of replication study samples
To verify whether an estimated CNV genotype using cohort samples is true or not, we carried out quantitative PCR (qPCR) using the TaqMan Copy number Assay (Life Technologies, Foster City, CA, USA) according to the manufacturer’s guidelines. A predesigned TaqMan probe (Assay ID: Hs04965547_cn) was used to validate the existence of the CNV. All experiments were replicated three times to enhance the validation accuracy. The samples used for validation were randomly selected from each genotype (Additional file 1: Figure S1 (B)). Copy number genotype for each sample was calculated by Copy caller v2.0 (Applied Biosystems, Calif., USA) using the manufacturer’s guideline.
Results
Evaluation with simulated data
Clustering
With the simulated data we evaluated the accuracy of estimating M and the estimated copy numbers for each individual when the true M was assumed to be known. The results from the proposed method were compared with CNVtools [10]. The probe intensity measurements were generated under the three different scenarios: BSC, MSC and WSC. For each individual, we calculated from the probe intensity measurements the mean, the first PC score and the fewest PC scores that explain more than 90% of the variation; they are denoted by mean, PC1 and PC.9 respectively. In addition to the original probe intensity measurements (RAW), we used the mean, PC1 and PC.9, for the proposed method and the results were compared with CNVtools. For CNVtools, the mean, PC1 and the onedimensional canonical correlation transformed vector of the probe intensity measurements were used.
Additional file 1: Tables S3 and S4 show the accuracy of the estimated value of M from 1000 replicates using PedCNV and CNVtools. The proposed method using PC1 was always the most accurate, followed by the proposed method using PC.9. The results from the proposed method performed better than CNVtools. CNVtools had a tendency to choose a larger number of clusters, and the results were rarely consistent, even when the clusters were well separated. CNVtools selects the number of clusters using a Bayesian information criterion [10], while the proposed method selects it with a silhouette score, which appears to be a better choice. In Additional file 1: Table S5, M was set to be the true value 3 for all methods, and the relative proportions of individuals for whom the estimated copy number was consistent with the true copy number were calculated from 1000 replicates under the null and alternative hypotheses. Additional file 1: Table S5 shows that the proposed method based on PC1 was the most accurate, followed by the proposed method based on PC.9. Therefore we conclude that basing our method on PC1 may be a reasonable choice.
Association analysis
Empirical type 1 error estimates (M = 3)
Significance Level  

.005  .05  .1  .2  
BSC  T1  0.0060 ± 0.0021  0.0504 ± 0.0061  0.1018 ± 0.0084  0.2082 ± 0.0113 
T2  0.0056 ± 0.0021  0.0550 ± 0.0063  0.1008 ± 0.0086  0.2072 ± 0.0112  
MSC  T1  0.0048 ± 0.0019  0.0486 ± 0.0060  0.1006 ± 0.0083  0.2104 ± 0.0113 
T2  0.0048 ± 0.0019  0.0472 ± 0.0059  0.0956 ± 0.0082  0.1884 ± 0.0108  
WSC  T1  0.0056 ± 0.0021  0.0516 ± 0.0061  0.0962 ± 0.0082  0.2006 ± 0.0111 
T2  0.0048 ± 0.0019  0.0498 ± 0.0060  0.0968 ± 0.0082  0.1922 ± 0.0109 
Empirical power estimates (M = 3)
Significance Level  β  

.1  .2  .3  .4  .5  .6  
.001  BSC  T1  0.0135  0.1390  0.4830  0.8410  0.9830  0.9985 
T2  0.0065  0.0510  0.2370  0.5745  0.8860  0.9795  
FBAT  5e4  0.0060  0.0340  0.1300  0.3570  0.6005  
MSC  T1  0.0160  0.1570  0.5530  0.8740  0.9885  1.0000  
T2  0.0085  0.0685  0.3200  0.6900  0.9290  0.9945  
FBAT  0.0000  0.0100  0.0695  0.2505  0.5575  0.8385  
WSC  T1  0.0195  0.1615  0.5375  0.8935  0.9910  0.9980  
T2  0.0075  0.0815  0.3300  0.7240  0.9545  0.9970  
FBAT  0.0010  0.0115  0.0915  0.3240  0.6460  0.8955  
.01  BSC  T1  0.0710  0.3585  0.7510  0.9605  0.9990  1.0000 
T2  0.0265  0.1780  0.4630  0.7900  0.9685  0.9950  
FBAT  0.0155  0.0455  0.1540  0.3775  0.6620  0.8430  
MSC  T1  0.0725  0.3805  0.8070  0.9690  0.9990  1.0000  
T2  0.0340  0.2100  0.5640  0.8660  0.9790  0.9980  
FBAT  0.0150  0.0550  0.2400  0.5385  0.8155  0.9645  
WSC  T1  0.0800  0.3795  0.7925  0.9740  0.9985  1.0000  
T2  0.0370  0.2115  0.5730  0.8855  0.9920  0.9995  
FBAT  0.0175  0.0710  0.2740  0.6350  0.8775  0.9740  
.05  BSC  T1  0.1905  0.5930  0.9075  0.9920  1.0000  1.0000 
T2  0.0975  0.3505  0.6880  0.9080  0.9910  0.9990  
FBAT  0.0575  0.9610  0.3660  0.6310  0.8555  0.9525  
MSC  T1  0.2050  0.6190  0.9295  0.9900  1.0000  1.0000  
T2  0.1110  0.3915  0.7650  0.9495  0.9970  1.0000  
FBAT  0.0725  0.2080  0.4990  0.7585  0.9305  0.993  
WSC  T1  0.2095  0.6145  0.9260  0.9950  0.9995  1.0000  
T2  0.1255  0.3995  0.7650  0.9530  0.9990  1.0000  
FBAT  0.0790  0.2240  0.5460  0.8415  0.9705  0.9960 
Results of real data analysis
CNV association
The most significant results of T _{ 1 } and T _{ 2 } from analyzing the family data
Chr  Position  0/1/2  T _{1}  T _{2} 

8  94141469–94142527  42/226/253  1.38e03  4.67e02 
5  147534018–147534337  119/265/137  1.19e02  3.92e03 
CNV validation
Among 500 multiclass CNVs, the CNV region (chr7:81279592–81280418) was randomly selected for evaluation of CNV genotype estimation. In total, 41 subjects were selected from each CNV cluster and a PCR experiment was conducted for them. Among these samples, 38 subjects (92.7%) had the same copy numbers as the estimates from the proposed methods (Additional file 1: Figures S3 and S4).
Discussion
Even though CNV has been expected to be an important genetic factor for many diseases, CNV association analysis has often been limited because of uncertainty of the copy number, and several statistical methods [8, 27] have been proposed to handle this uncertainty. However, even though some of the existing methods are relatively accurate, the estimated copy numbers are not accurate in some situations, which might cause a power loss for CNV association analysis. In this report, we propose new statistical methods for CNV association analysis with familybased samples. With extensive simulations, we showed that the proposed methods perform much better than the existing approaches. The proposed method was implemented in the R package, pedCNV and the main function in our R package was implemented with C++. We found that association analyses of 300 trios were completed within one minute using an Intel (R) Xeon (R) E52620 0 CPU at 2.00GHz, with a single node and 80 gigabyte memory.
Furthermore, the proposed method is flexible and can be extended to various scenarios. First, the proposed methods consist of T _{1} and T _{2}. The former is based on the estimated copy number and the latter is on the probe intensity measurements. Our simulation studies show that the most efficient statistic is always the statistic with the expected copy numbers. However, if the accuracy of the estimated copy numbers is not clear and there is a systematic bias, the statistical power of T _{1} can be substantially affected, and some modification can be made to the proposed methods. For instance, the minimum of the pvalues for T _{1} and T _{2} could be considered as a test statistic and permutationbased pvalues could be calculated. Alternatively, the posterior probabilities for each copy number estimated from the E step in T _{1} can be utilized as classified copy numbers. These modifications are computationally feasible and may provide less sensitive results compared to T _{1} and T _{2}. Second, the presence of population substructure has been known to be a factor that leads to violation of the assumptions underlying statistical association analysis. In our real data analysis, the genomic control approach [28] was adopted, but the linear mixed model is known to be the most efficient if the polygenic effects are substantial [29]. The correlations between individuals can be estimated with largescale genetic data such as genomewide SNPs, and this can be incorporated into the phenotype model in the proposed method. Third, the proposed methods can be simply extended to the sequencing data with a minor modification even though it only applied to aCGH data in this report. This will be investigated in our future work.
However, in spite of the practical advantage of the proposed methods, there exist some limitations, and further investigation is necessary. First, the incorporation of Mendelian transmission into the signal model induces a substantial computational burden for large families. In our PedCNV package, Mendelian transmission for a signal model is considered, but only for nuclear families. We found with simulation studies that the drop of accuracy is not substantial when Mendelian transmission is not considered, but its effect can be substantial if only a few large families are available. A peeling algorithm [30] has been developed that minimizes the computation of likelihoods for large families and it will be implemented in the PedCNV package. Second, the proposed method assumes that there is no de novo mutation and recombination. In such cases, the statistic T _{2} may be a better choice. Third, it has been observed that the bias in CNV calls can be different between parents and offspring, and our first statistic, T _{1}, can suffer from this differential bias. Our simulation studies do not examine any such violation of statistical assumptions, but its effect on T _{1} could be substantial in CNV association analysis with large families. Third, copy numbers for each individual were identified by calculating the expectation of copy numbers using the posterior probability and the expected copy numbers were utilized as λ in T _{1}. Although this maximum likelihood approach for classification can yield inconsistent estimators of parameters [31, 32], the simulation studies show that the accuracy of this method is higher. Thus we continued adopting this method in spite of its deficiencies.
In recent decades various types of genetic data have been used to detect the genetic factors underlying many diseases and many disease susceptibility loci have been found. Even though CNVs have been expected to be an important genetic factor, the findings of CNV association analysis have been limited and the proposed methods may bridge this gap by alleviating the issue of copy number uncertainty.
Conclusion
PedCNV presents a computationally efficient R package that provides two statistics for familybased CNV association analysis: first, the copy number is estimated by maximum likelihood and association of the estimated copy number with the phenotype is tested; second, the observed probe intensity measurements is directly used to detect association of CNV with the phenotypes of. The simulation studies showed that the proposed methods outperform other existing approaches. In particular, we found that statistical analysis with the expected copy number is more powerful than the statistic with the probe intensity measurements regardless of the accuracy of the estimation of copy numbers.
Abbreviations
 aCGH:

array comparative genomic hybridization
 BSC:

Badly separated clusters
 CNV:

Copy number variation
 GADA:

Genome Alteration Detection Analysis
 Hct:

Hematocrit
 KoGES:

Korean Genome Epidemiology Study
 LCL:

Lymphoblastoid cell line
 MSC:

Moderately separated clusters
 PC:

Principle component
 qPCR:

quantitative PCR
 QQ plot:

Quantile quantile plot
 REML:

Restricted maximum likelihood
 WSC:

Well separated clusters
Declarations
Acknowledgements
We do appreciate that Korean Genome Analysis Project (4845301) provided the genome data, and that all three anonymous reviewers and editor for valuable suggestions and constructive comments, which greatly help us to improve our manuscript.
Funding
This study was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by National Research Foundation of Korea Grant by the Korean Government (NRF2014S1A2A2028559) and the Ministry of Education (NRF2016R1D1A1A09919610), and an intramural grant from the Korea National Institute of Health (20127300400). The funding body was not involved in the design or conclusion of our study.
Availability of data and materials
Our R package, PedCNV, can be downloaded from https://cran.rproject.org/web/packages/PedCNV/index.html. Simulated data with 2000 replicates for different values of β under WSC when three copy number clusters are considered, can be downloaded as Additional file 2. In the real data analysis, we used the same DNA samples as were used in Lee et al. [23].
Authors’ contributions
Conceived and designed the experiments: ML and SW. Performed the experiments: ML and LW. Analyzed and interpreted the data: ML, SM, LW and SW. Drafted the manuscript: ML, SM, LW and SW. Edited the manuscript: SK, YJK, MYH, YJK, RCE and BJK. All authors read and approved the final version of the manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
All subjects from this study provided written informed consent and the institutional ethics committees of participating institutions approved the experimental protocols (approved IRB number: 201108CON10P).
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Sharp AJ, Cheng Z, Eichler EE. Structural variation of the human genome. Annu Rev Genomics Hum Genet. 2006;7:407–42.View ArticlePubMedGoogle Scholar
 Feuk L, Carson AR, Scherer SW. Structural variation in the human genome. Nat Rev Genet. 2006;7(2):85–97.View ArticlePubMedGoogle Scholar
 Lupski J. Genomic Disorders: The Genomic Basis of Disease. J Med Genet. 2008;45:S32.Google Scholar
 McCarroll SA, Altshuler DM. Copynumber variation and association studies of human disease. Nat Genet. 2007;39:S37–42.View ArticlePubMedGoogle Scholar
 Wang K, Li MY, Hadley D, Liu R, Glessner J, Grant SFA, Hakonarson H, Bucan M. PennCNV: An integrated hidden Markov model designed for highresolution copy number variation detection in wholegenome SNP genotyping data. Genome Res. 2007;17(11):1665–74.View ArticlePubMedPubMed CentralGoogle Scholar
 Colella S, Yau C, Taylor JM, Mirza G, Butler H, Clouston P, Bassett AS, Seller A, Holmes CC, Ragoussis J. QuantiSNP: an Objective Bayes HiddenMarkov Model to detect and accurately map copy number variation using SNP genotyping data. Nucleic Acids Res. 2007;35(6):2013–25.View ArticlePubMedPubMed CentralGoogle Scholar
 Li C. Automating dChip: toward reproducible sharing of microarray data analysis. BMC Bioinformatics. 2008;9:231.View ArticlePubMedPubMed CentralGoogle Scholar
 Korn JM, Kuruvilla FG, McCarroll SA, Wysoker A, Nemesh J, Cawley S, Hubbell E, Veitch J, Collins PJ, Darvishi K, et al. Integrated genotype calling and association analysis of SNPs, common copy number polymorphisms and rare CNVs. Nat Genet. 2008;40(10):1253–60.View ArticlePubMedPubMed CentralGoogle Scholar
 Kim JH, Hu HJ, Yim SH, Bae JS, Kim SY, Chung YJ. CNVRuler: a copy number variationbased casecontrol association analysis tool. Bioinformatics. 2012;28(13):1790–2.View ArticlePubMedGoogle Scholar
 Barnes C, Plagnol V, Fitzgerald T, Redon R, Marchini J, Clayton D, Hurles ME. A robust statistical method for casecontrol association testing with copy number variation. Nat Genet. 2008;40(10):1245–52.View ArticlePubMedPubMed CentralGoogle Scholar
 IonitaLaza I, Perry GH, Raby BA, Klanderman B, Lee C, Laird NM, Weiss ST, Lange C. On the analysis of copynumber variations in genomewide association studies: a translation of the familybased association test. Genet Epidemiol. 2008;32(3):273–84.View ArticlePubMedGoogle Scholar
 Bansal V, Libiger O, Torkamani A, Schork NJ. Statistical analysis strategies for association studies involving rare variants. Nat Rev Genet. 2010;11(11):773–85.View ArticlePubMedPubMed CentralGoogle Scholar
 Shi G, Rao DC. Optimum Designs for NextGeneration Sequencing to Discover Rare Variants for Common Complex Disease. Genet Epidemiol. 2011;35(6):572–9.PubMedPubMed CentralGoogle Scholar
 Laird NM, Horvath S, Xu X. Implementing a unified approach to familybased tests of association. Genet Epidemiol. 2000;19 Suppl 1:S36–42.View ArticlePubMedGoogle Scholar
 Murphy A, Won S, Rogers A, Chu JH, Raby BA, Lange C. On the genomewide analysis of copy number variants in familybased designs: methods for combining familybased and populationbased information for testing dichotomous or quantitative traits, or completely ascertained samples. Genet Epidemiol. 2010;34(6):582–90.View ArticlePubMedPubMed CentralGoogle Scholar
 Won S, Wilk JB, Mathias RA, O’Donnell CJ, Silverman EK, Barnes K, O’Connor GT, Weiss ST, Lange C. On the analysis of genomewide association studies in familybased designs: a universal, robust analysis approach and an application to four genomewide association studies. PLoS Genet. 2009;5(11):e1000741.View ArticlePubMedPubMed CentralGoogle Scholar
 Thornton T, McPeek MS. ROADTRIPS: Casecontrol Association Testing with Partially or Completely Unknown Population and Pedigree Structure. Am J Hum Genet. 2010;86(2):172–84.View ArticlePubMedPubMed CentralGoogle Scholar
 Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genomewide association studies. Nat Genet. 2006;38(8):904–9.View ArticlePubMedGoogle Scholar
 Hartigan JA. Clustering algorithms. New York: Wiley; 1975.Google Scholar
 McPeek MS, Wu XD, Ober C. Best linear unbiased allelefrequency estimation in complex pedigrees. Biometrics. 2004;60(2):359–67.View ArticlePubMedGoogle Scholar
 Rousseeuw PJ. Silhouettes–a Graphical Aid to the Interpretation and Validation of ClusterAnalysis. J Comput Appl Math. 1987;20:53–65.View ArticleGoogle Scholar
 Rao CR, Mitra SK. Generalized inverse of matrices and its applications. New York: Wiley; 1971.Google Scholar
 Lee EJ, Kim KJ, Kim HN, Bok J, Jung SC, Kim EK, Lee JY, Kim HL. Genomewide scan of granular corneal dystrophy, type II: confirmation of chromosome 5q31 and identification of new cosegregated loci on chromosome 3q26.3. Exp Mol Med. 2011;43(7):393–400.View ArticlePubMedPubMed CentralGoogle Scholar
 Craddock N, Hurles ME, Cardin N, Pearson RD, Plagnol V, Robson S, Vukcevic D, Barnes C, Conrad DF, Giannoulatou E, et al. Genomewide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls. Nature. 2010;464(7289):713–20.View ArticlePubMedGoogle Scholar
 Moon S, Jung KS, Kim YJ, Hwang MY, Han K, Lee JY, Park K, Kim BJ. KGVDB: a populationbased genomic map of CNVs tagged by SNPs in Koreans. Bioinformatics. 2013;29(11):1481–3.View ArticlePubMedPubMed CentralGoogle Scholar
 Marioni JC, Thorne NP, Valsesia A, Fitzgerald T, Redon R, Fiegler H, Andrews TD, Stranger BE, Lynch AG, Dermitzakis ET, et al. Breaking the waves: improved detection of copy number variation from microarraybased comparative genomic hybridization. Genome Biol. 2007;8(10):R228.View ArticlePubMedPubMed CentralGoogle Scholar
 Fiegler H, Redon R, Andrews D, Scott C, Andrews R, Carder C, Clark R, Dovey O, Ellis P, Feuk L, et al. Accurate and reliable highthroughput detection of copy number variation in the human genome. Genome Res. 2006;16(12):1566–74.View ArticlePubMedPubMed CentralGoogle Scholar
 Devlin B, Roeder K. Genomic control for association studies. Biometrics. 1999;55(4):997–1004.View ArticlePubMedGoogle Scholar
 Kang HM, Zaitlen NA, Wade CM, Kirby A, Heckerman D, Daly MJ, Eskin E. Efficient control of population structure in model organism association mapping. Genetics. 2008;178(3):1709–23.View ArticlePubMedPubMed CentralGoogle Scholar
 Elston RC, Stewart J. A general model for the genetic analysis of pedigree data. Hum Hered. 1971;21(6):523–42.View ArticlePubMedGoogle Scholar
 Bryant BG. Largesample results for optimizationbased clustering methods. J Classif. 1991;8(1):31–44.View ArticleGoogle Scholar
 Bryant PG, Williamson JA. Asymptotic Behaviour of Classification Maximum Likelihood Estimates. Biometrika. 1978;65:273–81.View ArticleGoogle Scholar