BACKGROUND
linkage disequilibrium  describes the condition that occurs when alleles at different loci are non-randomly associated in a given population.
under ld the frequency  of a haplotype  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  <cit> .
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  <cit> .
for unrelated individuals the estimation of ld relies on the estimation of haplotype frequencies.
in a  <dig> ×  <dig> table for a biallelic marker the haplotype phase of all individuals is known with the exception of the centre cell .
the estimated frequency, f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@, of the haplotype h <dig> is described by a cubic equation of the form
af^113+bf^112+cf^11+d=0
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaemyyaemafmozaymbakaadaqhaawcbagaegymaejaegymaedabagaeg4mamdaaogaey4kasiaemoyaimafmozaymbakaadaqhaawcbagaegymaejaegymaedabagaegomaidaaogaey4kasiaem4yammafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabgucariabdsgakjabg2da9iabicdawaaa@424c@
that is adapted from hill's equation   <cit>  with the constants defined under methods.
with f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and the allele frequencies, all four haplotype frequencies can be calculated, thus estimating the unknown proportions of the middle cell.
several approaches exist for solving equation , the solution of which enables estimation of haplotype frequencies and ld coefficients.
the first approach uses iteration-based algorithms.
an initial estimate of haplotype frequency  is substituted into the equation, providing a new estimate.
this is then fed back into the equation and the expectation-maximisation  process repeated until the values converge.
this is the basis both of the algorithm described by hill in  <dig> for the estimation of pairwise haplotype frequencies  <cit> , and of other em algorithms that enable the estimation of multilocus haplotype frequencies.
many programs exist that utilise variations on this approach, including: gold  <cit> , goldsurfer  <cit> , midas  <cit> , haploview  <cit>  and many others reviewed in  <cit> .
the potential problem for these approaches is that algorithms may converge on one of the alternative roots of the cubic equation .
other approaches include parsimony, eg hapar  <cit>  and bayesian algorithms, eg phase  <cit> .
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  <cit> .
an alternative approach would be exact solution, such as cardan's solution  <cit>  of the generalized cubic equation  is an example).
this provides all roots to the cubic equation, from which we can select those that are both real  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 non-iterative 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   <cit> .
rearranging terms for consequent diplotype frequency expectations for two biallelic loci luo and suhai  <cit>  obtained equation  <dig> given in the introduction : ax <dig> + bx <dig> + cx + d =  <dig>  where a = 4n; b = 2n  -  <dig> - n22; c = 2npq -  - n22; d = -pq; n = number of subjects; p = common allele freq of locus 1; q = common allele freq of locus 2; n <dig> is the number of subjects who are homozygous for the commoner allele at both loci; n <dig> are common homozygous at locus  <dig> and heterozygous at locus 2; n <dig> are heterozygous at locus  <dig> and common homozygous at locus 2; n <dig> are heterozygous at both loci  <cit> .
equation  <dig> can be solved exactly for x .
we have adopted the nickalls treatment of the cardan solution of the generalized cubic equation  <cit> , and written a python  <cit>  program "cubex" to solve equation  <dig> exactly.
in cubex, after calculation of constants a-d from diplotypic data the following are calculated:
xn = -b/; δ <dig> = /9a2; h <dig> = 4a2δ6; yn = axn <dig> + bxn <dig> + cxn + d.
the discriminant Δ <dig> = yn <dig> - h <dig> is then used to determine the outcome in real roots , with three possible outcomes:
outcome 1: if yn <dig> > h <dig> there will be only one real root  given by
α=xn+12a3+12a3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8xsdemaeyypa0jaemieag3aasbaasqaaiabd6eaobqabagccqghrawkdagcbaqaakqbaoaalaaabagaegymaedabagaegomaijaemyyaegaaowaaewaaeaacqghsislcqwg5bqedawgaawcbagaemota4eabeaakiabgucarmaakaaabagaemyeak3aa0baasqaaiabd6eaobqaaiabikdayaaakiabgkhitiabdigaonaacaaaleqabagaegomaidaaaqabaaakiaawicacaglpaaaasqaaiabiodazaaakiabgucarmaakeaabaqcfa4aasaaaeaacqaixaqmaeaacqaiyagmcqwghbqyaagcdaqadaqaaiabgkhitiabdmha5naabaaaleaacqwgobgtaeqaaogaeyoei0yaaoaaaeaacqwg5bqedaqhaawcbagaemota4eabagaegomaidaaogaeyoei0iaemiaag2aawbaasqabeaacqaiyagmaaaabeaaaogaayjkaiaawmcaaawcbagaeg4mamdaaaaa@582e@
outcome 2: if yn <dig> = h <dig> there are three real roots  and α and β are equal.
for a value of μ=yn2a3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8hvd0maeyypa0zaaoqaaeaajuagdawcaaqaaiabdmha5naabaaabagaemota4eabeaaaeaacqaiyagmcqwghbqyaaaaleaacqaizawmaaaaaa@3541@:
α = xn + μ
β = xn + μ
γ = xn - 2μ
outcome 3: if yn <dig> <h <dig> there are three real roots .
where θ=arccos⁡3
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8hudenaeyypa0tcfa4aasaaaeaacygghbqycqggybgccqggjbwycqggjbwycqggvbwbcqggzbwccqggoaakcqghsislcqwg5bqedawgaaqaaiabd6eaobqabagaei4la8iaemiaagmaeiykakcabagaeg4mamdaaaaa@3fef@:
α = xn + 2δcosθ
β = xn + 2δcos
γ = xn + 2δcos
values for d' and r <dig> are calculated as previously described  <cit> :
d=−
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaeeiraqkaeyypa0jaeiikagiafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabgena0kqbdagamzaajawaasbaasqaaiabikdayiabikdayaqabagccqggpaqkcqghsislcqggoaakcuwgmbgzgaqcamaabaaaleaacqaixaqmcqaiyagmaeqaaogaey41aqrafmozaymbakaadawgaawcbagaegomaijaegymaedabeaakiabcmcapaaa@44a9@
dmax = min  if d >  <dig> or dmax = min  if d < 0
d' = d/dmax
r <dig> = d2/q)
diplotype frequencies based on the estimated haplotype frequencies are compared to the input diplotype frequencies by a χ <dig> 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  minus four estimated parameters which are either three haplotypes  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  <dig> ×  <dig> the number of degrees of freedom is reduced by one.
if the user knows there are only three haplotypes present  then there are only three estimated parameters  and  <dig> 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 non-zero 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 f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and the derived f^12
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegomaidabeaaaaa@2f46@, f^21
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegomaijaegymaedabeaaaaa@2f46@ and f^22
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegomaijaegomaidabeaaaaa@2f48@ all fall within the range  <dig> to  <dig>  and add up to  <dig>  this constraint is tighter than those described elsewhere  <cit>  as it relies on the inherent assumption of representative sampling and hwe, an extreme chance distortion of which could lead to three solutions at snp allele frequencies of  <dig>  in sample data drawn from a population .
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  <dig> and  <dig>  at intervals of 1/ <dig>  that add up to  <dig>  these haplotype frequencies were then converted to diplotype frequencies according to hardy-weinberg equilibrium.
the results are plotted in figure  <dig>  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  <dig> subjects allows more than one biologically possible solution at a wide range of allele frequencies .
with  <dig> individuals a broad range of allele frequencies is still affected  – this has implications for analyses based on the hapmap ceu dataset of  <dig> unrelated individuals  <cit> .
at  <dig> individuals  the problem is limited to allele frequencies below 15% , while the plot for  <dig> individuals shows no condition under which there is more than one biologically possible solution .
this last observation is because under perfect sample hwe  the number of biologically possible solutions is always  <dig>  despite the number of real solutions exceeding  <dig> at lower allele frequencies .
number of solutions to the cubic equation with real data
we have also calculated the number of solutions to equation  <dig> for a set of real data from the hapmap project  <cit> .
these were a selection of snps from the ace-gh <dig> region of chromosome  <dig> for the ceu population .
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  <dig> samples , 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 near-perfect 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' .
note that d' and r <dig> behave quite differently in this respect, and r <dig> is much less affected.
however, as a |d'| of  <dig> indicates the existence of three or less haplotypes , |d'| is a good indicator of haplotype block structure, with a value of exactly  <dig> suggesting little or no recombination between two loci, and a value less than  <dig> supporting a break-down of ld.
in fact cubex provides both d' and r <dig>  allowing the user to select their measure of preference.
figure  <dig> 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 <dig> value, but the key point is that a |d'| of  <dig> indicates complete ld  despite a low r <dig>
comparison of the cubic exact solution with other approaches
for the purposes of comparison we have analysed two datasets with phase  <cit> , midas  <cit>   and cubex.
the first is a dataset of directly haplotyped samples comprising  <dig> subjects from  <dig> ethnic groups  for apoe  <cit> .
although all but one snp was in hardy-weinberg 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 multi-locus phased data from hapmap ceu samples  <cit>  for the igf <dig> gene region.
although these have not been directly haplotyped, the multi-locus 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 pair-wise subsets of these data.
for the apoe  <cit>  dataset the data are presented in additional file  <dig>  with a selected summary in table  <dig>  the subset in table  <dig> 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  <cit>  and midas  <cit>   both matched the real counts in  <dig> of  <dig> snp pairs, while cubex matched real counts in  <dig> of  <dig> snp pairs .
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 .
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 .
a selection of comparisons using direct haplotyped apoe data  <cit> .
full data are present as a additional table.
for haplotype numbers  incorrect answers are marked '*', correct answers are unmarked.
examples:  phase, midas and cubex  give correct answer.
only cubex gives the correct answer as one of its two solutions.
midas and cubex give the correct answer, phase and the other cubex solution are wrong.
all three approaches are wrong.
phase and cubex give the correct answer, midas and the other cubex solution are wrong.
for the hapmap  <cit>  igf <dig> region data  cubex gives only one solution in all cases, and there is little difference between the outcome of the three approaches .
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.
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  <dig> ×  <dig> table of nine diplotypes, generates cubic exact solutions to equation  <dig> and generates output in the format shown in figure  <dig>  the number of possible solutions is shown, followed by haplotype frequencies and ld statistics for those solutions.
below that a duplicate of the  <dig> ×  <dig> 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 χ <dig> 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  <dig>  although solution γ exhibits a slightly worse χ <dig> fit than solution β, the former is consistent with a prior hypothesis of only three of the four haplotypes existing , which is biologically likely in the absence of recombination between any two loci.
in fact, in all tested cases in figure  <dig> 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  <dig> ×  <dig> .
this suggests that the principal issue is whether three or four haplotypes exist, and in these cases the prior hypothesis  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 f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and allele frequencies.
this program provides a convenient utility for researchers to both analyse data for haplotype frequencies and ld statistics and to check previous analyses for potential problems caused by multiple solutions.
under perfect sample hwe the frequencies of all haplotypes can be directly inferred from the corresponding corner diplotypes of the  <dig> ×  <dig>  for example: n11=n f^112
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaemoba42aasbaasqaaiabigdaxiabigdaxaqabagccqgh9aqpcqqgubgbcqqggaaicuwgmbgzgaqcamaadaaaleaacqaixaqmcqaixaqmaeaacqaiyagmaaaaaa@36e2@, so f^11=n11n
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabg2da9maakaaabaqcfa4aasaaaeaacqwgubgbdawgaaqaaiabigdaxiabigdaxaqabaaabagaemoba4gaaawcbeaaaaa@35d8@.
that being the case there are only two possible values for f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@, one positive and one negative, the latter being biologically impossible.
perfect sample hwe therefore results in only a single biologically possible solution to the cubic equation.
in the case of extreme sample hwd where all samples fall within the middle cell of the  <dig> ×  <dig>  f^11
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ can contribute either a half, a quarter or none of the haplotypes to the middle cell.
there are therefore three biologically possible solutions under conditions of extreme sample hwd.
the results from real data confirm that in some cases more than one biologically possible solution to the cubic equation for haplotype frequency can exist.
the simulations suggest that this occurs where small sample size, sampling errors or non-random mating result in a distortion of sample hwe, and demonstrates the importance of testing hwe before haplotype analyses.
the greater the distortion of sample hwe the higher the allele frequency at which more than one solution can occur .
in these cases the cubic exact algorithm gives all possible solutions and a test of hwe, while an iteration-based method would only give one.
this supports the hypothesis that the cubic exact approach is superior to iteration-based methods in real-world datasets where sample data rarely fit exactly to hwe .
this is particularly important in the analysis of low frequency snps and paucimorphisms  <cit> , for which different solutions can significantly distort d' results, despite the relatively similar solutions giving similar r <dig> results.
in all the observed data with two solutions there were no occasions in which r <dig> exceeded  <dig>  for any biologically possible solution, and in most cases there is only a small difference in r <dig> between biologically possible solutions.
the largest effect is on d'.
on the basis of empirical data and using different approaches to inference wong et al showed that coding snps with minor allele frequencies < <dig>  are likely to be of functional importance  <cit> , and rarer alleles, haplotypes and diplotypes of causal importance have emerged in numerous disease contexts .
in addition to being applicable and giving exact evaluation for d' analysis of common snps, the cubic exact solution may prove of particular value for evaluating "post-hapmap" and "post-dbsnp" rarer haplotypes, for fully evaluating d' estimates from datasets with greater deviations from the random mating and hwe assumptions and for fully evaluating ld in small datasets.
finally, we have demonstrated by comparison with phase  <cit>  and midas  <cit>   that in certain situations  the cubic exact approach can perform better for pair-wise 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.
CONCLUSIONS
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  <dig> 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 genome-wide analyses .
finally, we present a program "cubex", freely available as an online program, which provides each of the biologically possible cubic exact solution to equation  <dig> 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:
operating system: platform independent
programming language: python
licence: cubex licence available from
any restrictions to use by non-academics: royalty-free use allowed within terms of licence
abbreviations
em – expectation-maximisation
hwe – hardy-weinberg 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.
supplementary material
additional file 1
comparisons of phase, midas and cubex on apoe data .
a comparison of phase, midas and cubex for pairwise analysis of genotype data derived from directly observed multi-locus haplotypes.
click here for file
additional file 2
comparisons of phase, midas and cubex on hapmap igf <dig> region data .
a comparison of phase, midas and cubex for pairwise analysis of genotype data derived from statistically inferred long-range multi-locus haplotypes.
click here for file
acknowledgements
trg is funded by a bhf  intermediate fellowship , sr by a hope  fellowship and work in our laboratory by the medical research council  .
we thank an anonymous reviewer for their suggestion of a comparison with phase on the apoe dataset  <cit> .

 linkage disequilibrium describes the condition that occurs when alleles at different loci are non-randomly associated in a given population. under ld the frequency of a haplotype 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. 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. for unrelated individuals the estimation of ld relies on the estimation of haplotype frequencies. in a × table for a biallelic marker the haplotype phase of all individuals is known with the exception of the centre cell. the estimated frequency, f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@, of the haplotype h is described by a cubic equation of the form af^113+bf^112+cf^11+d=0 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaemyyaemafmozaymbakaadaqhaawcbagaegymaejaegymaedabagaeg4mamdaaogaey4kasiaemoyaimafmozaymbakaadaqhaawcbagaegymaejaegymaedabagaegomaidaaogaey4kasiaem4yammafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabgucariabdsgakjabg2da9iabicdawaaa@424c@ that is adapted from hill's equation with the constants defined under methods. with f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and the allele frequencies, all four haplotype frequencies can be calculated, thus estimating the unknown proportions of the middle cell. several approaches exist for solving equation, the solution of which enables estimation of haplotype frequencies and ld coefficients. the first approach uses iteration-based algorithms. an initial estimate of haplotype frequency is substituted into the equation, providing a new estimate. this is then fed back into the equation and the expectation-maximisation process repeated until the values converge. this is the basis both of the algorithm described by hill in for the estimation of pairwise haplotype frequencies, and of other em algorithms that enable the estimation of multilocus haplotype frequencies. many programs exist that utilise variations on this approach, including: gold, goldsurfer, midas, haploview and many others reviewed in. the potential problem for these approaches is that algorithms may converge on one of the alternative roots of the cubic equation. other approaches include parsimony, eg hapar and bayesian algorithms, eg phase. 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. an alternative approach would be exact solution, such as cardan's solution of the generalized cubic equation is an example). this provides all roots to the cubic equation, from which we can select those that are both real 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 non-iterative 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. rearranging terms for consequent diplotype frequency expectations for two biallelic loci luo and suhai obtained equation given in the introduction: ax + bx + cx + d = where a = 4n; b = 2n - - n22; c = 2npq - - n22; d = -pq; n = number of subjects; p = common allele freq of locus 1; q = common allele freq of locus 2; n is the number of subjects who are homozygous for the commoner allele at both loci; n are common homozygous at locus and heterozygous at locus 2; n are heterozygous at locus and common homozygous at locus 2; n are heterozygous at both loci. equation can be solved exactly for x. we have adopted the nickalls treatment of the cardan solution of the generalized cubic equation, and written a python program "cubex" to solve equation exactly. in cubex, after calculation of constants a-d from diplotypic data the following are calculated: xn = -b/; δ = /9a2; h = 4a2δ6; yn = axn + bxn + cxn + d. the discriminant Δ = yn - h is then used to determine the outcome in real roots, with three possible outcomes: outcome 1: if yn > h there will be only one real root given by α=xn+12a3+12a3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8xsdemaeyypa0jaemieag3aasbaasqaaiabd6eaobqabagccqghrawkdagcbaqaakqbaoaalaaabagaegymaedabagaegomaijaemyyaegaaowaaewaaeaacqghsislcqwg5bqedawgaawcbagaemota4eabeaakiabgucarmaakaaabagaemyeak3aa0baasqaaiabd6eaobqaaiabikdayaaakiabgkhitiabdigaonaacaaaleqabagaegomaidaaaqabaaakiaawicacaglpaaaasqaaiabiodazaaakiabgucarmaakeaabaqcfa4aasaaaeaacqaixaqmaeaacqaiyagmcqwghbqyaagcdaqadaqaaiabgkhitiabdmha5naabaaaleaacqwgobgtaeqaaogaeyoei0yaaoaaaeaacqwg5bqedaqhaawcbagaemota4eabagaegomaidaaogaeyoei0iaemiaag2aawbaasqabeaacqaiyagmaaaabeaaaogaayjkaiaawmcaaawcbagaeg4mamdaaaaa@582e@ outcome 2: if yn = h there are three real roots and α and β are equal. for a value of μ=yn2a3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8hvd0maeyypa0zaaoqaaeaajuagdawcaaqaaiabdmha5naabaaabagaemota4eabeaaaeaacqaiyagmcqwghbqyaaaaleaacqaizawmaaaaaa@3541@: α = xn + μ β = xn + μ γ = xn - 2μ outcome 3: if yn <h there are three real roots. where θ=arccos⁡3 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbaaccigae8hudenaeyypa0tcfa4aasaaaeaacygghbqycqggybgccqggjbwycqggjbwycqggvbwbcqggzbwccqggoaakcqghsislcqwg5bqedawgaaqaaiabd6eaobqabagaei4la8iaemiaagmaeiykakcabagaeg4mamdaaaaa@3fef@: α = xn + 2δcosθ β = xn + 2δcos γ = xn + 2δcos values for d' and r are calculated as previously described: d=− mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xi8qivkypfjydhavhbbf9v8qqaqfr0xc9vqfj0dxdbba91qpepei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaeeiraqkaeyypa0jaeiikagiafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabgena0kqbdagamzaajawaasbaasqaaiabikdayiabikdayaqabagccqggpaqkcqghsislcqggoaakcuwgmbgzgaqcamaabaaaleaacqaixaqmcqaiyagmaeqaaogaey41aqrafmozaymbakaadawgaawcbagaegomaijaegymaedabeaakiabcmcapaaa@44a9@ dmax = min if d > or dmax = min if d < 0 d' = d/dmax r = d2/q) diplotype frequencies based on the estimated haplotype frequencies are compared to the input diplotype frequencies by a χ 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 minus four estimated parameters which are either three haplotypes 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 × the number of degrees of freedom is reduced by one. if the user knows there are only three haplotypes present then there are only three estimated parameters and 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 non-zero 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.  solutions are considered biologically possible when f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and the derived f^12 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegomaidabeaaaaa@2f46@, f^21 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegomaijaegymaedabeaaaaa@2f46@ and f^22 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegomaijaegomaidabeaaaaa@2f48@ all fall within the range to and add up to this constraint is tighter than those described elsewhere as it relies on the inherent assumption of representative sampling and hwe, an extreme chance distortion of which could lead to three solutions at snp allele frequencies of in sample data drawn from a population. 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 and at intervals of 1/ that add up to these haplotype frequencies were then converted to diplotype frequencies according to hardy-weinberg equilibrium. the results are plotted in figure 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 subjects allows more than one biologically possible solution at a wide range of allele frequencies. with individuals a broad range of allele frequencies is still affected – this has implications for analyses based on the hapmap ceu dataset of unrelated individuals. at individuals the problem is limited to allele frequencies below 15%, while the plot for individuals shows no condition under which there is more than one biologically possible solution. this last observation is because under perfect sample hwe the number of biologically possible solutions is always despite the number of real solutions exceeding at lower allele frequencies. number of solutions to the cubic equation with real data we have also calculated the number of solutions to equation for a set of real data from the hapmap project. these were a selection of snps from the ace-gh region of chromosome for the ceu population. 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 samples, 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 near-perfect 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'. note that d' and r behave quite differently in this respect, and r is much less affected. however, as a |d'| of indicates the existence of three or less haplotypes, |d'| is a good indicator of haplotype block structure, with a value of exactly suggesting little or no recombination between two loci, and a value less than supporting a break-down of ld. in fact cubex provides both d' and r allowing the user to select their measure of preference. figure 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 value, but the key point is that a |d'| of indicates complete ld despite a low r comparison of the cubic exact solution with other approaches for the purposes of comparison we have analysed two datasets with phase, midas and cubex. the first is a dataset of directly haplotyped samples comprising subjects from ethnic groups for apoe. although all but one snp was in hardy-weinberg 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 multi-locus phased data from hapmap ceu samples for the igf gene region. although these have not been directly haplotyped, the multi-locus 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 pair-wise subsets of these data. for the apoe dataset the data are presented in additional file with a selected summary in table the subset in table 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 and midas both matched the real counts in of snp pairs, while cubex matched real counts in of snp pairs. 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. 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. a selection of comparisons using direct haplotyped apoe data. full data are present as a additional table. for haplotype numbers incorrect answers are marked'*', correct answers are unmarked. examples: phase, midas and cubex give correct answer. only cubex gives the correct answer as one of its two solutions. midas and cubex give the correct answer, phase and the other cubex solution are wrong. all three approaches are wrong. phase and cubex give the correct answer, midas and the other cubex solution are wrong. for the hapmap igf region data cubex gives only one solution in all cases, and there is little difference between the outcome of the three approaches. 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. 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 × table of nine diplotypes, generates cubic exact solutions to equation and generates output in the format shown in figure the number of possible solutions is shown, followed by haplotype frequencies and ld statistics for those solutions. below that a duplicate of the × 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 χ 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 although solution γ exhibits a slightly worse χ fit than solution β, the former is consistent with a prior hypothesis of only three of the four haplotypes existing, which is biologically likely in the absence of recombination between any two loci. in fact, in all tested cases in figure 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 ×. this suggests that the principal issue is whether three or four haplotypes exist, and in these cases the prior hypothesis 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 f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ and allele frequencies. this program provides a convenient utility for researchers to both analyse data for haplotype frequencies and ld statistics and to check previous analyses for potential problems caused by multiple solutions. under perfect sample hwe the frequencies of all haplotypes can be directly inferred from the corresponding corner diplotypes of the × for example: n11=n f^112 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagaemoba42aasbaasqaaiabigdaxiabigdaxaqabagccqgh9aqpcqqgubgbcqqggaaicuwgmbgzgaqcamaadaaaleaacqaixaqmcqaixaqmaeaacqaiyagmaaaaaa@36e2@, so f^11=n11n mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaakiabg2da9maakaaabaqcfa4aasaaaeaacqwgubgbdawgaaqaaiabigdaxiabigdaxaqabaaabagaemoba4gaaawcbeaaaaa@35d8@. that being the case there are only two possible values for f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@, one positive and one negative, the latter being biologically impossible. perfect sample hwe therefore results in only a single biologically possible solution to the cubic equation. in the case of extreme sample hwd where all samples fall within the middle cell of the × f^11 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaacpc6xni=xh8vivgi8gi=heeeu0xxdbba9frfj0xb9qqpg0dxdb9aspei8k8fii+fsy=rqgqvepae9pg0db9vqaivgfr0xfr=xfr=xc9adbaqaaegacagaaiaabeqaaeqabiwaaagcbagafmozaymbakaadawgaawcbagaegymaejaegymaedabeaaaaa@2f44@ can contribute either a half, a quarter or none of the haplotypes to the middle cell. there are therefore three biologically possible solutions under conditions of extreme sample hwd. the results from real data confirm that in some cases more than one biologically possible solution to the cubic equation for haplotype frequency can exist. the simulations suggest that this occurs where small sample size, sampling errors or non-random mating result in a distortion of sample hwe, and demonstrates the importance of testing hwe before haplotype analyses. the greater the distortion of sample hwe the higher the allele frequency at which more than one solution can occur. in these cases the cubic exact algorithm gives all possible solutions and a test of hwe, while an iteration-based method would only give one. this supports the hypothesis that the cubic exact approach is superior to iteration-based methods in real-world datasets where sample data rarely fit exactly to hwe. this is particularly important in the analysis of low frequency snps and paucimorphisms, for which different solutions can significantly distort d' results, despite the relatively similar solutions giving similar r results. in all the observed data with two solutions there were no occasions in which r exceeded for any biologically possible solution, and in most cases there is only a small difference in r between biologically possible solutions. the largest effect is on d'. on the basis of empirical data and using different approaches to inference wong et al showed that coding snps with minor allele frequencies < are likely to be of functional importance, and rarer alleles, haplotypes and diplotypes of causal importance have emerged in numerous disease contexts. in addition to being applicable and giving exact evaluation for d' analysis of common snps, the cubic exact solution may prove of particular value for evaluating "post-hapmap" and "post-dbsnp" rarer haplotypes, for fully evaluating d' estimates from datasets with greater deviations from the random mating and hwe assumptions and for fully evaluating ld in small datasets. finally, we have demonstrated by comparison with phase and midas that in certain situations the cubic exact approach can perform better for pair-wise 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.  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 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 genome-wide analyses. finally, we present a program "cubex", freely available as an online program, which provides each of the biologically possible cubic exact solution to equation 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: operating system: platform independent programming language: python licence: cubex licence available from any restrictions to use by non-academics: royalty-free use allowed within terms of licence abbreviations em – expectation-maximisation hwe – hardy-weinberg 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. supplementary material additional file 1 comparisons of phase, midas and cubex on apoe data. a comparison of phase, midas and cubex for pairwise analysis of genotype data derived from directly observed multi-locus haplotypes. click here for file additional file 2 comparisons of phase, midas and cubex on hapmap igf region data. a comparison of phase, midas and cubex for pairwise analysis of genotype data derived from statistically inferred long-range multi-locus haplotypes. click here for file acknowledgements trg is funded by a bhf intermediate fellowship, sr by a hope fellowship and work in our laboratory by the medical research council. we thank an anonymous reviewer for their suggestion of a comparison with phase on the apoe dataset.
