Process, software arrangement and computer-accessible medium for obtaining information associated with a haplotype

[0001] The present application claims priority from U.S. patent application Ser. No. 60/557,768, filed Mar. 30, 2004, entitled “System, Method and Software Arrangement for bi-allele haplotype phasing,” the entire disclosure of which is incorporated herein by reference.

[0002] 1. Field of the Invention

[0003] The present invention relates to a system, method, and software arrangement for determining co-associations of allele types across consecutive genetic loci, permitting the haplotyping of genetic samples at one or more contiguous loci. The system, method, and software arrangement described herein can employ genotype data generated from a wide range of mapping methods to determine chromosomal haplotypes. As such, the system, method, and software arrangement of the present invention may be useful as an aid to the diagnosis and treatment of any disease which has a genetic component.

[0004] 2. Background Information

[0005] Diploid organisms are those whose somatic cells contain two copies of each chromosome. Each of these two copies of a particular chromosome may be distinguished by the presence, within the DNA which comprises the chromosome, of certain genetic variations, which may include restriction fragment length polymorphisms (RFLPs), single nucleotide polymorphisms (SNPs), sequence tag sites (STSs), microinsertions, microdeletions, or variable numbers of tandemly repeated elements (VNTRs). Based on the presence or absence of particular polymorphisms at specific loci, chromosomes may be assigned to one of two “haplotypes,” the name used to refer to the collection of identifiable genetic features present on one of the two haploid chromosomes that are contained within the diploid set. In certain situations, particular haplotypes may be associated with the presence or absence of a particular mutation or other functional variation in specific genetic loci. Because these genetic variations or genotypes may be associated with certain disease states or even with predisposition to disease, determination of relationships between haplotypes and genotypes are of intense interest in genetic research.

[0006] One of the most difficult problems in determining haplotypes in diploid organisms is establishing the proper assignment of multiple polymorphic markers to the same chromosome. Thus, distribution of two different variants (or alleles) at two different genetic loci, for example A/a and B/b, could generate haplotypes AB, Ab, aB or ab depending upon the distribution of the two alleles between the two chromosomes. The problems of inferring diploid haplotypes through the use of population data have been extensively investigated and widely acknowledged. See Clark, 1990, Mol. Biol. Evol. 7:111-122; Excoffier and Slatkin, 1995, Mol. Biol. Evol. 12:921-927; . Ma et al. 2000, Neural Computation 12:2881-2907; Gusfield, 2001, J. Computational Biology 8:305-323; Stephens et al., 2001, Am. J. Hum. Genet. 68:978-989; and Niu et al., 2002, Am. J. Hum. Genet. 70:156-169.

[0007] The present invention relates to a system, method, and software arrangement for determining co-associations of allele types across consecutive genetic loci, permitting the haplotyping of genetic samples at one or more contiguous loci. The system, method, and software arrangement described herein can employ genotype data generated from a wide range of mapping methods to determine chromosomal haplotypes. According to one exemplary embodiment of the present invention the system, method, and software arrangement for haplotype reconstruction of the instant invention detects polymorphic marker types at one or more contiguous genetic loci and then establishes the maximum likelihood that an arbitrary assignment of polymorphic markers to a particular haplotype adequately account for the particular genotypes observed at any given locus. The system, method, and software arrangement described herein may be useful as an aid to the diagnosis and treatment of any disease that has a genetic component associate therewith.

[0008] In contrast to these previous approaches, the present invention addresses the use of multiple independent mapping techniques (for example, those performed on a collection of large DNA fragments) as a source of the base data used to infer haplotypes. Such single molecule methods and technologies, which include optical mapping and “polony” (See Mitra et al., 1999, “In Situ Localized Amplification and Contact Replication of Many Individual DNA Molecules,” Nucleic Acids Research 27:e34-e34), may permit the application of high-throughput methodologies to the haplotyping of a diploid individual in a population. Moreover, unlike the previous haplotyping approaches described above, the present invention determines the haplotype of an individual without reliance upon the haplotype statistics of the population to which the individual belongs. A further advantage of the instant method is its potential to determine an individual's haplotype with accuracy and resolution far beyond methods based on population studies, both in the worst-case and in the average-case scenarios.

[0009] Thus, in accordance with an exemplary embodiment of the present invention, a system, process and software arrangement are provided for obtaining information associated with a haplotype of one or more genetic samples from genotype data. In particular, the genotype data is received, and polymorphic genetic markers are identified from the genotype data. Further, at least one association of the polymorphic genetic markers across genetic loci (which can be consecutive) may be determined to obtain the information.

[0010] In one further embodiment of the present invention, the genotype data can be obtained from the corresponding one or more genetic samples, the genotype data being contained in at least one dataset. The association can be determined at one or more contiguous loci. The determined association may constitute the haplotype of the one or more genetic samples. The polymorphic genetic markers may be restriction fragment length polymorphisms (“RFLPs”), single nucleotide polymorphisms (“SNPs”), sequence tag sites (“STSs”), insertions, microinsertions, deletions, microdeletions, variable numbers of tandemly repeated elements (“VNTRs”), microsatellites, expanded repeats of variable repeat number, or any combination thereof.

[0011] In still another exemplary embodiment of the present invention, the genotype data can be obtained from single DNA molecules with locations of polymorphic genetic markers. The polymorphic genetic markers are determined using at least one of a Maximum Likelihood Estimator (“MLE”) algorithm and an Expectation Maximization (“EM”) algorithm. The polymorphic genetic markers may be defined as events, and

[0012] $X\left(D\right)=\{\begin{array}{cc}1& \mathrm{if}\hat{\Phi }\left(D\right)\text{:}{\hat{\mu }}_1-{\hat{\mu }}_2-\delta >0\\ 0& o.w.\end{array},$
where {circumflex over (Φ)} (D) denotes the limit of the EM-algorithm with data set D at the loci j.

[0013] For example, the associations of the polymorphic genetic markers across consecutive genetic loci at one or more contiguous loci are determined through the formulation of a maximum likelihood problem. The likelihood function of the maximum likelihood problem can be given by:

[0014] $L\left(\theta \right)=P\left(D|\theta \right)=\frac{\Gamma \left(N\right)}{\prod _{\rho \in M}\Gamma \left(N{\alpha }_p\right)}\prod _{\rho \in M}{\theta }_{\rho }^{{\alpha }_{p}N},$
and the maximum likelihood estimation can include a finding
ρ ε A so that L(ρ)≧L ω) ∀ω ε A.

[0015] For a set of genetic loci {j₁, j₂, . . . ,j_v}, a likelihood function L_{M[j1,j2, . . . jv]} can be defined as that most likely to produce posterior α_{j1, j2. . . jv} over the space M[j₁, j₂, . . . j_v]. The associations of the polymorphic genetic markers across the genetic loci at one or more contiguous loci are determined through maximizing
Π_ρε{0,1}M-1θ_ρ^αρN over θ ε
or minimizing

[0016] $\sum _{j\in \left[1\dots M\right]}\frac{{\left({\alpha }_j-{\theta }_j\right)}^2}{{\alpha }_j}\mathrm{over}\theta \in A.$

[0017] The associations of the polymorphic genetic markers at one or more loci, thereby producing a contig, can be determined through an application of a VERIFY-PHASE function. The VERIFY-PHASE FUNCTION may determine one or more selected phasing criteria. The phasing criteria may be a statistically-significant rejection of Hardy-Weinberg Equilibrium. The association of the polymorphic genetic markers over the contig can be computed through application of a MLE-COLLAPSE function. Further, the association of polymorphic genetic markers at one or more contiguous loci, thereby producing a contig, can be determined through an application of a COMPUTE-PHASE function or a JOIN function.

