Abstract
Background
The frequency of a haplotype comprising one allele at each of two loci can be expressed as a cubic equation (the 'Hill equation'), the solution of which gives that frequency. Most haplotype and linkage disequilibrium analysis programs use iterationbased algorithms which substitute an estimate of haplotype frequency into the equation, producing a new estimate which is repeatedly fed back into the equation until the values converge to a maximum likelihood estimate (expectationmaximisation).
Results
We present a program, "CubeX", which calculates the biologically possible exact solution(s) and provides estimated haplotype frequencies, D', r^{2 }and χ^{2 }values for each. CubeX provides a "complete" analysis of haplotype frequencies and linkage disequilibrium for a pair of biallelic markers under situations where sampling variation and genotyping errors distort sample HardyWeinberg equilibrium, potentially causing more than one biologically possible solution. We also present an analysis of simulations and real data using the algebraically exact solution, which indicates that under perfect sample HardyWeinberg equilibrium there is only one biologically possible solution, but that under other conditions there may be more.
Conclusion
Our analyses demonstrate that lower allele frequencies, lower sample numbers, population stratification and a possible D' value of 1 are particularly susceptible to distortion of sample HardyWeinberg equilibrium, which has significant implications for calculation of linkage disequilibrium in small sample sizes (eg HapMap) and rarer alleles (eg paucimorphisms, q < 0.05) that may have particular disease relevance and require improved approaches for meaningful evaluation.
Background
Linkage disequilibrium (LD) describes the condition that occurs when alleles at different loci are nonrandomly associated in a given population. Under LD the frequency (f_{11}) of a haplotype (h_{11}) representing the "1" allele at two loci is significantly more or less than the product of the respective allele frequencies. Characterisation of LD is important in medical genetics, influencing association mapping of trait loci and providing information on interactions between genes [1,2]. LD is the result of a shared history of mutation and recombination, and other factors including: genetic drift, population growth, admixture, population structure, the ages of the polymorphisms, the physical distance separating them and the effects of selective pressure [3].
For unrelated individuals the estimation of LD relies on the estimation of haplotype
frequencies. In a 3 × 3 table for a biallelic marker the haplotype phase of all individuals
is known with the exception of the centre cell (representing individuals heterozygous
at both loci). The estimated frequency,
that is adapted from Hill's equation (4) [4] with the constants defined under Methods. With
Several approaches exist for solving equation (1), the solution of which enables estimation of haplotype frequencies and LD coefficients. The first approach uses iterationbased algorithms. An initial estimate of haplotype frequency (either random, or based on the known haplotype numbers) is substituted into the equation, providing a new estimate. This is then fed back into the equation and the expectationmaximisation (EM) process repeated until the values converge. This is the basis both of the algorithm described by Hill in 1974 for the estimation of pairwise haplotype frequencies [4], and of other EM algorithms that enable the estimation of multilocus haplotype frequencies. Many programs exist that utilise variations on this approach, including: GOLD [5], GOLDsurfer [6], MIDAS [7], Haploview [8] and many others reviewed in [912]. The potential problem for these approaches is that algorithms may converge on one of the alternative roots of the cubic equation (a local maximum rather than the global maximum).
Other approaches include parsimony, eg HAPAR [13] and Bayesian algorithms, eg PHASE [1416]. Parsimony and Bayesian methods are both better suited to estimating individual haplotypes than EM approaches, while Bayesian and EM methods are useful for estimating population frequencies [11].
An alternative approach would be exact solution, such as Cardan's solution [17] of the generalized cubic equation (of which equation (1) is an example). This provides all roots to the cubic equation, from which we can select those that are both real (i.e. not a complex number) and biologically possible. If more than one solution exists then the likelihoods of the different solutions can be compared and an informed evaluation made of the result. Theoretically, the noniterative approach may be computationally less intensive and more accurate, but computational efficiency and accuracy will be software and platform dependent.
Implementation
Hill assumed random mating and Hardy Weinberg Equilibrium (HWE) [4]. Rearranging terms for consequent diplotype frequency expectations for two biallelic
loci Luo and Suhai [18] obtained equation 1 given in the introduction (here redefining
We have adopted the Nickalls treatment of the Cardan solution of the generalized cubic equation [17], and written a Python [19] program "CubeX" to solve equation 1 exactly. In CubeX, after calculation of constants ad from diplotypic data the following are calculated:
The discriminant Δ_{3 }= y_{N}^{2 } h^{2 }is then used to determine the outcome in real roots (without having to go through complex number intermediates or ambiguities), with three possible outcomes:
Outcome 1: if y_{N}^{2 }> h^{2 }there will be only one real root (α) given by
Outcome 2: if y_{N}^{2 }= h^{2 }there are three real roots (α, β and γ) and α and β are equal. For a value of
Outcome 3: if y_{N}^{2 }<h^{2 }there are three real roots (α, β and γ). Where
Values for D' and r^{2 }are calculated as previously described [20,21]:
D_{max }= min [p(1q),(1p)q] if D > 0 or D_{max }= min [pq, (1p)(1q)] if D < 0
Diplotype frequencies based on the estimated haplotype frequencies are compared to the input diplotype frequencies by a χ^{2 }test, which effectively tests sample deviation from the null hypothesis of HWE for the diplotypes formed of the four haplotypes. The number of degrees of freedom is equal to the number of observations (diplotype counts) minus four estimated parameters which are either three haplotypes (the fourth can be inferred) and D, or one haplotype, two allele frequencies and D. If nine different diplotypes are observed the number of degrees of freedom is therefore five. For each empty cell in the 3 × 3 the number of degrees of freedom is reduced by one. If the user knows there are only three haplotypes present (and therefore six diplotypes) then there are only three estimated parameters (D is inferred by the three haplotype frequencies) and 3 df. It is important to note that in the latter case neither cubic solution nor iteration is necessary as the haplotype frequencies can be directly counted from the diplotype data. If the user believes that there are only three alleles and hence six diplotypes, but there are nonzero values for any of the other three possible diplotypes, then reconsideration of the technical veracity of the data and of the homogeneity of the population sample would be wise.
Results
Solutions are considered biologically possible when
Number of solutions to the cubic equation with simulated data
We have calculated the number of possible solutions to the cubic equation for genotypes of simulated pairs of SNPs with a range of allele frequencies for a range of sample sizes. The genotype numbers were calculated assuming HWE with a wide range of LD situations for the two SNPs. This was achieved by simulating all combinations of haplotype frequencies between 0 and 1, at intervals of 1/55, that add up to 1. These haplotype frequencies were then converted to diplotype frequencies according to HardyWeinberg equilibrium. The results are plotted in Figure 1. Small samples result in minor deviations from sample HWE, allowing more than one solution. The smaller the sample size, the greater the range of allele frequencies over which this occurs. A sample of 10 subjects allows more than one biologically possible solution at a wide range of allele frequencies (Figure 1A). With 60 individuals a broad range of allele frequencies is still affected (Figure 1B) – this has implications for analyses based on the HapMap CEU dataset of 60 unrelated individuals [23,24]. At 100 individuals (Figure 1C) the problem is limited to allele frequencies below 15% (Figure 1C), while the plot for 1000 individuals shows no condition under which there is more than one biologically possible solution (Figure 1D). This last observation is because under perfect sample HWE (infinite samples) the number of biologically possible solutions is always 1, despite the number of real solutions exceeding 1 at lower allele frequencies (data not shown).
Figure 1. Simulated data in which HWE is observed to the limit of rounding errors (whole number values for counts of individuals). (A) Number of biologically possible solutions to the cubic equation in (A) 10 individuals; (B) 60 individuals; (C) 100 individuals (D) 1000 individuals. xaxis: allele frequency of SNP1, yaxis: allele frequency of SNP2. Black = more than one solution. Grey = one solution.
Number of solutions to the cubic equation with real data
We have also calculated the number of solutions to equation 1 for a set of real data from the HapMap project [23,24]. These were a selection of SNPs from the ACEGH1 region of chromosome 17 for the CEU population (60 unrelated individuals). Figure 2A shows that at the lower allele frequencies the possibility of more than one real solution to the cubic equation begins to arise. This is consistent with the simulated data for 60 samples (Figure 1B), except that a broader range of allele frequencies is affected. This is probably due to the inherent errors of real data increasing the deviations from HWE relative to nearperfect simulated data. In most cases of multiple solutions only two of the three real roots are biologically possible. Figure 2B compares these two values, indicating that in most cases the differences in estimated haplotype are small. In the minority of cases with three solutions these fit the same pattern. However, this can have major consequences for the calculation of D' (as illustrated in Figure 3). Note that D' and r^{2 }behave quite differently in this respect, and r^{2 }is much less affected. However, as a D' of 1 indicates the existence of three or less haplotypes (r^{2 }of 1 indicates two haplotypes), D' is a good indicator of haplotype block structure, with a value of exactly 1 suggesting little or no recombination between two loci, and a value less than 1 supporting a breakdown of LD. In fact CubeX provides both D' and r^{2}, allowing the user to select their measure of preference. Figure 4 illustrates the relationship between these two measures in the simulated and real datasets, which clarifies how a large D' value can be observed with a low r^{2 }value, but the key point is that a D' of 1 indicates complete LD (i.e. three or less haplotypes) despite a low r^{2}.
Figure 2. Evaluation of number of solutions for real data. (A) Number of biologically possible solutions over a range of allele frequencies using a large sample of SNP data (Chr. 17:60 to 60.5 MB, 121 SNPs) from the HapMap project [23,24]. xaxis: allele frequency of SNP1, yaxis: allele frequency of SNP2. Black = more than one solution. Grey = one solution. (B) Comparison of two solutions within the dataset. xaxis: higher value solution, yaxis: lower value solution.
Figure 3. Screenshot of results screen from CubeX online analysis program. In this example there are two biologically possible solutions. Results for both are shown (upper table), and observed (input values) and expected diplotype frequencies (for the two solutions) displayed for comparison (lower table).
Figure 4. The range of LD in datasets using the CubeX tool to calculate r2 and D'. (A) Simulated data. D' on xaxis, r^{2 }on y axis. (B) Real SNP data (Chr. 17:60 to 60.5 MB, 121 SNPs) from the HapMap project [23,24]. D' on xaxis, r^{2 }on y axis.
Comparison of the cubic exact solution with other approaches
For the purposes of comparison we have analysed two datasets with PHASE [16], MIDAS [7] (Hill EM) and CubeX. The first is a dataset of directly haplotyped samples comprising 80 subjects from 3 ethnic groups (Asian, African and Caucasian) for APOE [25]. Although all but one SNP was in HardyWeinberg equilibrium, this dataset has the potential to invalidate some of the assumptions of the programs due to the mixture of ethnicity. However, this provides a useful substrate on which to test the influence of stratification on the outcome of the cubic exact solution. The second dataset is a set of multilocus phased data from HapMap CEU samples [23,24] for the IGF2 gene region. Although these have not been directly haplotyped, the multilocus phased haplotypes are expected to be very accurate, and this dataset comprises Caucasians, so will not suffer from the same stratification issues. We tested the programs on pairwise subsets of these data.
For the APOE [25] dataset the data are presented in Additional File 1, with a selected summary in Table 1. The subset in Table 1 demonstrate the advantage of being provided with all possible solutions by CubeX, but also demonstrates that all three approaches can be wrong. To summarise the outcome, PHASE [16] and MIDAS [7] (Hill EM) both matched the real counts in 28 of 36 SNP pairs, while CubeX matched real counts in 33 of 36 SNP pairs (for one of its solutions). However, in five of those cases the user would need to determine which of the two CubeX solutions to use based on their prior knowledge of the LD structure in the region (i.e. do they expect three or four haplotypes). This comparison confirms the risk of EM finding a local maximum when there is more than one biologically possible solution, and suggests that CubeX may offer advantages in stratified datasets or datasets with low SNP minor allele frequencies (confirming the results from simulated data above).
Additional file 1. Comparisons of PHASE, MIDAS and CubeX on APOE data (from [25]). A comparison of PHASE, MIDAS and CubeX for pairwise analysis of genotype data derived from directly observed multilocus haplotypes.
Format: PDF Size: 82KB Download file
This file can be viewed with: Adobe Acrobat Reader
Table 1. Illustrative examples of comparison of CubeX with PHASE [16] and MIDAS [7] (Hill EM).
For the HapMap [23,24] IGF2 region data (comprising SNPs rs3802971, rs734351, rs3213221, rs4244808, rs1003483, rs3741208, rs1004446, rs4320932 and rs7924316) CubeX gives only one solution in all cases, and there is little difference between the outcome of the three approaches (Additional File 2). This confirms that in situations of higher allele frequencies there is less of an issue with multiple biologically possible solutions to the cubic equation, and iterative approaches are completely acceptable.
Additional file 2. Comparisons of PHASE, MIDAS and CubeX on HapMap IGF2 region data (from http://www.hapmap.org webcite, [23,24]). A comparison of PHASE, MIDAS and CubeX for pairwise analysis of genotype data derived from statistically inferred longrange multilocus haplotypes.
Format: PDF Size: 81KB Download file
This file can be viewed with: Adobe Acrobat Reader
Discussion
We have written an online program, "CubeX", to enable simple analysis of the biologically
possible estimated haplotypes for pairs of biallelic markers. This program takes data
from a pair of markers as a standard 3 × 3 table of nine diplotypes, generates cubic
exact solutions to equation 1 and generates output in the format shown in Figure 3. The number of possible solutions is shown, followed by haplotype frequencies and
LD statistics for those solutions. Below that a duplicate of the 3 × 3 input table
is displayed with the addition of expected absolute diplotype frequencies calculated
from the haplotype frequencies. The difference between these and the input data are
subjected to a χ^{2 }test, which effectively tests sample deviation from the null hypothesis of HWE for
the diplotypes formed of the four haplotypes. However, the interpretation of solutions
depends on the prior hypothesis. In the example in Figure 3, although solution γ exhibits a slightly worse χ^{2 }fit than solution β, the former is consistent with a prior hypothesis of only three of the four haplotypes
existing (see Figure 5 in reference [7]), which is biologically likely in the absence of recombination between any two loci.
In fact, in all tested cases in Figure 2 generating more than one solution, the diplotype data included zero values in at
least one corner cell and the two adjacent edge cells of the 3 × 3 (i.e. where one
possible solution has a D' = 1, although it should be noted that more than one solution
can occur without zero values if double heterozygotes are greatly overrepresented).
This suggests that the principal issue is whether three or four haplotypes exist,
and in these cases the prior hypothesis (based on distance and recombination rates)
is of utmost importance. If input data for individual SNPs are significantly out of
HWE a warning message is given at the top of the page. For completeness, the biologically
impossible real number solutions are displayed at the bottom, along with minimum and
maximum biologically possible values for
Under perfect sample HWE the frequencies of all haplotypes can be directly inferred
from the corresponding corner diplotypes of the 3 × 3. For example:
Finally, we have demonstrated by comparison with PHASE [16] and MIDAS [7] (Hill EM) that in certain situations (low minor allele frequency, population stratification) the cubic exact approach can perform better for pairwise analyses than alternative approaches by indicating the existence of multiple solutions. However, our findings confirm that in most other situations iterative approaches are robust and accurate.
Conclusion
We present a comprehensive analysis of the consequences of different variables on the number of solutions to the cubic equation for haplotype frequency. Our analyses demonstrate that lower allele frequencies, lower sample numbers and a possible D' value of 1 can result in more than one solution. This has significant implications for the calculation of LD in small sample sizes and with rarer alleles that may have particular disease relevance. This evaluation provides essential information for an understanding of the limitations of LD estimation, which is particularly relevant for genomewide analyses (where sample sizes and allele frequencies can be low). Finally, we present a program "CubeX", freely available as an online program, which provides each of the biologically possible cubic exact solution(s) to equation 1 for haplotype frequency, enabling the user to identify the solution that best fits their prior hypothesis for number of haplotypes.
Availability and Requirements
Project name: CubeX
Project home page: http://www.oege.org/software/cubex webcite
Operating system(s): Platform independent (webbased)
Programming language: Python http://www.python.org webcite
Licence: CubeX licence available from http://www.oege.org/software/cubex webcite
Any restrictions to use by nonacademics: royaltyfree use allowed within terms of licence
Abbreviations
EM – ExpectationMaximisation
HWE – HardyWeinberg Equilibrium
LD – Linkage Disequilibrium
Authors' contributions
TRG wrote the CubeX program, ran the simulations and analyses and drafted the manuscript. SR advised on LD calculation and output format, tested the program and contributed to the manuscript. INMD drafted the solution to the cubic equation, advised on methods, tested the program and contributed to the manuscript. All authors read and approved the final manuscript.
Acknowledgements
TRG is funded by a BHF (British Heart Foundation) Intermediate Fellowship (FS/05/065/19497), SR by a HOPE (Wessex Medical Trust) fellowship and work in our laboratory by the Medical Research Council (UK) (Programme Grant G9800748). We thank an anonymous reviewer for their suggestion of a comparison with PHASE on the APOE dataset [25].
References

Weiss KM, Clark AG: Linkage disequilibrium and the mapping of complex human traits.
Trends in Genetics 2002, 18:1924. PubMed Abstract  Publisher Full Text

Palmer LJ, Cardon LR: Shaking the tree: mapping complex disease genes with linkage disequilibrium.
The Lancet 2005, 366:12231234. Publisher Full Text

Ardlie KG, Kruglyak L, Seielstad M: Patterns of linkage disequilibrium in the human genome.
Nat Rev Genet 2002, 3:299309. PubMed Abstract  Publisher Full Text

Hill WG: Estimation of linkage disequilibrium in randomly mating populations.
Heredity 1974, 33:229239. PubMed Abstract

Abecasis GR, Cookson WO: GOLDgraphical overview of linkage disequilibrium.
Bioinformatics 2000, 16:182183. PubMed Abstract  Publisher Full Text

Pettersson F, Jonsson O, Cardon LR: GOLDsurfer: three dimensional display of linkage disequilibrium.
Bioinformatics 2004, 20:32413243. PubMed Abstract  Publisher Full Text

Gaunt TR, Rodriguez S, Zapata C, Day IN: MIDAS: software for analysis and visualisation of interallelic disequilibrium between multiallelic markers.
BMC Bioinformatics 2006, 7:227. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Barrett JC, Fry B, Maller J, Daly MJ: Haploview: analysis and visualization of LD and haplotype maps.
Bioinformatics 2005, 21:263265. PubMed Abstract  Publisher Full Text

Jorde LB: Linkage disequilibrium and the search for complex disease genes.
Genome Res 2000, 10:14351444. PubMed Abstract  Publisher Full Text

Mueller JC: Linkage disequilibrium for different scales and applications.
Brief Bioinform 2004, 5:355364. PubMed Abstract  Publisher Full Text

Weale ME: A survey of current software for haplotype phase inference.
Hum Genomics 2004, 1:141144. PubMed Abstract

Salem RM, Wessel J, Schork NJ: A comprehensive literature review of haplotyping software and methods for use with unrelated individuals.
Human Genomics 2005, 2:3966. PubMed Abstract

Wang L, Xu Y: Haplotype inference by maximum parsimony.
Bioinformatics 2003, 19:17731780. PubMed Abstract  Publisher Full Text

Stephens M, Scheet P: Accounting for decay of linkage disequilibrium in haplotype inference and missingdata imputation.
Am J Hum Genet 2005, 76:449462. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stephens M, Donnelly P: A comparison of bayesian methods for haplotype reconstruction from population genotype data.
Am J Hum Genet 2003, 73:11621169. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data.
Am J Hum Genet 2001, 68:978989. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nickalls RWD: A new approach to solving the cubic: Cardan's solution revealed.
The Mathematical Gazette 1993, 77:354359. Publisher Full Text

Luo ZW, Suhai S: Estimating Linkage Disequilibrium Between a Polymorphic Marker Locus and a Trait Locus in Natural Populations.
Genetics 1999, 151:359371. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Foundation PS: The Python Programming Language. [http://www.python.org] webcite
2006.

Lewontin RC: The interaction of selection and linkage. I. General considerations; heterotic models.
Genetics 1964, 49:4967. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hill WG, Robertson A: Linkage disequilibrium in finite populations.

Mano S, Yasuda N, Katoh T, Tounai K, Inoko H, Imanishi T, Tamiya G, Gojobori T: Notes on the Maximum Likelihood Estimation of Haplotype Frequencies.
Annals of Human Genetics 2004, 68:257264. PubMed Abstract  Publisher Full Text

The International HapMap Project
Nature 2003, 426:789796. PubMed Abstract  Publisher Full Text

Consortium TIHM: A haplotype map of the human genome.
Nature 2005, 437:12991320. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Orzack SH, Gusfield D, Olson J, Nesbitt S, Subrahmanyan L, Stanton VP Jr.: Analysis and Exploration of the Use of RuleBased Algorithms and Consensus Methods for the Inferral of Haplotypes.
Genetics 2003, 165:915928. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Day INM, Alharbi KK, Smith MJ, Aldahmesh MA, Chen X, Lotery AJ, PantedeSousa G, Hou G, Ye S, Eccles DM, Cross NCP, Fox KR, Rodriguez S: Paucimorphic Alleles versus Polymorphic Alleles and Rare Mutations in Disease Causation: Theory, Observation and Detection.
Current Genomics 2004, 5:431438. Publisher Full Text

Alharbi KK, Aldahmesh MA, Spanakis E, Haddad L, Whittall RA, Chen X, Rassoulian H, Smith MJ, Sillibourne J, Ball NJ, Graham NJ, Briggs PJ, Simpson IA, Phillips DIW, Lawlor DA, Ye S, Humphries SE, Cooper C, Smith GD, Ebrahim S, Eccles DM, Day INM: Mutation scanning by meltMADGE: Validations using BRCA1 and LDLR, and demonstration of the potential to identify severe, moderate, silent, rare, and paucimorphic mutations in the general population.
Genome Res 2005, 15:967977. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Alharbi KK, Spanakis E, Tan K, Smith MJ, Aldahmesh MA, O'Dell SD, Sayer AA, Lawlor DA, Ebrahim S, Davey Smith G, O'Rahilly S, Farooqi S, Cooper C, Phillips DI, Day IN: Prevalence and functionality of paucimorphic and private MC4R mutations in a large, unselected European British population, scanned by meltMADGE.
Hum Mutat 2007, 28(3):294302. PubMed Abstract  Publisher Full Text

Wong GKS, Yang Z, Passey DA, Kibukawa M, Paddock M, Liu CR, Bolund L, Yu J: A Population Threshold for Functional Polymorphisms.
Genome Res 2003, 13:18731879. PubMed Abstract  Publisher Full Text  PubMed Central Full Text