the frequency of a haplotype comprising one allele at each of two loci can be expressed as a cubic equation, the solution of which gives that frequency. we present a program, "cubex", which calculates the biologically possible exact solution and provides estimated haplotype frequencies, d', r and χ values for each. we also present an analysis of simulations and real data using the algebraically exact solution, which indicates that under perfect sample hardy-weinberg equilibrium there is only one biologically possible solution, but that under other conditions there may be more. most haplotype and linkage disequilibrium analysis programs use iteration-based 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. 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 hardy-weinberg equilibrium, potentially causing more than one biologically possible solution. our analyses demonstrate that lower allele frequencies, lower sample numbers, population stratification and a possible |d'| value of are particularly susceptible to distortion of sample hardy-weinberg equilibrium, which has significant implications for calculation of linkage disequilibrium in small sample sizes and rarer alleles that may have particular disease relevance and require improved approaches for meaningful evaluation.
 the frequency of a haplotype comprising one allele at each of two loci can be expressed as a cubic equation, the solution of which gives that frequency. most haplotype and linkage disequilibrium analysis programs use iteration-based 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.  we present a program, "cubex", which calculates the biologically possible exact solution and provides estimated haplotype frequencies, d', r and χ 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 hardy-weinberg 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 hardy-weinberg equilibrium there is only one biologically possible solution, but that under other conditions there may be more.  our analyses demonstrate that lower allele frequencies, lower sample numbers, population stratification and a possible |d'| value of are particularly susceptible to distortion of sample hardy-weinberg equilibrium, which has significant implications for calculation of linkage disequilibrium in small sample sizes and rarer alleles that may have particular disease relevance and require improved approaches for meaningful evaluation. 