[0018] According to yet another exemplary embodiment of the present invention, a linear number of the polymorphic genetic markers can be examined. The information associated with the haplotype data may be determined without a need to obtain at least one of first data associated with a pedigree and second data associated with a sub-population of individuals. In addition, further information associated with a confidence level of accuracy of the information can be obtained. A resolution and the further information may be proportional to an amount of effort associated with the determination of the information.

[0025] The problem of reconstructing two haplotypes from genotype data generated by general mapping techniques can be considered according to the exemplary embodiment of the present invention, with a focus on single molecule methods. The genotype data is a set of observations D=<d_i>_{iε[1 . . . N]}. Each observation is derived from one of the two distinct but unknown haplotypes. Each observation di=<d_ij>_{jε[1 . . . M]} is a set of observations over the loci index j with d_ij ε R^r.

[0026] Mapping processes are subject to noise, for which a Gaussian model d_ij˜N(μ,σ) with parameter μ depending on the underlying haplotype of d_i can be assumed. Mapping processes are designed to discriminate the polymorphic allele types in the data space for each locus. Hence, the set of observation points <d_ij>_{iε[1 . . . N]} can be derived from a mixed distribution, which may display bimodal characteristics in the presence of a polymorphic feature. By estimating the parameters of the distribution, a posteriori distribution that a particular point in R^r may be derived from an allele type can be assigned.

[0027] Since the mapping errors for d_ij and d_ij, can be assumed to be independent, determining and computing the posteriori distribution for haplotypes with product allele types is likely straightforward, and is one of the advantages of utilizing single molecule methods in association studies.

[0028] One of the problems of phasing is to determine which haplotypes are most likely to account for the observed genotype data. It is preferable to establish the phase by inferring the most likely parameter correlations across the loci index accounting for the posteriori distribution.

Mapping Techniques

[0029] The system, method and software arrangement according to an exemplary embodiment of the present invention is applicable to datasets generated by a wide spectrum of mapping techniques. In this manner, a large number of polymorphic markers of different types (e.g. SNPs, RFLPs, micro-insertions and deletions, satellite copy numbers) can be used in an association procedure in accordance with the present invention. An exemplary embodiment of system, method and software arrangement presented herein can utilize mapping techniques capable of a) discriminating alleles at polymorphic loci, and b) providing haplotype data at multiple loci. Other known polymorphic genetic markers may be used by the system, method, and software arrangement according to the present invention.

[0030] A mapping technique designed for association studies should preferably be discriminating. In particular, for each polymorphic loci, data points in the data space R^r which are derived from separate allele types should preferably form distinct clusters in the data. One of the exemplary techniques which allows an observation of a single haplotype over multiple loci may preferably be used for an efficient phasing procedure that can be used. For example, single molecule methods may be of exemplary interest. The models and analysis presented herein relates to and is effected by such methods' applicability to association studies.

[0031] As an example, the length between two restriction fragments may be considered. The observable x may be modeled as a random variable depending on the actual distance μ.

[0032] $P\left(x|\mu \right)=\frac{1}{\sqrt{2\pi \mu }}\exp \left(\frac{-{\left(x-\mu \right)}^2}{2\mu }\right)$

[0033] Then, it is possible to isolate a specific pair of restriction sites on one of the haplotypes H₁, and let the distance between them be provided by μ₁. The distance between the homologous pair on the second haplotype H₂ can be provided by μ₂. An observation x from the genotype data can then be either derived from H₁ or H₂, denoted x˜H₁ and x˜H₂, respectively.

[0034] $\begin{array}{c}P\left(x\right)=P\left(x|x\sim H_1\right)P\left(x\sim H_1\right)+P\left(x|x\sim H_2\right)P\left(x\sim H_2\right)\\ =\frac{1}{\sqrt{2\pi {\mu }_1}}\exp \left(\frac{-{\left(x-{\mu }_1\right)}^2}{2{\mu }_1}\right)P\left(x\sim H_1\right)+\\ \frac{1}{\sqrt{2\pi {\mu }_2}}\exp \left(\frac{-{\left(x-{\mu }_2\right)}^2}{2{\mu }_2}\right)P\left(x\sim H_2\right).\end{array}$

[0035] Using the RFLP sizing mapping technique, observable d_ij, can have independent error sources depending on loci-specific parameters. The set {d_ij, i ε [1 . . . N]} can provide points in R which may be discriminated using, e.g., a Gaussian Mixture model. Due to the uncertainty of mapping and underlying haplotypes, posteriori distribution α(x)=[P(x˜H₁),P(x˜H₂)] can be selected to model the data rather than determined allele types.

[0036] The additional description of the exemplary embodiment of the present invention is organized into the following four sections: Section 1 describes the EM-Algorithm procedure used with the exemplary embodiment of the present invention, Section 2 discusses the phasing problem, addressed by the present invention, Section 3 describes a procedure implementation and examples thereof, and Section 4 describes results and applications of the present invention.

[0037] Section 1. EM-Algorithm for Detection of Bi-Allelic Polymorphisms of the Present Invention

[0038] The use of the EM-Algorithm for inferring parameters of a Gaussian mixture model is a well known technique (see Dempster et al., 1977, J. Roy. Stat. Soc. 39:1-38; and Roweis et al. 1999, Neural Computation 11:305-345), and, as described herein, also useful in the detection of biallelic polymorphisms. In the presence of polymorphisms at loci j, informative mapping data can provide a bimodal distribution in the data space R^r. Detailed exemplary computations for E-Step and M-Step of the EM-Algorithm are described herein in the Examples section. For each locus j, the EM-algorithm can be executed in accordance with the present invention until convergence occurs, the result being: <α_k(x),{circumflex over (Φ)}=<{circumflex over (μ)}₁,μ₂, . . . ,μ_K,σ>>. In this example, α is a posteriori probability that data point x is derived from allele type k ε [1,2, . . . , K].

[0039] Criteria for Polymorphisms. Let {circumflex over (Φ)} (D) denote the limit of the EM-algorithm with data set D at the loci j. The following issue may be raised: when will a locus exhibit two specific allele variations? In the setting when K=2, as in the remainder of this description (hence {circumflex over (Φ)}=<μ₁,μ₂,σ>), polymorphic loci are defined as events:

