Optimal Step Length EM Algorithm (OSLEM) for the estimation of haplotype frequency and its application in lipoprotein lipase genotyping

Background Haplotype based linkage disequilibrium (LD) mapping has become a powerful and cost-effective method for performing genetic association studies, particularly in the search for genetic markers in linkage disequilibrium with complex disease loci. Various methods (e.g. Monte-Carlo (Gibbs sampling); EM (expectation maximization); and Clark's method) have been used to estimate haplotype frequencies from routine genotyping data. Results These algorithms can be very slow for large number of SNPs. In order to speed them up, we have developed a new algorithm using numerical analysis technology, a so-called optimal step length EM (OSLEM) that accelerates the calculation. By optimizing approximately the step length of the EM algorithm, OSLEM can run at about twice the speed of EM. This algorithm has been used for lipoprotein lipase (LPL) genotyping analysis. Conclusions This new optimal step length EM (OSLEM) algorithm can accelerate the calculation for haplotype frequency estimation for genotyping data without pedigree information. An OSLEM on-line server is available, as well as a free downloadable version.


Background
Estimation of haplotype frequencies from routine genotyping data plays an important role in LD analysis, and can be achieved by varied methods, including Monte-Carlo (Gibbs sampling [1], PHASE method [2]), EM [3][4][5] and Subtraction [6]. At least two studies have tested and compared some of these programs. Xu and his collaborators [7] have empirically evaluated and compared the accuracy of the Subtraction method [6], the expectation-maximization (EM) method, and the PHASE method [2] for estimating haplotype frequency and for predicting haplotype phase. Summarizing from the studies of Xu and his collaborators in [7]: "Where there was near complete linkage disequilibrium (LD) between SNPs (the NAT2 gene), all three methods provided effective and accurate estimates for haplotype frequencies and individual haplotype phases. For a genomic region in which marked LD was not maintained (the chromosome X locus), the computational methods were adequate in estimating overall haplotype frequencies. However, none of the methods were accurate in predicting individual haplotype phases. The EM and the PHASE methods provided better estimates for overall haplotype frequencies relative to the Subtraction method for both genomic regions." The PHASE algorithm is extremely slow. A comparison of run-times was reported for five SNPs from the NAT2 gene [7]: Subtraction method, 0.01 second; EM method, 13.48 seconds; and PHASE method, over 128 minutes. Zhang and his collaborators [8] pointed out that, " The PHASE method did not yield significantly different results from a simple maximumlikelihood procedure." From the two comparative studies [7,8] and from one simulation study [9], it was shown convincingly that the expectation maximization (EM) algorithm is accurate for estimation of haplotype frequencies.
The EM algorithm is much faster than Monte-Carlo based algorithms. Whereas it is fast with single runs for a relative small number of SNPs, it can be slow with multiple runs and for large number of SNPs. The EM algorithm is an optimization algorithm. In order to obtain a global maximization, EM should run many times from varied starting points. This can be very time consuming. Although in the standard EM algorithm procedure, the step length is the optimal length under the conditional expectation, the step length is not optimal in general. In order to run faster and more accurately, we have developed a new algorithm using numerical analysis technology, a so-called optimal step length EM (OSLEM) to accelerate the calculation. By optimizing approximately the step length of the EM (expectation maximization) algorithm, OSLEM can run at about twice the speed of EM.

Algorithm
We start with the same premises and notations as Stephen et al [2]. Given n diploid individuals from a population, let G = (G 1 ,..., G n ) denote the (known) genotypes for the individuals, let H = (H 1 ,..., H n ) denote the (unknown) corresponding haplotype pairs, let F = (F 1 ,...,F M ) denote the set of (unknown) population haplotype frequencies (the M possible haplotypes are arbitrarily labeled 1,...,M).
Here H i , a random variable depends on F. Let i be the set of all (ordered) haplotype pairs consistent with the multilocus genotype G i , and suppose the distribution of H i on i will follow the Hardy-Weinberg equilibrium.
The EM and OSLEM algorithms attempt to find the haplotype F that maximizes the likelihood.
Here, where i is the set of all (ordered) haplotype pairs consistent with the multilocus genotype G i . Note that this likelihood is just the probability of observing the sample genotypes, as a function of the population haplotype frequencies, under the assumption of Hardy-Weinberg equilibrium (HWE).
Before running the iteration, for each genotype, find all possible haplotype pairs that are consistent with the genotype. Given k markers, there are 2 k-1 possible haplotype pairs per genotype. In our implementation, if the haplotype was already generated, we will create a link to connect the haplotype and the genotype. Otherwise, a new haplotype will be generated and linked to the genotype.
We outline the EM and OSLEM algorithms as follows: Step1: Obtain an initial distribution for genotype observed to corresponding haplotype pairs. For example, equal distribution is commonly used but random generation is also possible.
Step3: Recalculate distributions for genotypes by Hardy-Weinberg equilibrium condition

D preN
Step 4: Recalculate distribution by optimal step length: Step 5: Go to step 2 until step size becomes less then a given small value (precision).
Where D preN , D N and D N-1 are array variables and λ is a constant.
The EM algorithm jumps over Step 4, so it always takes λ = 1. In order to generate global maximization, the EM (or OSLEM) procedures are usually repeated 100 or more times for different initial distributions. In our implementation, we generate the initial distributions randomly from the first one as the equal distributions.
We use the following procedure to calculate λ in step 4. (where C i is the counting number of genotype G i , S i is the sum of the products of haplotype counts for all haplotype pairs (b k , b l ) that belongs to G i , A k is the haplotype count for b k and A l is the haplotype count for b l ).
After the above two nested loops, calculate the average of λ (k,l) . If the average is bigger than the minimum of the bounds B (k,l) , take the minimum upper bound as λ. We use 2 as an upper bound for λ. If the average is less than 1, we set the average of λ (k,l) = 1.
This iteration procedure can be viewed as searching for a fixed-point. By trying to solve the fixed-point equation approximately, we obtain an almost optimal step length formula for λ.

Results
By our tests using real data and tailored data, this new algorithm runs about twice as fast and obtains the same results as the EM algorithm. To test the performance, we generated the data in Table 1 for a single run (equal initial distribution) and Table 2 for multiple runs (initial distri-bution generated randomly). The whole run procedure includes three steps: the input/output step, data manipulation step, and the haplotype frequency estimation step. In the following tables, we only consider the haplotype frequency estimation step. The tailored data is edited from our epidemiological data. The data set is available on our website.
We have applied OSLEM to reconstruct haplotypes for epidemiological data. Lipoprotein lipase (LPL) is a glycoprotein involved in the transformation of dietary lipids into sources of energy for peripheral tissues (e.g., heart, muscle, adipose tissue) [10]. We performed an exhaustive analysis of genotypes and haplotypes spanning the LPL gene in 186 subjects whose blood lipid levels conferred a high risk of atherosclerosis (hereby referred to as "cases") and in 185 controls with non-atherogenic blood lipid profiles. Those subjects, ages 35 to 74, are representative of the general population of Geneva, Switzerland, in 1999 and 2000. Lipoprotein lipase sequence variants were surveyed by first re-sequencing its 10 exons and introns/ flanking regions in a selected subgroup of the case-control sample, followed by measurement of the most common SNPs in all cases and controls. Haplotypes were reconstructed from the individual SNPs separately for cases, controls, and the total sample. The relative frequencies of the estimated haplotypes in cases and controls are shown in table 3.

Discussion
By optimizing the step length of the EM (expectation maximization) algorithm, we have developed an accurate and faster algorithm for haplotype frequency estimation. This algorithm has been used successfully for lipoprotein lipase (LPL) genotyping analysis. The genetic analysis of lipoprotein lipase (LPL) gene-variants and their relation to population based variance in lipid profiles is been published separately [13].
The theoretical analysis of global optimization is, in general, a challenging mathematical target. Thus, a rigorous analysis of the rate of convergence for OSLEM may be quite difficult, but is very important. For this reason, there is a need for both analytic work and further computer simulation work.
It may be the case that there are multiple local maximization points for the mathematical formulation of haplotype frequency estimation. In this case, it would be possible to devise an algorithm to determine the extent of local maximization, which, in turn, would allow one to determine whether the globe optimal solution has been obtained.
We have set up a web server to provide haplotype frequency estimation service. The URL is http://genome3.cpmc.columbia.edu/~genome/HDL/.

Author's contributions
CG and AM conceived the Lipoprotein Lipase Genotyping project. PZ developed the new optimal step length EM   Haplotype (LPL exons)* Frequencies within subgroups (OSLEM) algorithm. PZ and HS coded the algorithm and developed the web server. All authors read and approved the final manuscript.