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

Gilliam, T. Conrad

Morabia, Alfredo

Sheng, Huitao

Zhang, Peisen

PubMed

Columbia University Academic Commons

BioMed CentralBMC BioinformaticsssOpen AcceMethodology articleOptimal Step Length EM Algorithm (OSLEM) for the estimation of haplotype frequency and its application in lipoprotein lipase genotypingPeisen Zhang*1, Huitao Sheng1, Alfredo Morabia3 and T Conrad Gilliam1,2Address: 1Columbia Genome Center, Columbia University, New York NY 10032, USA, 2Departments of Genetics & Development, and Psychiatry, Columbia University, New York NY 10032, USA and 3Division d'Epidémiologie Clinique, Hôpitaux Universitaires de Genève, Geneva, SwitzerlandEmail: Peisen Zhang* - pz6@columbia.edu; Huitao Sheng - hs734@columbia.edu; Alfredo Morabia - Alfredo.Morabia@hcuge.ch; T Conrad Gilliam - tcg1@columbia.edu* Corresponding author    Haplotype frequency estimationExpectation Maximization AlgorithmGenotypeLipoprotein LipaseAbstractBackground: Haplotype based linkage disequilibrium (LD) mapping has become a powerful andcost-effective method for performing genetic association studies, particularly in the search forgenetic 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 toestimate 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 steplength EM (OSLEM) that accelerates the calculation. By optimizing approximately the step length ofthe EM algorithm, OSLEM can run at about twice the speed of EM. This algorithm has been usedfor lipoprotein lipase (LPL) genotyping analysis.Conclusions: This new optimal step length EM (OSLEM) algorithm can accelerate the calculationfor haplotype frequency estimation for genotyping data without pedigree information. An OSLEMon-line server is available, as well as a free downloadable version.BackgroundEstimation of haplotype frequencies from routine geno-typing data plays an important role in LD analysis, andcan be achieved by varied methods, including  Monte-Car-lo (Gibbs sampling [1], PHASE method [2]), EM [3-5] andSubtraction [6]. At least two studies have tested and com-pared some of these programs. Xu and his collaborators[7] have empirically evaluated and compared the accuracyof the Subtraction method [6], the expectation-maximiza-tion (EM) method, and the PHASE method [2] for esti-mating haplotype frequency and for predicting haplotypephase. Summarizing from the studies of Xu and his col-laborators in [7]: "Where there was near complete linkagedisequilibrium (LD) between SNPs (the NAT2 gene), allthree methods provided effective and accurate estimatesfor haplotype frequencies and individual haplotype phas-es. For a genomic region in which marked LD was notmaintained (the chromosome X locus), the computation-al methods were adequate in estimating overall haplotypefrequencies. However, none of the methods were accuratein predicting individual haplotype phases. The EM andthe PHASE methods provided better estimates for overallPublished: 15 January 2003BMC Bioinformatics 2003, 4:3Received: 16 September 2002Accepted: 15 January 2003This article is available from: http://www.biomedcentral.com/1471-2105/4/3© 2003 Zhang et al; licensee 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.Page 1 of 5(page number not for citation purposes)BMC Bioinformatics 2003, 4 http://www.biomedcentral.com/1471-2105/4/3haplotype frequencies relative to the Subtraction methodfor both genomic regions." The PHASE algorithm is ex-tremely slow. A comparison of run-times was reported forfive SNPs from the NAT2 gene [7]: Subtraction method,0.01 second; EM method, 13.48 seconds; and PHASEmethod, over 128 minutes. Zhang and his collaborators[8] pointed out that, " The PHASE method did not yieldsignificantly different results from a simple maximum-likelihood procedure." From the two comparative studies[7,8] and from one simulation study [9], it was shownconvincingly that the expectation maximization (EM) al-gorithm is accurate for estimation of haplotype frequen-cies.The EM algorithm is much faster than Monte-Carlo basedalgorithms. Whereas it is fast with single runs for a relativesmall number of SNPs, it can be slow with multiple runsand for large number of SNPs. The EM algorithm is an op-timization algorithm. In order to obtain a global maximi-zation, EM should run many times from varied startingpoints. This can be very time consuming. Although in thestandard EM algorithm procedure, the step length is theoptimal length under the conditional expectation, thestep length is not optimal in general. In order to run fasterand more accurately, we have developed a new algorithmusing numerical analysis technology, a so-called optimalstep length EM (OSLEM) to accelerate the calculation. Byoptimizing approximately the step length of the EM (ex-pectation maximization) algorithm, OSLEM can run atabout twice the speed of EM.AlgorithmWe start with the same premises and notations as Stephenet al [2]. Given n diploid individuals from a population,let G = (G1,..., Gn) denote the (known) genotypes for theindividuals, let H = (H1,..., Hn) denote the (unknown)corresponding haplotype pairs, let F = (F1,...,FM) denotethe set of (unknown) population haplotype frequencies(the M possible haplotypes are arbitrarily labeled 1,...,M).Here Hi, a random variable depends on F. Let i be theset of all (ordered) haplotype pairs consistent with themultilocus genotype Gi, and suppose the distribution ofHi on i will follow the Hardy-Weinberg equilibrium.The EM and OSLEM algorithms attempt to find the haplo-type F that maximizes the likelihood.Here,where i is the set of all (ordered) haplotype pairs con-sistent with the multilocus genotype Gi. Note that thislikelihood is just the probability of observing the samplegenotypes, as a function of the population haplotype fre-quencies, under the assumption of Hardy-Weinberg equi-librium (HWE).Before running the iteration, for each genotype, find allpossible haplotype pairs that are consistent with the gen-otype. Given k markers, there are 2k-1 possible haplotypepairs per genotype. In our implementation, if the haplo-type was already generated, we will create a link to con-nect the haplotype and the genotype. Otherwise, a newhaplotype will be generated and linked to the genotype.We outline the EM and OSLEM algorithms as follows:Step1: Obtain an initial distribution for genotype ob-served to corresponding haplotype pairs. For example,equal distribution is commonly used but random genera-tion is also possible.Step2: Gene-Counting [11,12] calculating haplotype fre-quencies from the haplotype pair distribution.Step3: Recalculate distributions for genotypes by Hardy-Weinberg equilibrium conditionDpreNStep 4: Recalculate distribution by optimal step length:DN = DN-1 + λ (DpreN - DN-1) where λ ≥ 1Step 5: Go to step 2 until step size becomes less then a giv-en small value (precision).Where DpreN, DN and DN-1 are array variables and λ is aconstant.The EM algorithm jumps over Step 4, so it always takes λ= 1. In order to generate global maximization, the EM (orOSLEM) procedures are usually repeated 100 or moretimes for different initial distributions. In our implemen-tation, we generate the initial distributions randomlyfrom the first one as the equal distributions.We use the following procedure to calculate λ in step 4.HHL F G F G Fiin( ) = ( ) = ( )=∏Pr | Pr | .1Pr | ,,G F F Fi b bb b i( ) =( )∈∑ 1 21 2 HHPage 2 of 5(page number not for citation purposes)BMC Bioinformatics 2003, 4 http://www.biomedcentral.com/1471-2105/4/3Do loop for every ambiguous genotype Gi (genotypeswith more than one heterozygous locus),Do loop for every haplotype pair (bk, bl) that belongs toGi,If (DpreN (k,l) - DN-1(k,l)) < 0Calculate the upper bound for B(k,l) = -DN-1(k,l) /(DpreN(k,l) - DN-1(k,l))elseλ(k,l) = 1 + Ci/[Si/(Ak + Al) - Ci]if λ(k,l) < 0λ(k,l) = 1,            (where Ci is the counting number of genotype Gi,Si is the sum of the products of haplotype counts for allhaplotype pairs (bk, bl) that belongs to Gi, Ak is the haplo-type count for bk and Al is the haplotype count for bl).          After the above two nested loops, calculate the average ofλ(k,l). If the average is bigger than the minimum of thebounds B(k,l), take the minimum upper bound as λ. Weuse 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 afixed-point. By trying to solve the fixed-point equation ap-proximately, we obtain an almost optimal step length for-mula for λ.ResultsBy our tests using real data and tailored data, this new al-gorithm runs about twice as fast and obtains the same re-sults as the EM algorithm. To test the performance, wegenerated the data in Table 1 for a single run (equal initialdistribution) and Table 2 for multiple runs (initial distri-bution generated randomly). The whole run procedure in-cludes three steps: the input/output step, datamanipulation step, and the haplotype frequency estima-tion step. In the following tables, we only consider thehaplotype frequency estimation step. The tailored data isedited from our epidemiological data. The data set isavailable on our website.We have applied OSLEM to reconstruct haplotypes for ep-idemiological data. Lipoprotein lipase (LPL) is a glycopro-tein involved in the transformation of dietary lipids intosources of energy for peripheral tissues (e.g., heart, mus-cle, adipose tissue) [10]. We performed an exhaustiveanalysis of genotypes and haplotypes spanning the LPLgene in 186 subjects whose blood lipid levels conferred ahigh risk of atherosclerosis (hereby referred to as "cases")and in 185 controls with non-atherogenic blood lipidprofiles. Those subjects, ages 35 to 74, are representativeof the general population of Geneva, Switzerland, in 1999and 2000. Lipoprotein lipase sequence variants were sur-veyed by first re-sequencing its 10 exons and introns/flanking regions in a selected subgroup of the case-controlsample, followed by measurement of the most commonSNPs in all cases and controls. Haplotypes were recon-structed from the individual SNPs separately for cases,controls, and the total sample. The relative frequencies ofthe estimated haplotypes in cases and controls are shownin table 3.DiscussionBy optimizing the step length of the EM (expectation max-imization) algorithm, we have developed an accurate andfaster algorithm for haplotype frequency estimation. Thisalgorithm has been used successfully for lipoprotein li-pase (LPL) genotyping analysis. The genetic analysis of li-poprotein lipase (LPL) gene-variants and their relation topopulation based variance in lipid profiles is been pub-lished separately [13].The theoretical analysis of global optimization is, in gen-eral, a challenging mathematical target.  Thus, a rigorousTable 1: Single Run Comparisons of OSLEM and EM:Number of Loops CPU time# SNPs Precision OSLEM EM OSLEM: EM OSLEM EM OSLEM: EM12 1e-9 305 621 0.49 0.45 0.80 0.5613 1e-7 305 579 0.53 0.95 1.69 0.5614 1e-7 183 343 0.53 0.78 1.62 0.4815 1e-7 182 342 0.53 0.87 1.78 0.4916 1e-9 523 1149 0.45 4.92 11.28 0.44EM = Expectation maximization algorithm, OSLEM = Optimal Step Length EM algorithm The unit for CPU time is millisecond. The precision is the sum of the absolute differences of the haplotype frequencies between two loops.Page 3 of 5(page number not for citation purposes)BMC Bioinformatics 2003, 4 http://www.biomedcentral.com/1471-2105/4/3analysis of the rate of convergence for OSLEM may bequite difficult, but is very important. For this reason, thereis a need for both analytic work and further computer sim-ulation work.It may be the case that there are multiple local maximiza-tion points for the mathematical formulation of haplo-type frequency estimation. In this case, it would bepossible to devise an algorithm to determine the extent oflocal maximization, which, in turn, would allow one todetermine whether the globe optimal solution has beenobtained.We have set up a web server to provide haplotype frequen-cy estimation service. The URL is http://genome3.cp-mc.columbia.edu/~genome/HDL/.Author's contributionsCG and AM conceived the Lipoprotein Lipase Genotypingproject. PZ developed the new optimal step length EMTable 2: Multiple-Run Comparisons of OSLEM and EM:Table 2.1 Multiple-Run Maximum λ = 2 Precision: 1e-7 Run: 1000 timesSNPs # OSLEM EM OSLEM : EM12 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 timesSNPs # OSLEM EM OSLEM : EM12 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%Table 3: OSLEM-reconstructed haplotypes:Haplotype (LPL exons)* Frequencies within subgroups3 4 5  6  8  9  10 All Cases Controls0: 0 0 000 000 00 00 00 0.4582 0.4933 0.44591: 0 0 000 0v0 00 00 00 0.1191 0.1224 0.10002: 0 0 vvv v0v 0v vv v0 0.0767 0.0474 0.10773: 0 v 000 000 vv v0 0v 0.0389 0.0361 0.03634: 0 0 000 000 vv v0 0v 0.0320 0.0324 0.02715: 0 0 000 000 0v v0 0v 0.0363 0.0279 0.04436: 0 0 vvv v0v vv v0 0v 0.0239 0.0250 0.02237: 0 0 000 000 0v vv v0 0.0240 0.0200 0.02488: v 0 000 v00 00 00 00 0.0200 0.0186 0.01529: 0 0 vvv v0v 00 00 00 0.0211 0.0172 0.0225Totals: 0.8502 0.8403 0.8461On 14 Single Nucleotide Polymorphisms of the lipoprotein lipase gene in 163 atherogenic cases and 157 non-atherogenic controls. Geneva, Swit-zerland, 1999–2000. * SNPs are divided into subgroups according to their locations on exons, for example the first SNP is located on exon3, the second on exon 4, the third, fourth, and fifth are on exon 5, and so on. In the haplotype sequences, v means minor variation and 0 means major var-iation.Page 4 of 5(page number not for citation purposes)BMC Bioinformatics 2003, 4 http://www.biomedcentral.com/1471-2105/4/3Publish with BioMed Central   and  every scientist can read your work free of charge"BioMed Central will be the most significant development for disseminating the results of biomedical research in our lifetime."Sir Paul Nurse, Cancer Research UKYour research papers will be:available free of charge to the entire biomedical communitypeer reviewed and published immediately upon acceptancecited in PubMed and archived on PubMed Central yours — you keep the copyrightSubmit your manuscript here:http://www.biomedcentral.com/info/publishing_adv.aspBioMedcentral(OSLEM) algorithm. PZ and HS coded the algorithm anddeveloped the web server. All authors read and approvedthe final manuscript.AcknowledgementsWe 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.References1. Niu T, Qin ZS, Xu X and Liu JS Bayesian haplotype inference formultiple linked single-nucleotide polymorphisms. Am J HumGenet 2002, 70(1):157-692. Stephens M, Smith NJ and Donnelly P A new statistical methodfor haplotype reconstruction from population data. Am J HumGenet 2001, 68:978-9893. Excoffier L and Slatkin M Maximum-likelihood estimation ofmolecular haplotype frequencies in a diploid population. MolBiol Evol 1995, 12:921-9274. Hawley M and Kidd K HAPLO: a program using the EM algo-rithm to estimate the frequencies of multi-site haplotypes. JHered 1995, 86:409-4115. Long JC, Williams RC and Urbanek M An E-M algorithm and test-ing strategy for multiple locus haplotypes. Am J Hum Genet1995, 56:799-81036. Clark AG Inference of haplotypes from PCR-amplified sam-ples of diploid populations. Mol Biol Evol 1990, 7:111-1227. Xu CF, Karen Lewis, Kathryn Cantone L, Parveen Khan, ChristineDonnelly, Nicola White, Nikki Crocker, Pete Boyd R, Dmitri ZaykinV and Ian Purvis J Effectiveness of computational methods inhaplotype prediction Hum Genet 2002, 110(2):148-1568. Zhang S, Andrew Pakstis J, Kenneth Kidd K and Hongyu Zhao Com-parisons of Two Methods for Haplotype Reconstruction andHaplotype Frequency Estimation from Population Data Am JHum Genet 2001, 69:906-9129. Fallin D and Schork NJ Accuracy of haplotype frequency esti-mation for biallelic loci, via the expectation-maximization al-gorithm for unphased diploid genotype data. Am J Hum Genet2000, 67(4):947-5910. Murthy V, Julien P and Gagne C Molecular pathobiology of thehuman lipoprotein lipase gene. Pharmacol Ther 1996, 70:101-13511. Ceppellini RM, Siniscalco M and Smith CAB The estimation ofgene frequencies in a random mating population. Ann HumGenet 1955, 20:97-11512. Smith CAB Counting methods in genetical statistics. Ann HumGenet 1957, 21:254-27613. Morabia A, Cayanis E, Costanza MC, Ross BM, Bernstein MS, FlahertyMS, Alvin GB, Das K, Morris MA, Penchaszadeh GK, Zhang P and Gil-liam TC  Association between the lipoprotein lipase (LPL)gene and blood lipids: a common variant for a common trait.Genetic Epidemiol 2003, Page 5 of 5(page number not for citation purposes)

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

https://academiccommons.columbia.edu/doi/10.7916/D8RV1036/download

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

Abstract

Similar works

Full text

Available Versions

Columbia University Academic Commons