[0040] $X\left(D\right)=\{\begin{array}{cc}1& \mathrm{if}\hat{\Phi }\left(D\right)\text{:}{\hat{\mu }}_1-{\hat{\mu }}_2-\delta >0\\ 0& o.w.\end{array}$

[0041] Effectiveness of the EM-Algorithm. Mapping techniques may contain errors that are Gaussian across a diverse set of technologies. Genetic markers may be associated or linked to allele types in the population. The mixture model/technique treated with EM Algorithm can operate effectively, and possibly distinguishing fits beyond visual accuracy. The constraint for a single value of α can force the EM-Algorithm to result in one of two steady states, e.g., μ₁≠μ₂ or μ₁=μ₂ (a single Gaussian). Although the EM-Algorithm estimates are slightly biased, the estimators are consistent, and the bias is known to diminish with larger data sets.

[0042] The individual experimental data {d_ij:i ε[1 . .. N]} can be mapped to posteriori probability that measures over the allele classes. Thus, a probability function a (y) reflecting a confidence (in the presence of mapping error) that point y corresponds to one of our allele types can be produced. For polymorphism assignments, false positives are unlikely to disturb the phasing, while false negatives may affect the size of phased contigs.

[0043] Section 2. Phasing Genotype Data

[0044] “Phasing” is the problem of determining the association of alleles, due to a linkage on the same haplotype. Letting Λ_j be the allele space at loci j, a haplotype may be considered an element of the set: Π_{j ε[1,2, . . . ,M]}Λ_j.

[0045] In phasing polymorphic alleles for an individual's genotype data (a mix of two haplotypes), it can be assumed that about half of the data can be derived from each of the underlying haplotypes H₁ and H₂. In this context, haplotypes have a complementary structure in that the individual's genotype should be heterozygous at each polymorphic locus.

[0046] A haplotype space can be defined in accordance with the present invention, and methods for estimating the probability that an observation di is derived from a particular haplotype over a set of loci are described. The maximum likelihood problem for haplotype inference can then be formulated, and this formulation may be the proposed solution to the phasing problem.

[0047] Haplotype Space and Joint Distributions. The full space of haplotypes is the product over all allele spaces {1, 2, . . . ,M}. In general, haplotype space is in one-to-one correspondence with M={−1,l}^M. The discrete-measure space (M,2^M) can be used to denote the haplotypes, while M[j₁, j₂, . . . , j_v] denotes the haplotypes over the range of loci j₁, j₂, . . . , j_v. The result of phasing genotype data may be a probability measure on the space (M, 2^M). Noiseless data may result in a measure assigning ½ to each of the complementary haplotypes, and 0 to all others. This uniform measure over complements corresponds to perfect knowledge of what the haplotypes are. The procedure that can be used by an exemplary embodiment of the system, method and software arrangement of the present invention is consistent in that the correct result is achieved for suitably large data sets.

[0048] For example, let Λ_j be the allele set for the polymorphic loci j. Two bi-allelic loci j and j′ can be used. For clarity, we will assume that Λ_j={A,a} while Λ_j′={B,b}. A data observation d_i may be derived from one of the four classes: AB, Ab, aB, ab. Because the mapping noise at loci j and j′ are independent, the probability (based on the loci a posteriori) that the observation can be derived from the following four classes:
P(d_i˜AB)=α_jA(d_ij)α_j′B(d_ij′)
P(d_i˜Ab)=α_jA(d_ij)α_j′b(d_ij′)=αjA(d_ij)(1−α_j′B(d_ij′))
P(d_i˜aB)=α_ja(d_ij)α_j′B(d_ij′)=(1−α_jA(d_ij))α_j′B(d_ij′)
P(d_i˜ab)=α_ja(d_ij)α_j′b(d_ij′)=(1−α_jA(d_ij))(1−α_j′B(d_ij′))

[0049] α_jj′⁽ⁱ⁾ can be defined as the estimated probability distribution for observation i on haplotypes over the loci j, j′:
α_jj′⁽ⁱ⁾=[α_jj′AB(d_i),α_jj′Ab(d_i), α_jj′aB(d_i), α_jj′ab(d_i)]

[0050] α_jj′ is defined as the estimated probability distribution over the data set on haplotypes over the loci j, j′:

[0051] ${\alpha }_{j{j}^{\prime }}\left(D\right)=\frac{1}{N}\sum _{i=1}^N{\alpha }_{j{j}^{\prime }}^{\left(i\right)}$

[0052] For ρ ε M[j₁, j₂ . . . , j_M] and α_jwρw (d_i)=Prob(d_i˜ρ_w) with ρ_w ε Λ_jw, the estimates can be extended to any set of indices producing:

[0053] $\begin{array}{c}{\alpha }_{j_1{j}_2\dots j_w}^{\left(i\right)}={\left[\prod _{w\in \left[1\dots v\right]}^{}{\alpha }_{j_w{\rho }_w\left(d_i\right)}\right]}_{\rho \in M\left[j_1,j_2,\dots ,j_v\right]}\\ {\alpha }_{j_1{j}_2\dots j_v}=\frac{1}{N}\sum _i{\alpha }_{j_1{j}_2\dots j_v}^{\left(i\right)}\end{array}$

[0054] Complementarity. In phasing the diploid genotype data into two haplotypes ρ₁, ρ₂ ε M, there may be a special property present, e.g., haplotype ρ₂ can be complementary to haplotype ρ₁, denoted ρ₂=ρ₁. The complementary pair of haplotypes may be represented by a change of variables ω ε {−1,1}^M-1, and the transformation to the haplotypes may be given by the map:

[0055] $\begin{array}{c}{\rho }_1\left(b\right)=\{\begin{array}{cc}-1& \mathrm{if}b=1\\ -1& \prod _{j=1:\left(b-1\right)}^{}w\left(j\right)\mathrm{for}b\in \left[2\dots M\right]\end{array}\\ {\rho }_2\left(b\right)=\{\begin{array}{cc}1& \mathrm{if}b=1\\ 1& \prod _{j=1:\left(b-1\right)}^{}w\left(j\right)\mathrm{for}b\in \left[2\dots M\right]\end{array}\end{array}$
In evaluating the data, there may be a possible 2^M-1 complementary pairs of allele types to search.

[0056] The confidence of a set of complementary haplotypes can be modeled as a probability distribution on the discrete measure space (M, 2^M), which is the convex hull of the following set of extremal points which correspond to certain knowledge of complementary haplotypes.

[0057] $A=\{{\theta }_{\rho }:{\theta }_{\rho }\left(\delta \right)=\left\{\begin{array}{cc}\frac{1}{2}& \mathrm{if}\delta =\rho \\ \frac{1}{2}& \mathrm{if}\delta =\stackrel{\_}{\rho }\\ 0& o.w.\end{array}\mathrm{for}\delta \in M\right\}$

[0058] These values can represent the uniform distribution over complementary haplotypes and geometrically are vertices of a high dimensional hypercube. For example, let A[j₁, j₂, . . . , j_v] be the corresponding distribution over the haplotype space M[j₁, j₂, . . . , j_v].

[0059] Maximum Likelihood Problem. For every loci j, it can be assumed that the data {d_ij:i ε [1 . . . N]} contains an equal distribution of data from the underlying haplotypes H₁, H₂ that can be inferred. Using the estimated values α for the joint distribution over loci product spaces, the haplotypes most likely producing a can be computed. The corresponding maximum likelihood problem may be formulated as follows:

[0060] Let the likelihood function be given by:

[0061] $\begin{array}{c}L\left(\Theta \right)=P\left(D|\Theta \right)=\frac{\Gamma \left(N\right)}{\prod _{\rho \in M}\Gamma \left(N{\alpha }_p\right)}\prod _{\rho \in M}{\Theta }_{\rho }^{{\alpha }_{p}N}\\ \mathrm{MLE}1\mathrm{Find}\rho \in A\mathrm{so}\mathrm{that}L\left(\rho \right)\ge L\left(\omega \right)\forall \omega \in A.\end{array}$

[0062] Similarly, for any specified set of loci {j₁, j₂, . . . , j_v}, a likelihood function L_{M[j1,j2, . . . jv}] maybe defined as the most likely to produce posterior α_{j1, j2 . . . jv} over the space M[j₁, j₂, . . . j_v].

[0063] $\begin{array}{c}\mathrm{Lemma}1\mathrm{If}d\left(\alpha ,A\right)<\in \mathrm{for}\mathrm{some}\in \mathrm{small}\mathrm{enough},\\ \mathrm{and}d\left(\alpha ,A\right)={\min }_{\theta \in A}{\alpha -\theta }_2.\\ \mathrm{Maximizing}\prod _{\rho \in {\left\{0,1\right\}}^{M-1}}^{}{\Theta }_{\rho }^{{\alpha }_{p}N}\mathrm{over}\Theta \in A\\ \mathrm{is}\mathrm{equivatent}\mathrm{to}\mathrm{minimizing}\\ \sum _{j\in 1\dots M}\frac{{\left({\alpha }_j-{\Theta }_j\right)}^2}{{\alpha }_j}\mathrm{over}\Theta \in A.\end{array}$

[0064] The description provided below in the Examples section is derived from a Taylor-series expansion of the likelihood function. It demonstrates that the MLE result in set A is the vertex of a 2^M-1 hyper-cube closest to the estimated joint probability function α, measured by a modified L₂ norm.

[0065] With this result, the following function to be used in the algorithms presented later can be assumed:

[0066] $\begin{array}{c}\mathrm{Algorithm}1\\ \mathrm{MLE}\text{-}\mathrm{COLLAPSE}\left(j_1,j_2\dots ,j_v\right)\\ \mathrm{Compute}\rho \in A\left[j_1,j_2\dots ,j_v\right]\\ \mathrm{minimizing}\sum _{j\in j_1,j_2\dots ,j_v}\frac{{\left({\alpha }_j-{\Theta }_j\right)}^2}{{\alpha }_j}\mathrm{over}\Theta \in A\\ \mathrm{return}\rho \end{array}$

[0067] Section 3. Exemplary Procedures

[0068] Exemplary procedures generally focus on growing disjoint-phased contiguous sets of loci called contigs. For example, all loci can be assigned an arbitrary phase and begin as a singleton phased contig. A JOIN operation checks if these phased contigs may be phased relative to one another using a function called VERIFY-PHASE. VERIFY-PHASE can be designed to check a phasing criteria, for example refuting a hypothesis of Hardy-Weinberg Equilibria is discussed below in the Examples section. Other examples of suitable phasing criteria are known to those of ordinary skill in the art. Such criteria may, for instance, be based on the statistical distribution of haplotypes in the ambient population, on the perfect phylogeny hypothesis, or on the relation to genotypes of parents, siblings and other closely-related family members.

[0069] If a pair of phased contigs can be joined by passing the test, implied by VERIFY-PHASE function, then the disjoint sets are combined into a single phased contig and the joint distribution over the set is computed with the MLE-COLLAPSE function. After completion of a successful join operation, the resulting distribution function may be regarded as the most likely one among all haplotypes that can generate the observed data over the specified loci. Because the growth of contigs is monotonic and depends on local information available at the time of the operation, an ADJUST operation is also considered that fractures and rejoins contigs using a larger locality of data than what was available during the JOIN.

[0070] The operations are described in detail, the results are analyzed, and methods for avoidance of incorrect operations are indicated.

[0071] Collapse. The collapse of the phase operation can be described as the MLE-COLLAPSE function, the example of which is shown in FIG. 1 which provides a collapse of the phase of two polymorphic loci 30 into a single haplotype 40. It may be used to update a joint probability distribution over a set of contigs, and has the effect of keeping the contig structures bound to haplotype states which simplifies the computing of a phase. FIG. 2 shows an exemplary illustration of the phase of three polymorphic loci 70 provided into two complementary haplotypes 80 in accordance with an exemplary embodiment of the present invention

[0072] Join: Let K be a parameter denoting neighborhood size. Let C₁={j₁, j₂, . . . , j_v} and C₂+{j′₁, j′₂, . . . j′_w}; then the join operation is as follows:

[0073] Given joint-probability functions α_j1, α_j2, . . . α_jv, α_j′1, α_j′2, . . . α_j′w, compute the joint probability function α_{C1, C2} with formula
α_CB,C2=α_{(jv)(j′1j′2 . . . j′ω)}=ω_jv,j′1α_{(jω)(j′1j′2 . . . j′ω)}+ω_jvj′α_{(jv)(j′1j′2 . . . j′w)}+. . . +ω_jv,j′Kα_{(jv)(j′1j′2 . . . j′K . . . j′ω)}
with

[0074] ${\alpha }_{\left(j_v\right)(\left(j_1^{\prime }{j}_2^{\prime }\dots j_{\kappa }^{\prime }\dots j_w^{\prime }\right)}=\sum _{i=1:N}{\alpha }_{j_v{j}_x^{\prime }}^{\left(i\right)}=\sum _{i=1:N}{\alpha }_{j_u}^{\left(i\right)}\left(d_{{\mathrm{ij}}_u}\right){\alpha }_{j_x^r}^{\left(i\right)}\left(d_{{\mathrm{ij}}_x^{\prime }}\right)$ ${\omega }_{j_v,j_x^{\prime }}=\kappa \frac{1}{d\left(j_v,j_x^{\prime }\right)}$ $\mathrm{Here},\kappa =\frac{1}{\frac{1}{d\left(j_v,j_1^{\prime }\right)}+\frac{1}{d\left(j_v,j_2^{\prime }\right)}+\dots +\frac{1}{d\left(j_v,j_{\kappa }^{\prime }\right)}}\mathrm{and}d\left(j_v,j_x^{\prime }\right)\mathrm{is}\mathrm{proportional}\mathrm{to}\mathrm{genomic}\mathrm{distance}\mathrm{between}\mathrm{loci}{j}_v\mathrm{and}{j}_x^{\prime }.\mathrm{Algorithm}2\mathrm{COMPUTE}\text{-}\mathrm{PHASE}\left(C_1=\left\{j_1,j_2,\dots ,j_v\right\},C_2=\left\{j_1^{\prime },j_2^{\prime },\dots ,j_w^{\prime }\right\},\kappa \right)\mathrm{assume}{j}_v\mathrm{in}{C}_1\mathrm{is}\mathrm{such}\mathrm{that}d\left(j_v,C_2\right)\le d\left(j_1{C}_2\right)\forall j\in C_1;\mathrm{Compute}{\alpha }_{\left(j_v\right)\left(j_1^{\prime },j_2^{\prime },\dots ,j_w^{\prime }\right)}\mathrm{using}\mathrm{parameter}\kappa .\mathrm{return}{\alpha }_{\left(j_v\right)\left(j_1^{\prime },j_2^{\prime },\dots ,j_w^{\prime }\right)}\mathrm{Algorithm}3\mathrm{JOIN}\left(C_1=\left\{j_1,j_2,\dots ,j_v\right\},C_2=\left\{j_1^{\prime },j_2^{\prime },\dots ,j_w^{\prime }\right\},\kappa \right)\mathrm{COMPUTE}\text{-}\mathrm{PHASE}\left(C_1=\left\{j_1,j_2,\dots ,j_v\right\},C_2=\left\{j_1^{\prime },j_2^{\prime },\dots ,j_w^{\prime }\right\},\kappa \right)\mathrm{if}\left(\mathrm{VERIFY}\text{-}\mathrm{PHASE}\left({\alpha }_{j_1{j}_2\dots j_v}\right)\right)\mathrm{then}$ ${\alpha }_{j_1{j}_2\dots j_v}\leftarrow \mathrm{MLE}-\mathrm{COLLAPSE}\left(j_1{j}_2\dots j_v\right);$

[0075] The method, system and software arrangement of an exemplary embodiment according to the present invention can estimate the haplotypes by solving an ordered set of local MLE problems.

[0076] Implementation. Input. The input is a set of data points {d_ij ε R^r: i ε [1 . . . N], j ε [1 . . . M]}. The following assumptions are made about the input:

• • [0077] For each j the points d_1j,d_2j, . . . d_Nj are derived from the Gaussian mixture model corresponding to mapping data at polymorphic loci j.
  • [0078] For each i points d_i1, d_i2, . . . , d_iM are independent random variables with parameters associated to underlying haplotypes.

[0079] With the knowledge of the mapping order of polymorphic loci, the positions of the genome can be assumed to be χ₁, χ₂, . . . , χ_M.

[0080] Implementation. Pre-Process. The EM-algorithm procedure can be executed for each locus: {d_ij: i ε [1 . . . N] observable}→{{circumflex over (Φ)}_j: α_j}∀j ε [1, . . . , M].

[0081] The result is a set of estimates for bi-allelic loci, {{circumflex over (Φ)}₁, {circumflex over (Φ)}₂, . . . , {circumflex over (Φ)}_M}, as well as a set of functions estimating the probability that any data point derives from the distinct alleles {α₁, α₂, . . . , α_M}.

[0082] Next, a join schedule can be constructed. Letting β_j=χ_j+1−χ_j, the results are sorted into an index array giving an increasing sequence: j₁, j₂, . . . j_v, . . . j_m-1.

[0083] Implementation. Main Algorithm and Data Structure. Contigs can be maintained in a modified union-find data structure designed to encode a collection of disjoint, unordered sets of loci, which may be merged at any time. Union-find supports two operations, UNION and FIND (see Taijan, 1983, Data Structures and Network Algorithms, CBMS 44, SIAM, Philadelphia). For example, UNION can merge two sets into one larger set, and FIND can identify the set containing a particular element. Loci j may be represented by the estimated distribution α_J, and can reference its left and right neighbor. At any instant, a phased contig may be represented by:

• • [0084] An MLE distribution or haplotype assignment for the range of loci in the contig (if one can be evaluated).
  • [0085] Boundary loci: Each contig has a reference to left- and right-most locus.

[0086] In the vth step of the procedure, consider the set of loci determined by β_v, {j_v, j_v+1}: If FIND (j_v) and FIND (j_v+1) are in distinct contigs C_p and C_q, then (a) attempt to UNION C_p and C_q, by use of the JOIN operation, and (b) update the MLE distribution and boundary loci at the top level if the JOIN is successful.

[0087] Implementation. Output. Output can be a disjointed collection of sets, each of which is a phased contig. It represents the most likely haplotypes over that particular region.

[0088] Implementation. Time Complexity. The preprocess may involve using the EM-algorithm/procedure once for each locus. The convergence rate of the EM-algorithm procedure has been analyzed (see Ma et al., 2000, Neural Computation 12:2881-2907) and depends on the amount of overlap in the mixture of distributions. For moderate-sized data sets, no difficulties with convergence of the EM-algorithm procedure have been observed.

[0089] The time complexity of the main exemplary procedure can be estimated by implementing the K-neighbor version. For each β_jv there may be two find operations. The number of union operations preferably does not exceed the cardinality of the set {β_j: j ε [j₁, j₂, . . . j_M-1]}, as contigs grow monotonically. The time-cost of a single “find” operation is preferably at most γ(M), where γ is the inverse of Ackermann's function. Therefore, the time cost of all union-find operations is preferably at most O(Mγ(M)). The join operation, on the other hand, uses the execution of the K-neighbor optimization routine, at a cost of O(K). Thus, the main exemplary procedure has a worst-case time complexity of
O(M(γ(M)+K))=O(Mγ(M))
and may be regarded as approximately linear in the number of markers, M for all practical purposes, since K is likely a small constant.

[0090] The method, system and software arrangement of the exemplary embodiment of the present invention can be demonstrated on two simulated data sets composed of ordered restriction fragment lengths subject to sizing error. FIG. 3, which shows exemplary haplotypes, obtained using such exemplary embodiment is presented in the following bands:

• • [0091] The band 100 nearest the bottom in the layout is the simulated haplotype.
  • [0092] The second band 110 from the bottom is the haplotype molecule map for a diploid organism. These molecules (which are sorted into two haplotype classes in the layout) can be mixed and made available to the procedure of the present invention as a single set of genotype data.
  • [0093] The third band 120 from the bottom shows the results of the EM-algorithm and the set of markers that are determined to have polymorphic alleles.
  • [0094] The fourth band 130 in the layout provides the history of contig operations. From this tree, it is possible to view: 1) the developing k-neighborhoods, and 2) the distinct phased contigs.
  • [0095] The top band 140 in the layout provide the algorithmic output for this problem, including phased-in subsets that span the distance indicated by the bars above and below the loci markers. The areas where phase structure overlaps but cannot extend are regions that are of interest to target with more specific sequences, in order to extend the phasing.

[0096] Parameters of the simulations are summarized in Table 1.

[0097]

TABLE 1
Parameters* employed in performing
the simulations for Datasets 1 and 2.
Parameter	Symbol	Data Set 1	Data Set 2
Number of molecules	M	80	150
Number of fragments RFLP	F	20	100
and non RFLP
Size of the genome	G	12000	50000
Expected molecule size	EMS	2000	2000
Variance in molecule size	VMS	50	500
Variance in fragment length size	VFS	1	20
P-value that any given fragment	P-BIMODE	.5	.3
is an RFLP
Expected separation of means	ERFLPSEP	10	50
for RFLP
Variance in the separation of	VRFLPSEP	.01	6
means for RFLP
*Any parameter with both an expectation and variance can be generated with a normal distribution.

[0098] A simple VERIFY-PHASE function which checked that the posteriori distribution C_a, C_b is separated by a distance of C>0 from the point

[0099] $\left[\frac{1}{2},\frac{1}{2}\right]$
may be used as an example. In practice, it was discovered that the parameter C should preferably depend on the local coverage.

[0100] For the first simulation on dataset I, shown in FIG. 3, a relatively small set was chosen so that the potential limitations of the procedure may be revealed. Here the neighborhood size was set to k=5. False positive RFLP detections were not guarded against, yet phasings are computed. It is clear that mistakes provided therein were likely due to the low coverage library.

[0101] In the second simulation on dataset II, seen in FIG. 4, the result 200 shows that good phasing results may be achieved even on large, sparse data sets.

[0102] $\mathrm{MLE}\mathrm{Estimate}$ $\mathrm{If}d\left(\alpha ,A\right)<\in \mathrm{for}\mathrm{some}\in \mathrm{small}\mathrm{enough}.\mathrm{Maximizing}$ $\prod _{\rho \in {\left\{0,1\right\}}^{M-1}}^{}{\Theta }_{\rho }^{{\alpha }_{p}N}\mathrm{over}\Theta \in A\mathrm{is}\mathrm{equivatent}\mathrm{to}\mathrm{minimizing}$ $\sum _{j\in 1\dots M}\frac{{\left({\alpha }_j-{\Theta }_j\right)}^2}{{\alpha }_i}\mathrm{over}\Theta \in A.\mathrm{Proof}.\mathrm{Let}F\left\{\Theta \right\}=\frac{n!}{{\Pi }_{\mathrm{jw1}:k}{n}_j}{\Pi }_{\mathrm{jw1}:k}{\Theta }_j^{n_j}.\mathrm{Computing}\mathrm{the}\mathrm{second}\mathrm{variation}\text{:}$ $F^{″}\left(\theta \right)=\left(\left[\begin{array}{cccc}\frac{n_1^2}{{\theta }_1^2}& \frac{n_1{n}_2}{{\theta }_1{\theta }_2}& \cdots & \frac{n_1{n}_k}{{\theta }_1{\theta }_{\kappa }}\\ \frac{n_2{n}_1}{{\theta }_2{\theta }_1}& \frac{n_2^2}{{\theta }_2^2}& \cdots & \frac{n_2{n}_k}{{\theta }_2{\theta }_{\kappa }}\\ ⋮& ⋮& \cdots & ⋮\\ \frac{n_k{n}_1}{{\theta }_{\kappa }{\theta }_1}& \frac{n_k{n}_2}{{\theta }_{\kappa }{\theta }_2}& \cdots & \frac{n_k^2}{{\theta }_{\kappa }^2}\end{array}\right]-\left[\begin{array}{cccc}\frac{n_1}{{\theta }_1^2}& 0& \cdots & 0\\ 0& \frac{n_2}{{\theta }_2^2}& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & \frac{n_k}{{\theta }_{\kappa }^2}\end{array}\right]\right)F\left(\theta \right)$ $\mathrm{Since}F\left(\theta \right)\mathrm{is}\mathrm{smooth}\mathrm{in}\theta ,\mathrm{Taylor}\text{'}s\mathrm{remainder}\mathrm{theorem}\mathrm{gives},F\left(\Theta \right)=F\left(\alpha \right)+\nabla F\left(\alpha \right)·\left(\Theta -\alpha \right)+{\left(\Theta -\alpha \right)}^T{F}^{″}\left(\alpha \right)\left(\Theta -\alpha \right)+o\left({\left(\Theta -\alpha \right)}_2\right)$ $\mathrm{When}\alpha =\frac{n_1}{n},\dots ,\frac{n_k}{n},\nabla F\left(\alpha \right)=0,\mathrm{this}\mathrm{is}a\mathrm{standard}\mathrm{MLE}$ $\mathrm{result}\mathrm{for}a\mathrm{multi}\text{-}\mathrm{nominal}\mathrm{distribution}.\mathrm{Computing}\mathrm{the}\mathrm{quadratic}\mathrm{function}\text{:}$ $\begin{array}{c}{\left(\Theta -\alpha \right)}^T{F}^{″}\left(\alpha \right)\left(\Theta -\alpha \right)={\left(\Theta -\alpha \right)}^T\\ \left(n^2\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& 1& \cdots & 1\\ ⋮& ⋮& \cdots & ⋮\\ 1& 1& \cdots & 1\end{array}\right]-n\left[\begin{array}{cccc}\frac{1}{{\alpha }_1}& 0& \cdots & 0\\ 0& \frac{1}{{\alpha }_2}& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & \frac{1}{{\alpha }_{\kappa }}\end{array}\right]\right)\\ F\left(\alpha \right)·\left(\Theta -\alpha \right)\\ =F\left(\alpha \right){n^2\left(\sum _j\left({\Theta }_j-{\alpha }_j\right)\right)}^2-F\left(\alpha \right)n\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}\\ =-F\left(\alpha \right)n\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}\end{array}$ $\mathrm{Thus}\mathrm{for}\mathrm{an}\theta \mathrm{very}\mathrm{near}\mathrm{to}\alpha \mathrm{the}\mathrm{level}\mathrm{curves}\mathrm{of}F\mathrm{are}$ $\mathrm{given}\mathrm{by}{\Theta }_{\delta }=\left\{\theta \text{:}F\left(\theta \right)=F\left(\alpha \right)-\delta \right\}\mathrm{are}\mathrm{approximately}\mathrm{ellipsoids}.\text{}F\left(\Theta \right)=F\left(\alpha \right)-F\left(\alpha \right)n\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}+o\left({\left(\Theta -\alpha \right)}_2\right)=F\left(\alpha \right)-F\left(\alpha \right)n\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}+o\left(\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}\right)$ $\mathrm{Let}{\left(\Theta -\alpha \right)}_{\alpha }^2=\sum _j\frac{{\left({\Theta }_j-{\alpha }_j\right)}^2}{{\alpha }_j}.\mathrm{Letting}{L}_1=L\left(\alpha \right)\mathrm{and}$ $\mathrm{assuming}\mathrm{there}\mathrm{is}a\mathrm{second}\mathrm{local}\mathrm{optima}\mathrm{for}\mathrm{the}\mathrm{likelihood}$ $\mathrm{function}\mathrm{value}\mathrm{at}{L}_2,\mathrm{let}V\left(\alpha \right)=\left\{\Theta \text{:}L\left(\Theta \right)>L_1-\frac{L_1-L_2}{2}\right\}.\mathrm{We}\mathrm{must}\mathrm{show}\mathrm{that}\mathrm{there}\mathrm{is}a\delta \mathrm{so}\mathrm{that}$ $\left\{\Theta \text{:}{\left(\Theta -\alpha \right)}_{\alpha }^2<\delta \right\}⋐V\left(\alpha \right).\mathrm{And}\mathrm{this}\mathrm{is}\mathrm{clear}\mathrm{from}\mathrm{the}\mathrm{inequality}$ $L_1-L_{1}n{\Theta -\alpha }_{\alpha }^2-o\left({\Theta -\alpha }_{\alpha }^2\right)<F\left(\Theta \right)<L_1-L_{1}n{\Theta -\alpha }_{\alpha }^2+o\left({\Theta -\alpha }_{\alpha }^2\right)$ $\mathrm{by}\mathrm{choosing}a\delta \mathrm{small}\mathrm{enough}\mathrm{that}{L}_{1}n\delta +o\left(\delta \right)<\frac{L_1-L_2}{2}.\mathrm{We}\mathrm{conclude}$ $\mathrm{that}\mathrm{if}\mathrm{there}\mathrm{is}a\mathrm{point}\mathrm{of}A\in \left\{\Theta :{\Theta -\alpha }_{\alpha }^2<\delta \right\}\mathrm{then}\mathrm{it}\mathrm{must}\mathrm{be}\mathrm{the}$ $\mathrm{unique}\mathrm{maxima}\mathrm{in}A\mathrm{for}\mathrm{our}\mathrm{likelihood}\mathrm{function}.$

[0103] The Chi-squared statistical test for determining whether allelic data at loci j and j′ display linkage disequilibrium, and hence are not in Hardy-Weinberg Equilibrium (HWE) has been reviewed in great detail. See Weir, 1996, Genetic Data Analysis II, Sinauer Associates, Sunderland, Mass. for details and a complete statistical treatment. The Chi-squared statistical test for gametic disequilibrium at two loci has been modified by using additive disequilibrium coefficients to adjust to our population model. The end result is a Chi-squared statistical test that allows the rejection of HWE from observed frequencies alone. Since determination of linkage is a prerequisite to phasing, or at least in finding structure in the joint distribution over allele spaces of adjacent loci, this statistical test is important. The boundaries of haplotype blocks (or phased contigs) are an interesting and important issue in understanding population dynamics.

[0104] For example, let D_ab denote the disequilibrium coefficient between alleles a at loci j and b at loci ^t:
D_ab=p_ab−p_ap_b

• • [0105] Where p_ab, p_a, p_b are the population frequencies for allele type: ab, a, b respectively. In the presence of HWE D_ab can be expected to be zero. Letting {circumflex over (D)}_ab denote an estimate from estimate frequencies:
    {circumflex over (D)}_ab={tilde over (p)}_ab−{tilde over (p)}_a{tilde over (p)}_b
  • [0106] with:
    p_a=1/N Σ_i=1:Nα_aj(d_ij), p_b=1/N Σ_i=1:N α_bj′(d_ij′)

[0107] Computation of Expectation and Variance:

[0108] $E\left({\hat{D}}_{\mathrm{ab}}\right)=\frac{N-1}{N}{D}_{\mathrm{ab}}$ $V\left({\hat{D}}_{\mathrm{ab}}\right)\approx \frac{1}{N}\left[p_a{q}_a{p}_b{q}_b+\left(1-2p_a\right)\left(1-2q_a\right)D_{\mathrm{ab}}-D_{\mathrm{ab}}^2\right]$

[0109] The variance can be computed using Fisher's approximate variance formula. Under the assumption that loci j and j′ are in HWE, D_ab=0 and:

[0110] $E_{\mathrm{HWE}}^{\prime }\left({\hat{D}}_{\mathrm{ab}}\right)=0$ $V_{\mathrm{HWE}}\left({\hat{D}}_{\mathrm{ab}}\right)=\frac{1}{N}p_a{q}_a{p}_b{q}_b$

[0111] From this information, it is possible to construct a Chi-Squared test to evaluate the hypothesis that alleles a and b at loci j and j′ are acting as they would if they were in HWE.

[0112] ${\chi }_{\mathrm{ab}}^2=z^2=\frac{{\mathrm{ND}}_{\mathrm{ab}}^2}{{\tilde{p}}_a^A{\tilde{q}}_a^A{\tilde{p}}_b^B{\tilde{q}}_b^B}$

[0113] It is possible to reject the HWE hypothesis correctly 9 times in 10 by using a reference value of z²>2.71, or we may reject HWE correctly 99 times in 100 using reference values z²>6.63. If alleles are linked by a haplotype, this test may be used as the VERIFY-PHASE function mentioned previously in this text.

[0114] EM-Algorithm Analytic Results

[0115] A. An Example Using RFLP Markers

[0116] The data at loci j can refer to the observed distances between restriction sites j and j+1, as they are derived from two haplotypes H₁ and H₂ with underlying genome distances μ₁ and μ₂. The distribution of data points for loci j is given by:

[0117] $f_j\left(x\right)=\frac{1}{\sqrt{2{\mathrm{\pi \mu }}_{\mathrm{j1}}^2}}\exp \left(\frac{-{\left(x-{\mu }_{\mathrm{j1}}\right)}^2}{2{\mu }_{\mathrm{j1}}^2}\right){\alpha }_{\mathrm{j1}}\left(x\right)+\frac{1}{\sqrt{2{\mathrm{\pi \mu }}_{\mathrm{j2}}^2}}\exp \left(\frac{-{\left(x-{\mu }_{\mathrm{j2}}\right)}^2}{2{\mu }_{\mathrm{j2}}^2}\right){\alpha }_{\mathrm{j2}}\left(x\right)$

[0118] It is preferable to make a simplifying assumption that σ=½ (μj₁+μj₂) so that f_j may be closely approximated by:

[0119] $F_j\left(x\right)=\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}^2}}\exp \left(\frac{-{\left(x-{\mu }_{j1}\right)}^2}{2{\sigma }^2}\right){\alpha }_{j1}\left(x\right)+\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}^2}}\exp \left(\frac{-{\left(x-{\mu }_{j2}\right)}^2}{2{\sigma }^2}\right){\alpha }_{j2}\left(x\right)$

