Copy number variation in the Framingham Heart Study

One of the most important challenges facing the field of biology today is making sense of genetic variation. One form of variation that is gaining an increasing level of attention is copy number variation (CNV). The term CNV is used to encompass a variety of polymorphisms that occur at the sequence level and that affect the copy number of regions of genetic material. Examples include insertions, deletions, rearrangements, and duplications. While humans have CNV, the degree of that CNV is only now becoming appreciated [1, 2], and its influences on phenotypic variation is not yet well understood. However, while a number of associations between single-nucleotide polymorphisms (SNPs) and disease have been verified, it remains the case that the effect sizes of these polymorphisms appear relatively small, with reported odds ratios typically being below 1.5 (see Estivill and Armengol [3] for a recent review). While the reasons for this remain unclear, it is plausible that SNP data are not capturing all the relevant associations. Consequently, in this paper we implement population-based tests for association between CNV and phenotypic variation within the Framingham Heart Study samples provided as part of Genetic Analysis Workshop 16 Problem 2.

Our analysis currently focuses on a subset of Offspring Cohort individuals with diabetes. In order to protect against undiagnosed diabetes, we insisted that cases had a fasting plasma glucose measure of at least 126 mg/dl, while controls had a fasting measure less than 110 mg/dl. The controls were also frequency-matched to cases by 5-year age intervals, sex, and ever smoked status (age and smoking history taken at baseline). Some individuals are related, but the majority are not (at least by the pedigree information provided). This resulted in a final sample size of 194 cases and 213 controls, for which we analyzed the 500 k SNP data. We then called copy number. Space prohibits a display of estimated copy number here, but, in summary, we observed that previously known sites of copy number variations [2] were detected as such in this data set, and that, for each method, there is clear correlation between the locations at which copy number changes across individuals (despite the fact that the algorithm call copy number independently for each individual). As expected, the number of changes of copy number in a gene was directly proportional to its length (correlation coefficient = 0.89).

Table 1 presents a summary of the results for the GADA and CNAM analyses for each chromosome. No gene attains genome-wide significance when tested for association between copy number and diabetes. However, some interesting features are observed. For example, under the null hypothesis of no association, p-values should be uniformly distributed, with each gene having a probability of 0.05 of being reported as 'significant' on this basis. Thus, for each chromosome, the number of genes showing association at the level p < 0.05 will be binomially distributed. We used this to test whether there was an over- or under-representation of such associated genes on any chromosome. After a multiple comparison correction, several chromosomes showed evidence of excess, or lack, of associated genes. Furthermore, there is a clear excess of small p-values for the CNAM analysis. This can be seen in Figure 1, where we show a quartile-quartile (Q-Q) plot of observed vs. expected p-values genome-wide (we plot -log₁₀ of the p-values). The Q-Q plots reflect the relative over- or under-abundance of small p-values resulting from the CNAM or GADA analysis. In Figure 2 we show how the distribution of the genes with the smallest p-values when tested for association with diabetes varies along the genome for the GADA and CNAM analyses. There appears to be little agreement between the two methods. This raises the question of why the results of the methods show such poor agreement among these genes. Recall that the first step for both methods is the normalization of SNP intensities, but that the two methods employ different normalization techniques. We therefore examined the degree of agreement between the SNP intensities after normalization. We show this in Figure 3 (left side), a scatter plot of the intensities resulting from the two schemes. The results are striking: there is very poor correlation between the normalized intensities. Thus, even if the methods were using identical routines to call copy number, which they are not, we would expect to obtain widely different calls of copy number, with consequent large differences between p-values resulting from subsequent tests of association between copy number and phenotype. We show the poor agreement in copy number in Figure 3 (right side). It is important to note that both methods use somewhat arbitrary definitions of upper and lower cut-points (CU and CL, respectively), defined such that if the normalized intensity for a SNP is below CL it is determined to have a copy number less than two, whereas if normalized intensity is above CU the copy number is called as greater than two; otherwise the copy number is called as two. These arbitrary cut-points are different for the two methods and, combined with the differing results from normalization, result in wildly different calls of copy number. Most often both methods call copy number as two, but it is seldom the case that the methods simultaneously determine copy number to be other than two for any given SNP. This, clearly, is a concern.

			% Significant genesc
