Optimal Step Length EM Algorithm (OSLEM) for the estimation of haplotype frequency and its application in lipoprotein lipase genotyping
- Peisen Zhang^{1}Email author,
- Huitao Sheng^{1},
- Alfredo Morabia^{3} and
- T Conrad Gilliam^{1, 2}
DOI: 10.1186/1471-2105-4-3
© Zhang et al; licensee BioMed Central Ltd. 2003
Received: 16 September 2002
Accepted: 15 January 2003
Published: 15 January 2003
Abstract
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.
Keywords
Haplotype frequency estimation Expectation Maximization Algorithm Genotype Lipoprotein LipaseBackground
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-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 maximum-likelihood 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,
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.
Step2: Gene-Counting [11, 12] calculating haplotype frequencies from the haplotype pair distribution.
Step3: Recalculate distributions for genotypes by Hardy-Weinberg equilibrium condition
D _{ preN }
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 ) })
else
λ_{ ( 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 λ.
Results
Single Run Comparisons of OSLEM and EM:
Number of Loops | CPU time | ||||||
---|---|---|---|---|---|---|---|
# SNPs | Precision | OSLEM | EM | OSLEM: EM | OSLEM | EM | OSLEM: EM |
12 | 1e-9 | 305 | 621 | 0.49 | 0.45 | 0.80 | 0.56 |
13 | 1e-7 | 305 | 579 | 0.53 | 0.95 | 1.69 | 0.56 |
14 | 1e-7 | 183 | 343 | 0.53 | 0.78 | 1.62 | 0.48 |
15 | 1e-7 | 182 | 342 | 0.53 | 0.87 | 1.78 | 0.49 |
16 | 1e-9 | 523 | 1149 | 0.45 | 4.92 | 11.28 | 0.44 |
Multiple-Run Comparisons of OSLEM and EM:
Table 2.1 Multiple-Run Maximum λ = 2 Precision: 1e-7 Run: 1000 times | |||
---|---|---|---|
SNPs # | OSLEM | EM | OSLEM : EM |
12 | 3m10.30s | 5m4.43s | 62.51% |
13 | 8m20.61s | 14m33.71s | 57.30% |
14 | 11m29.68s | 20m38.47s | 55.69% |
15 | 12m50.12s | 23m38.03 | 54.31% |
16 | 37m58.03s | 1h17m41.10s | 48.87% |
Table 2.2 Multiple-Run Maximum λ = 2 Precision: 1e-7 Run: 100 times | |||
SNPs # | OSLEM | EM | OSLEM : EM |
12 | 19.59s | 28.51s | 65.21% |
13 | 54.18s | 88.05s | 61.53% |
14 | 65.65s | 118.96s | 55.19% |
15 | 75.65s | 130.90s | 57.79% |
16 | 228.39s | 426.47s | 53.55% |
OSLEM-reconstructed haplotypes:
Haplotype (LPL exons)* | Frequencies within subgroups | |||
---|---|---|---|---|
3 4 5 6 8 9 10 | All | Cases | Controls | |
0: | 0 0 000 000 00 00 00 | 0.4582 | 0.4933 | 0.4459 |
1: | 0 0 000 0v 0 00 00 00 | 0.1191 | 0.1224 | 0.1000 |
2: | 0 0 vvv v 0v 0v vv v 0 | 0.0767 | 0.0474 | 0.1077 |
3: | 0 v 000 000 vv v 0 0v | 0.0389 | 0.0361 | 0.0363 |
4: | 0 0 000 000 vv v 0 0v | 0.0320 | 0.0324 | 0.0271 |
5: | 0 0 000 000 0v v 0 0v | 0.0363 | 0.0279 | 0.0443 |
6: | 0 0 vvv v 0v vv v 0 0v | 0.0239 | 0.0250 | 0.0223 |
7: | 0 0 000 000 0v vv v 0 | 0.0240 | 0.0200 | 0.0248 |
8: | v 0 000 v 00 00 00 00 | 0.0200 | 0.0186 | 0.0152 |
9: | 0 0 vvv v 0v 00 00 00 | 0.0211 | 0.0172 | 0.0225 |
Totals: | 0.8502 | 0.8403 | 0.8461 |
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 (OSLEM) algorithm. PZ and HS coded the algorithm and developed the web server. All authors read and approved the final manuscript.
Declarations
Acknowledgements
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.
Authors’ Affiliations
References
- 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/338446PubMed CentralView ArticlePubMedGoogle Scholar
- 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/319501PubMed CentralView ArticlePubMedGoogle Scholar
- Excoffier L, Slatkin M: Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population. Mol Biol Evol 1995, 12: 921–927.PubMedGoogle Scholar
- 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.PubMedGoogle Scholar
- 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.PubMed CentralPubMedGoogle Scholar
- Clark AG: Inference of haplotypes from PCR-amplified samples of diploid populations. Mol Biol Evol 1990, 7: 111–122.PubMedGoogle Scholar
- 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-4View ArticlePubMedGoogle Scholar
- 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/323622PubMed CentralView ArticlePubMedGoogle Scholar
- 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/303069PubMed CentralView ArticlePubMedGoogle Scholar
- 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-8View ArticlePubMedGoogle Scholar
- Ceppellini RM, Siniscalco M, Smith CAB: The estimation of gene frequencies in a random mating population. Ann Hum Genet 1955, 20: 97–115.View ArticlePubMedGoogle Scholar
- Smith CAB: Counting methods in genetical statistics. Ann Hum Genet 1957, 21: 254–276.View ArticlePubMedGoogle Scholar
- 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.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.