[0120] For loci j the set of points {d_ij=i ε [1, 2, . . . N]} is data. It is preferable to infer the model parameters Φ={σ=μ₁, μ₂} and posteriori distribution or by use of the EM-algorithm. The subscript j can be dropped in the following equation, the objective being to iteratively optimize the function:

[0121] $H\left(\alpha ,\Phi \right)=\sum _{i\in 1:n}\sum _{k\in 1:2}\left({\alpha }_k\left(d_{\mathrm{ji}}\right)\ln G_k\left(d_{\mathrm{ji}}❘\Phi \right)-{\sigma }_k\left(d_{\mathrm{ji}}\right)\ln \left({\alpha }_k\left(d_{\mathrm{ji}}\right)\right)\right)$ $\mathrm{With}{G}_k\left(x❘\varphi \right)=\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}^2}}\exp \left(\frac{-{\left(x-{\mu }_k\right)}^2}{2{\sigma }^2}\right)\mathrm{the}k\mathrm{th}\mathrm{Gaussian}$ $\mathrm{kernel}.$

[0122] Optimization can be done in two steps:

• • [0123] 1. E-STEP Holding Φ fixed, optimize H(α, Φ) over α, letting Φ be the previous estimate of parameters.
  • [0124] The result for the argmax α is:

[0125] ${\alpha }_k\left(x\right)\leftarrow \frac{G_k\left(x❘\Phi \right)}{\sum _i{G}_l\left(x❘\Phi \right)}$