Chr	Affymetrix genesa	Genotyped genesb	GADA	CNAM
1	1650	1500	5.12	12.60
2	1046	974	0.72	16.32
3	893	829	0.36	13.51
4	661	634	1.31	1.57
5	722	680	7.45	22.35
6	883	822	1.84	18.86
7	709	655	3.09	6.87
8	539	505	10.18	1.78
9	632	582	0.87	4.81
10	636	603	1.84	8.96
11	1028	907	1.68	4.85
12	871	799	4.80	15.14
13	295	287	8.81	10.45
14	527	474	4.25	12.66
15	507	474	1.92	9.92
16	580	483	5.64	3.52
17	851	712	8.09	0.14
18	241	241	4.17	1.66
19	912	713	1.51	0.14
20	470	439	2.75	5.47
21	201	187	0	9.09
22	365	326	3.12	0
Total/Average	15219	13826	3.61	8.62

^aOverall number of genes as defined by Affymetrix annotation files

^bNumber of those genes for which we had SNP intensity data

^cPercentage of those genes in which associated CNV was detected at the p = 0.05 value (uncorrected for multiple comparisons) using GADA and CNAM

The importance of CNV has only recently become appreciated. The relationship of CNV to phenotypic variation is even less well developed. While we hope the analysis presented in this paper is a useful step forward in this area, it merely scratches the surface of what is likely to be an extremely complex challenge. CNV lacks some of the 'neatness' of SNP data. It does not occur at well defined positions (i.e., the points at which CNV changes is often different across individuals). Furthermore, for a variety of reasons, the 500 k chip platforms are not ideal for detecting CNV when compared with more recent platforms. It is also the case that many functional mutations occur outside genes (in promoter regions, for example). As such, many regions of CNV may not be detected. It is also likely to be challenging to detect small CNVs using these technologies.

An approach in which we break the genome into regions of maximal length, such that copy number remains constant for each individual within each region, results in an extremely large number of regions (around 82,000 regions for the samples analyzed in this paper). While each region can be treated as if it were a (multi-allelic) locus, and marginal tests can then be performed, such tests are likely to be far from optimal.

The principal reason for this is the correlation between intervals. This is directly analogous to the situation with SNP analysis, but may well be even more complex in this new setting. An alternative approach might be to attempt to relate SNP variation to CNV, and adopt an approach akin to the tag-SNP idea.

In this paper we attempted to move past the perils of a strategy based upon marginal tests by conducting a gene-at-a-time analysis, relating mean copy number within a gene to phenotype. Such an approach seems reasonable, succeeds in reducing the analysis to a reasonable number of tests, and thereby avoids the worst excesses of multiple comparison corrections. However, it should be noted that the very poor agreement in results from the two methods explored in this paper indicates that substantial work remains to be done. It is this lack of agreement that represents the principal lesson to be drawn from the current study. Normalization is regarded as an important, but somewhat routine step in analyses such as these. However, our paper demonstrates that the particular method of normalization chosen can have a key influence on the results obtained. In our case, the two normalization methods are both widely used, and appear inherently sensible, but result in normalized intensities that are very poorly correlated across methods. Consequently, subsequent analyses will produce wildly different results. As the phrase "garbage in, garbage out" reminds us, it is important to ensure that such normalization routines are adding signal, rather than noise to the data. As such, there is an urgent need for a widespread comparison of normalization methods in order to better assess which of them is most effective.

Finally, it should also be noted that, while we look at all genes in the present study, there is little reason a priori to expect candidate genes chosen on the basis of a SNP study to also have a function due to CNV. It is entirely possible that genes that affect phenotype through CNV will be distinct from those that have an effect due to SNP polymorphism. In a recent study, Stranger et al. examined data from Phase 1 of HapMap and noted that SNPs and CNVs captured 83.6% and 17.7%, respectively, of the total detected genetic variation in the expression of around 15,000 genes, but that "the signals from the two types of variation had little overlap" [7]. Other studies take a different view (e.g., McCarroll et al. [8]).