Optimal Step Length EM Algorithm (OSLEM) for the estimation of haplotype frequency and its application in lipoprotein lipase genotyping
BMC Bioinformatics volume 4, Article number: 3 (2003)
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.
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.
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.
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 , PHASE method ), EM [3-5] and Subtraction . At least two studies have tested and compared some of these programs. Xu and his collaborators  have empirically evaluated and compared the accuracy of the Subtraction method , the expectation-maximization (EM) method, and the PHASE method  for estimating haplotype frequency and for predicting haplotype phase. Summarizing from the studies of Xu and his collaborators in : "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 : Subtraction method, 0.01 second; EM method, 13.48 seconds; and PHASE method, over 128 minutes. Zhang and his collaborators  pointed out that, " The PHASE method did not yield significantly different results from a simple maximum-likelihood procedure." From the two comparative studies [7,8] and from one simulation study , 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.
We start with the same premises and notations as Stephen et al . Given n diploid individuals from a population, let G = (G1,..., G n ) denote the (known) genotypes for the individuals, let H = (H1,..., H n ) denote the (unknown) corresponding haplotype pairs, let F = (F1,...,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.
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 2k-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
Step 4: Recalculate distribution by optimal step length:
D N = D N-1 + λ (D preN - D N-1 ) where λ ≥ 1
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.
Do loop for every ambiguous genotype G i (genotypes with more than one heterozygous locus),
Do loop for every haplotype pair (b k , b l ) that belongs to G i ,
If (D preN ( k,l ) - D N-1( k,l ) ) < 0
Calculate the upper bound for B ( k,l ) = -D N-1( k,l ) / (D preN( k,l ) - D N-1( k,l ) )
λ ( k,l ) = 1 + C i /[ S i /( A k + A l ) - C i ]
if λ ( k,l ) < 0
λ ( k,l ) = 1,
(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 λ.
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 distribution 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) . 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.
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 .
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/.
CG and AM conceived the Lipoprotein Lipase Genotyping project. PZ developed the new optimal step length EM (OSLEM) algorithm. PZ and HS coded the algorithm and developed the web server. All authors read and approved the final manuscript.
Niu T, Qin ZS, Xu X, Liu JS: Bayesian haplotype inference for multiple linked single-nucleotide polymorphisms. Am J Hum Genet 2002, 70(1):157–69. 10.1086/338446
Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data. Am J Hum Genet 2001, 68: 978–989. 10.1086/319501
Excoffier L, Slatkin M: Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population. Mol Biol Evol 1995, 12: 921–927.
Hawley M, Kidd K: HAPLO: a program using the EM algorithm to estimate the frequencies of multi-site haplotypes. J Hered 1995, 86: 409–411.
Long JC, Williams RC, Urbanek M: An E-M algorithm and testing strategy for multiple locus haplotypes. Am J Hum Genet 1995, 56: 799–8103.
Clark AG: Inference of haplotypes from PCR-amplified samples of diploid populations. Mol Biol Evol 1990, 7: 111–122.
Xu CF, Karen Lewis, Kathryn CantoneL, Parveen Khan, Christine Donnelly, Nicola White, Nikki Crocker, Pete Boyd, Dmitri Zaykin, Ian Purvis J: Effectiveness of computational methods in haplotype prediction. Hum Genet 2002, 110(2):148–156. 10.1007/s00439-001-0656-4
Zhang S, Andrew Pakstis J, Kenneth Kidd K, Hongyu Zhao: Comparisons of Two Methods for Haplotype Reconstruction and Haplotype Frequency Estimation from Population Data. Am J Hum Genet 2001, 69: 906–912. 10.1086/323622
Fallin D, Schork NJ: Accuracy of haplotype frequency estimation for biallelic loci, via the expectation-maximization algorithm for unphased diploid genotype data. Am J Hum Genet 2000, 67(4):947–59. 10.1086/303069
Murthy V, Julien P, Gagne C: Molecular pathobiology of the human lipoprotein lipase gene. Pharmacol Ther 1996, 70: 101–135. 10.1016/0163-7258(96)00005-8
Ceppellini RM, Siniscalco M, Smith CAB: The estimation of gene frequencies in a random mating population. Ann Hum Genet 1955, 20: 97–115.
Smith CAB: Counting methods in genetical statistics. Ann Hum Genet 1957, 21: 254–276.
Morabia A, Cayanis E, Costanza MC, Ross BM, Bernstein MS, Flaherty MS, Alvin GB, Das K, Morris MA, Penchaszadeh GK, Zhang P, Gilliam TC: Association between the lipoprotein lipase (LPL) gene and blood lipids: a common variant for a common trait. Genetic Epidemiol 2003, in press.
We would like to thank Dr. Joseph Terwilliger, Dr. Hank Juo and Dr. Haghighi for helpful discussion, and Dr Michael C Costanza for performing the application. We would like to thank the reviewers and the editors for their comments and suggestions.
About this article
Cite this article
Zhang, P., Sheng, H., Morabia, A. et al. Optimal Step Length EM Algorithm (OSLEM) for the estimation of haplotype frequency and its application in lipoprotein lipase genotyping. BMC Bioinformatics 4, 3 (2003). https://doi.org/10.1186/1471-2105-4-3