• • [0126] 2. M-STEP
  • [0127] Holding α fixed, optimize H(αΦ) over Φ (using the previous estimate of parameters on the hidden categories α (y), which depend on the previous estimate denoted {circumflex over (Φ)}=[{circumflex over (μ)}₁, {circumflex over (μ)}₂, ]). The result argmin H(&, Φ) is:

[0128] ${\mu }_k\leftarrow \frac{\sum _{i=1:n}{\delta }_k\left(d_{\mathrm{ij}}\right)d_{\mathrm{ij}}}{\sum _{i=1:n}{\delta }_k\left(d_{\mathrm{ij}}\right)}$ $\sigma \leftarrow \sqrt{\frac{1}{2}\sum _{j\ne 1:2}\frac{1}{N}\sum _{i=1:N}{\hat{\alpha }}_j\left(d_{\mathrm{ij}}\right){\left(d_{\mathrm{ij}}-{\hat{\mu }}_j\right)}^2}$

[0129] The EM-algorithm procedure can be executed until convergence in the parameter space occurs. Exemplary detailed computations for the E-Step and M-Step are provided in the appendix. Detailed proofs of each step are provided below.

[0130] E-Step

[0131] Proof. Consider the calculus problem of optimizing:

[0132] $f\left(\varphi \right)=\varphi \left(A_1-\log \left(B\varphi \right)\right)+\left(1-\varphi \right)\left(A_2-\log \left(B\left(1-\varphi \right)\right)\right)$ $f^{\prime }\left(\varphi \right)=0\Rightarrow {\varphi }^{\prime }=\left(\frac{1}{e^{A_1-A_2}+1}\right)$

[0133] Φε (0, 1). Apply this fact to the optimization problem of finding numbers Q=<α_iv>_1:N,vw1:2, so that the following function is optimized.

[0134] $\sum _{i=1:n}\sum _{\nu =1:2}\left({\alpha }_{i\nu }\left(A_{i\nu }-\log \left(B_{{\alpha }_{\omega }}\right)\right)\right)=\sum _{i=i:n}{\alpha }_{i\nu }\left(A_{i1}-\log \left(B_{{\alpha }_{\omega }}\right)\right)+\left(1-{\alpha }_{\omega }\right)\left(A_{i2}-\log \left(B\left(1-{\alpha }_{\omega }\right)\right)\right)$

• • [0135] Where

[0136] $A_1=\frac{{\left(\mathrm{dj}-\mu 1\right)}^2}{2\sigma 2}\mathrm{and}{A}_2=\frac{{\left(\mathrm{dj}-\mu 2\right)}^2}{2\sigma 2}\mathrm{and}B=\sqrt{2{\mathrm{\pi \sigma }}^2}.$

• • [0137]  It is shown that the answer is given by maximizing each summand and hence given by:

[0138] $\begin{array}{c}{\alpha }_1=\left(\frac{1}{e^{\left(\frac{{\left(d_i-{\mu }_1\right)}^2}{2{\sigma }^2}-\frac{{\left(d_i-{\mu }_2\right)}^2}{2{\sigma }^2}\right)}+1}\right)\\ =\frac{G_1\left(d_i\right)}{G_1\left(d_i\right)+G_2\left(d_i\right)}\end{array}$

• • [0139]  and similarly for σ₂.

[0140] M-Step

[0141] Proof. Consider the calculus problem of optimizing:

[0142] $\begin{array}{c}f\left({\mu }_1,{\mu }_2,\sigma \right)=\sum _{i=1:N}\sum _{\nu =1:2}{q}_{\nu i}\log \left(\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}^2}}\exp \frac{-{\left({\alpha }_i-{\mu }_{\nu }\right)}^2}{2{\sigma }^2}\right)+H\\ =\sum _{i=1:N}\sum _{\nu =1:2}{q}_{\nu i}\left(\frac{-{\left({\alpha }_i-{\mu }_{\nu }\right)}^2}{2{\sigma }^2}-\log \sqrt{2{\mathrm{\pi \sigma }}^2}\right)+H\end{array}$

• • [0143] Where H is constant in (μ₁, μ₂, σ). Consider the partial derivative of f with respect o μ_v:

[0144] $\begin{array}{c}\frac{\partial f}{\partial {\mu }_{\nu }}=\sum _{i=1:N}\frac{\partial }{\partial {\mu }_{\nu }}\left(\frac{q_{\nu i}}{2{\sigma }^2}\left(-{\alpha }_i^2+2{\alpha }_i{\mu }_{\nu }-{\mu }_{\nu }^2-2{\sigma }^2\log \sqrt{2{\mathrm{\pi \sigma }}^2}\right)\right)\\ =\frac{1}{{\sigma }^2}\sum _{i=1:N}{q}_{\nu i}{\alpha }_i-{\mu }_{\nu }{q}_{\nu i}\end{array}$

• • • [0145] Thus, the following is obtained:

[0146] $\frac{\partial f}{\partial {\mu }_{\nu }}=0\Leftarrow {\mu }_{\nu }=\frac{\sum _{i=1:N}{q}_{i\nu }{\alpha }_i}{\sum _{i=1:N}{q}_{i\nu }}$

• • • [0147] Now consider the partial off with respect to σ:

[0148] $\frac{\partial f}{\partial \sigma }=\sum _{i=1:N}\sum _{\nu =1:2}{q}_{\nu i}\left(\frac{{\left({\alpha }_i-{\mu }_{\nu }\right)}^2}{2{\sigma }^2}-1\right)$

[0149] Section 4. Conclusions

[0150] The simulation results described herein demonstrate that locally the phasing may be highly accurate. When local coverage derived from one haplotype is low, the detection of polymorphisms can become difficult. In the first dataset, a false negative detection may be found on the 8th marker from the left. This false negative was due to zero coverage from one of the haplotypes at that point. The ninth marker can be a false positive detection, and may be attributed to zero coverage from one haplotype and low coverage (two molecules) from the alternative haplotype. The false positive does not cause errors in the phase information for correctly detected polymorphic loci in the phased-contig achieved over marker index in the set {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19}. Designing a mapping experiment targeting a polymorphic marker in the set {6, 7, 8, 9, 10} can allow one to phase the two contigs into a single contig.

[0151] The results of the method for haplotype phasing of the present invention may be assessed in terms of two absolute quantities: resolution and accuracy. Not only does the method of the present invention provide conclusions, it also reports a confidence level associated with the conclusion. Since an individual's haplotype structure is singular and absolute, any method that makes conclusions about that structure should assess its own accuracy. The method according to the present invention provides this assessment in such a way that the resolution and accuracy of the conclusions drawn for an individual scale with the amount of effort (i.e. mapping experiments) expended on that individual. This feature of the method is important because (a) an accurate and high resolution determination of an individual's haplotype structure may be drawn in the absence of knowledge or information from her/his pedigree or sub-population, in contrast to population-based studies, which rely on mapping data from a set of closely related individuals, and (b) since resolution and accuracy scale with effort, the same scaling relationship is present for cost. In various applications of the method and system of the instant invention where different resolutions and accuracies are required, the present method provides a better estimate for the same cost when compared to currently available alternative methods.

[0152] FIG. 5 shows a block diagram of an exemplary embodiment of a system 300 according to the present invention which is capable of storing thereon the storage arrangement of the present invention, and operable to execute thereon the method according to the present invention. In particular, the system includes a processing arrangement 310 and a storage arrangement 320. The storage arrangement 320 can be one or more hard drives, memory (read-only memory, random access memory—“RAM”, DRAM), CD-ROMs, floppy disks, etc. and/or combination thereof, and may store thereon a software arrangement (e.g., a software program). The software program can be accessed by the processing arrangement 310 (e.g., a processor such as Intel Pentium® processor), and the software arrangement can make the processing arrangement operable to execute the exemplary embodiment of the present described herein.

[0153] For example, the system 300 can receive genotype data associated with one or more genetic sample, either from external sources or from the storage arrangement 320. The processing arrangement 310 (executing the software arrangement obtained from the storage arrangement 320) obtains the data, and performs the procedures in accordance with the exemplary embodiment of the present invention. After the processing is completed by the processing arrangement 310, the processing arrangement can obtain and output haplotype data for respective genetic sample (in block 340). In addition or as an alternative, the processing arrangement 310 can obtain a confidence level of accuracy of the obtained haplotype data (block 350), and/or results of examination of a linear number of polymorphic genetic markers (block 360).

[0154] FIG. 6 shows a flow diagram of a top level of one exemplary embodiment of the method 400 according to the present invention for establishing the haplotype of one or more genetic samples from genotype data obtained from the corresponding one or more genetic samples. This method can be executed by the system 300 of FIG. 5, or by any other arrangement or system that is capable of implementing the method. For example, in step 410, the genotype data (associated with one or more genetic samples) is obtained. Then, in step 420, the polymorphic genetic markers are identified from the genotype data. Further, one or more associations of the polymorphic genetic markers across genetic loci are determined to obtain the information associated with the haplotype data.

[0155] While the invention has been described in connecting with preferred embodiments, it will be understood by those of ordinary skill in the art that other variations and modifications of the preferred embodiments described above may be made without departing from the scope of the invention. Other embodiments will be apparent to those of ordinary skill in the art from a consideration of the specification or practice of the invention disclosed herein. For example, the exemplary embodiments of the present invention can also be applicable for polyploid organisms, as well as diploid organisms. It should be understood that the present invention is operable on one or more chromosome of an organism, and can be used even if certain marker information is ambiguous or missing altogether. It is intended that the specification and the described examples are considered as exemplary only, with the true scope and spirit of the invention indicated by the following claims. Additionally, all references cited herein are hereby incorporated by this reference as though set forth fully herein.